lcsim/src/org/lcsim/contrib/Cassell/recon/Cheat
diff -u -r1.3 -r1.4
--- CheatReconDriver.java 13 Sep 2007 14:41:29 -0000 1.3
+++ CheatReconDriver.java 16 Oct 2007 18:19:08 -0000 1.4
@@ -34,18 +34,26 @@
import java.util.Map;
import org.lcsim.util.lcio.LCIOConstants;
import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.BasicHep3Vector;
+import org.lcsim.util.aida.AIDA;
/**
*
* @author cassell
*/
public class CheatReconDriver extends Driver
{
+ private AIDA aida = AIDA.defaultInstance();
String inFSname = "GenFinalStateParticles";
String inTname = "RefinedCheatTracks";
String inCname = "RefinedCheatClusters";
String outRname;
+ String outVRname;
+ String outGenRname;
+ String outGenVRname;
String outFSname;
+ String outTrFSname;
String outTname;
+ String outFSTname;
String outCname;
String[] diginames = {"EcalBarrDigiHits", "EcalEndcapDigiHits","HcalBarrDigiHits", "HcalEndcapDigiHits"};
//
@@ -58,20 +66,21 @@
// Ask Jeremy how to do this
double emR = 1270.;
double emZ = 1680.;
-// Just keep it distance cut
- double jki = 100.;
// Too close cut
- double tccut = 400.;
-// Minimum length of neutral parent to consider for replacement
- double minneutlen = 400.;
+ double tccut = 200.;
Map<MCParticle,Track> tmap;
TraceOrigin to;
/** Creates a new instance of CheatSimReconParticleDriver */
public CheatReconDriver()
{
outRname = "ReconPerfectReconParticles";
+ outVRname = "ReconPerfectVisReconParticles";
+ outGenRname = "GenPerfectReconParticles";
+ outGenVRname = "GenPerfectVisReconParticles";
outFSname = "ReconFSParticles";
+ outTrFSname = "TrackedReconFSParticles";
outTname = "ReconTracks";
+ outFSTname = "FSReconTracks";
outCname = "ReconClusters";
//
// Make an initial list of MCParticles to be considered
@@ -96,29 +105,36 @@
{outCname = s;}
public void setCheatFSParticleOutputName(String s)
{outFSname = s;}
+ public void setCheatTrackedFSParticleOutputName(String s)
+ {outTrFSname = s;}
+ public void setCheatFSTrackOutputName(String s)
+ {outFSTname = s;}
public void setCheatReconstructedParticleOutputName(String s)
{outRname = s;}
- public void setJustKeepItCut(double x)
- {jki = x;}
+ public void setCheatVisReconstructedParticleOutputName(String s)
+ {outVRname = s;}
+ public void setGenReconstructedParticleOutputName(String s)
+ {outGenRname = s;}
+ public void setGenVisReconstructedParticleOutputName(String s)
+ {outGenVRname = s;}
public void setTooCloseCut(double x)
{tccut = x;}
- public void setMinNeutralLengthCut(double x)
- {minneutlen = x;}
public void process(EventHeader event)
{
super.process(event);
List<MCParticle> initfs = event.get(MCParticle.class,inFSname);
List<Track> inittracks = event.get(Track.class,inTname);
List<Cluster> initclusters = event.get(Cluster.class,inCname);
- List<MCParticle> add = new ArrayList<MCParticle>();
- List<MCParticle> remove = new ArrayList<MCParticle>();
- List<MCParticle> tempfs = new ArrayList<MCParticle>();
- List<MCParticle> untracked = new ArrayList<MCParticle>();
List<MCParticle> reconfs = new ArrayList<MCParticle>();
+ List<MCParticle> trreconfs = new ArrayList<MCParticle>();
List<Track> fstracks = new ArrayList<Track>();
+ List<Track> tracks = new ArrayList<Track>();
List<Cluster> fsclusters = new ArrayList<Cluster>();
List<Cluster> nonfsclusters = new ArrayList<Cluster>();
- List<Cluster> untrackedchargedclusters = new ArrayList<Cluster>();
+ List<ReconstructedParticle> perfect = new ArrayList<ReconstructedParticle>();
+ List<ReconstructedParticle> visperfect = new ArrayList<ReconstructedParticle>();
+ List<ReconstructedParticle> genperfect = new ArrayList<ReconstructedParticle>();
+ List<ReconstructedParticle> genvisperfect = new ArrayList<ReconstructedParticle>();
tmap = new HashMap<MCParticle,Track>();
//
// Make a map of trackable tracks
@@ -126,262 +142,140 @@
for(Track t:inittracks)
{
MCParticle p = ((BaseTrackMC) t).getMCParticle();
- if(isTrackable(p,t))tmap.put(p,t);
+ if(isTrackable(p,t))
+ {
+ tmap.put(p,t);
+ tracks.add(t);
+ }
}
//
// Start with the generator final state particles
//
+ double Einit = 0.;
+ double Eadd = 0.;
+ double Esub = 0.;
for(MCParticle p:initfs)
{
-//
-// If the particle has already been removed skip it
-//
- if(remove.contains(p))continue;
+ Einit += p.getEnergy();
//
// For the majority of generator fs particles, they originate near the ip
// and end near or beyond the Ecal. Just keep them
//
Hep3Vector vtx = p.getOrigin();
Hep3Vector ep = p.getEndPoint();
- if(vtx.magnitude() < jki)
+ if(vtx.magnitude() < tccut)
{
- if( (Math.sqrt(ep.x()*ep.x()+ep.y()*ep.y())>emR-jki)||(Math.abs(ep.z())>emZ-jki) )
+ if( (Math.sqrt(ep.x()*ep.x()+ep.y()*ep.y())>emR-tccut)||(Math.abs(ep.z())>emZ-tccut) )
{
+ reconfs.add(p);
continue;
}
}
//
-// Check if we really want the parent
-//
- MCParticle tp = getTrackableParent(p);
-//
-// No trackable parentage. Only keep parent if really close to Cal
+// Check if we really want the parent: Only if the vtx is close to Ecal
//
- if(tp == null)
+ if( (Math.sqrt(vtx.x()*vtx.x()+vtx.y()*vtx.y())>emR-tccut)||(Math.abs(vtx.z())>emZ-tccut) )
{
-//
-// If particle is trackable, we aren't going to remove it for parentage
-//
- if(!tmap.containsKey(p))
+ MCParticle par = p.getParents().get(0);
+ Hep3Vector pvtx = par.getOrigin();
+ while( (Math.sqrt(pvtx.x()*pvtx.x()+pvtx.y()*pvtx.y())>emR-tccut)||(Math.abs(pvtx.z())>emZ-tccut) )
{
- if(tooClose(vtx))
- {
-//
-// Vertex is close to Ecal. Consider replacing with a neutral parent
-//
- MCParticle par = getLongestNeutralParent(p);
- if(par != null)
- {
- if(!hasTrackableDaughters(par))
- {
- add.add(par);
- List<MCParticle> dl = new ArrayList<MCParticle>();
- for(MCParticle d:par.getDaughters())
- {
- dl.add(d);
- }
- while(dl.size() > 0)
- {
- remove.add(dl.get(0));
- for(MCParticle d:dl.get(0).getDaughters())
- {
- dl.add(d);
- }
- dl.remove(0);
- }
- continue;
- }
- }
- }
+ par = par.getParents().get(0);
+ pvtx = par.getOrigin();
}
- }
-//
-// Have trackable parentage. Only keep parent if no trackable offspring
-// of trackable parent.
-//
- else
- {
- if(!hasTrackableDaughters(tp))
+// System.out.println("Parent replacement: vtx close to Ecal");
+ if(!reconfs.contains(p))
{
-//
-// No trackable offspring, so keep trackable parent and remove all
-// offspring.
-//
- add.add(tp);
- List<MCParticle> dl = new ArrayList<MCParticle>();
- for(MCParticle d:tp.getDaughters())
- {
- dl.add(d);
- }
- while(dl.size() > 0)
- {
- remove.add(dl.get(0));
- for(MCParticle d:dl.get(0).getDaughters())
- {
- dl.add(d);
- }
- dl.remove(0);
- }
- continue;
+ reconfs.add(par);
+ Eadd += par.getEnergy();
+// System.out.println("Adding to FS list "+par.getType().getName()+": E = "+par.getEnergy());
}
+ Esub += p.getEnergy();
+// System.out.println("Subtracting from FS list "+p.getType().getName()+": E = "+p.getEnergy());
+ continue;
}
//
// Done checking for parentage replacement. Now check for offspring replacement
//
//
-// Start with a special case for photon conversions. Only replace the photon
-// if both electrons are trackable
+// If the endpoint is close to or beyond the Ecal, keep it
//
- if(p.getPDGID() == 22)
+ if( (Math.sqrt(ep.x()*ep.x()+ep.y()*ep.y())>emR-tccut)||(Math.abs(ep.z())>emZ-tccut) )
{
- if(hasTrackableDaughters(p))
- {
- boolean hasep = false;
- boolean hasem = false;
- for(MCParticle d:p.getDaughters())
- {
- if(d.getSimulatorStatus().isBackscatter())continue;
- if(d.getSimulatorStatus().vertexIsNotEndpointOfParent())continue;
- if(tmap.containsKey(d))
- {
- if(d.getPDGID() == 11)hasem = true;
- if(d.getPDGID() == -11)hasep = true;
- }
- }
- if(hasem&&hasep)
- {
- remove.add(p);
- for(MCParticle d:p.getDaughters())
- {
- if(d.getSimulatorStatus().isBackscatter())continue;
- if(d.getSimulatorStatus().vertexIsNotEndpointOfParent())continue;
- add.add(d);
- }
- }
- }
+ reconfs.add(p);
continue;
}
//
-// And on to the general case
+// If the particle is trackable and there are no trackable offspring, keep it.
+// Otherwise, replace it.
//
- if(hasTrackableDaughters(p))
- {
- remove.add(p);
- for(MCParticle d:p.getDaughters())
- {
- if(d.getSimulatorStatus().isBackscatter())continue;
- if(d.getSimulatorStatus().vertexIsNotEndpointOfParent())continue;
- add.add(d);
- }
- }
- else
+
+ if(tmap.containsKey(p)&&!hasTrackableDaughters(p))
{
- if(tooClose(ep))continue;
- if(tmap.containsKey(p))continue;
- remove.add(p);
- for(MCParticle d:p.getDaughters())
- {
- if(d.getSimulatorStatus().isBackscatter())continue;
- if(d.getSimulatorStatus().vertexIsNotEndpointOfParent())continue;
- add.add(d);
- }
+ reconfs.add(p);
+ continue;
}
- }
//
-// Done with pass over original FS particles. Now need to check added particles
-// for further offspring replacement
+// Now do the offspring replacement. Only keep trackable particles and
+// particles which reach the Ecal.
//
- List<MCParticle> newadd = new ArrayList<MCParticle>(add);
- while(newadd.size() > 0)
- {
- List<MCParticle> add2 = new ArrayList<MCParticle>();
- for(MCParticle p:newadd)
- {
- if(p.getPDGID() == 22)
- {
- if(hasTrackableDaughters(p))
- {
- boolean hasep = false;
- boolean hasem = false;
- for(MCParticle d:p.getDaughters())
- {
- if(d.getSimulatorStatus().isBackscatter())continue;
- if(d.getSimulatorStatus().vertexIsNotEndpointOfParent())continue;
- if(tmap.containsKey(d))
- {
- if(d.getPDGID() == 11)hasem = true;
- if(d.getPDGID() == -11)hasep = true;
- }
- }
- if(hasem&&hasep)
- {
- for(MCParticle d:p.getDaughters())
- {
- if(d.getSimulatorStatus().isBackscatter())continue;
- if(d.getSimulatorStatus().vertexIsNotEndpointOfParent())continue;
- add2.add(d);
- }
- continue;
- }
- }
- tempfs.add(p);
- continue;
- }
- if(hasTrackableDaughters(p))
- {
- for(MCParticle d:p.getDaughters())
- {
- if(d.getSimulatorStatus().isBackscatter())continue;
- if(d.getSimulatorStatus().vertexIsNotEndpointOfParent())continue;
- add2.add(d);
- }
+ Esub += p.getEnergy();
+// System.out.println("Offspring replacement");
+// System.out.println("Removing "+p.getType().getName()+": E = "+p.getEnergy());
+ List<MCParticle> repl = new ArrayList<MCParticle>(p.getDaughters());
+ while(repl.size() > 0)
+ {
+ MCParticle d = repl.get(0);
+ repl.remove(d);
+ if(d.getSimulatorStatus().isBackscatter())continue;
+ if(d.getSimulatorStatus().vertexIsNotEndpointOfParent())continue;
+ if(tmap.containsKey(d))
+ {
+ reconfs.add(d);
+ Eadd += d.getEnergy();
+// System.out.println("Adding "+d.getType().getName()+": E = "+d.getEnergy());
}
else
{
- if(tooClose(p.getEndPoint()))
- {
- tempfs.add(p);
- continue;
- }
- if(tmap.containsKey(p))
+ Hep3Vector dep = d.getEndPoint();
+ if( (Math.sqrt(dep.x()*dep.x()+dep.y()*dep.y())>emR)||(Math.abs(dep.z())>emZ) )
{
- tempfs.add(p);
- continue;
+ reconfs.add(d);
+ Eadd += d.getEnergy();
+// System.out.println("Adding "+d.getType().getName()+": E = "+d.getEnergy());
}
- for(MCParticle d:p.getDaughters())
+ else
{
- if(d.getSimulatorStatus().isBackscatter())continue;
- if(d.getSimulatorStatus().vertexIsNotEndpointOfParent())continue;
- add2.add(d);
+ repl.addAll(d.getDaughters());
}
}
}
- newadd = new ArrayList<MCParticle>(add2);
}
+// System.out.println("Done with replacement lists:");
+// System.out.println("E GenFS = "+Einit+": Added E = "+Eadd+": Subtracted E = "+Esub);
//
-// Now have remove and add lists for a temp fs particle list. Make it.
+// Now make the recon lists.
//
- for(MCParticle p:initfs)
- {
- if(!remove.contains(p))tempfs.add(p);
- }
- to = new TraceOrigin(tempfs);
//
-// Now make the recon lists. (Removes untrackable tracks, and invisible particles
+// First the tracks. All trackable and final state trackable
//
- Map<MCParticle,Cluster> fsclmap = new HashMap<MCParticle,Cluster>();
- for(Track t:inittracks)
+ for(Track t:tracks)
{
MCParticle p = ((BaseTrackMC) t).getMCParticle();
- if(tmap.containsKey(p))
+ if(reconfs.contains(p))
{
- if(tempfs.contains(p))
- {
- reconfs.add(p);
- fstracks.add(t);
- }
+ fstracks.add(t);
+ trreconfs.add(p);
}
}
+//
+// Now the clusters. Trace the cheated clusters back to the reconfs
+// list. Also trace to genfs list.
+//
+ to = new TraceOrigin(reconfs);
+ Map<MCParticle,Cluster> fsclmap = new HashMap<MCParticle,Cluster>();
+ TraceOrigin gento = new TraceOrigin(initfs);
+ Map<MCParticle,Cluster> genfsclmap = new HashMap<MCParticle,Cluster>();
for(Cluster c:initclusters)
{
CheatCluster cc = (CheatCluster) c;
@@ -402,43 +296,77 @@
CheatCluster bc = new CheatCluster(trmc);
bc.addCluster(c);
fsclmap.put(trmc, bc);
- if(!reconfs.contains(trmc))reconfs.add(trmc);
}
}
- }
- List<MCParticle> rm = new ArrayList<MCParticle>();
- for(MCParticle p:reconfs)
- {
- if(p.getCharge() != 0)
+ trmc = gento.traceit(cc.getMCParticle());
+ if(trmc != null)
{
- if(!tmap.containsKey(p))
+ if(genfsclmap.containsKey(trmc))
{
- rm.add(p);
- untrackedchargedclusters.add(fsclmap.get(p));
- fsclmap.remove(p);
+ CheatCluster bc = (CheatCluster) genfsclmap.get(trmc);
+ bc.addCluster(c);
+ }
+ else
+ {
+ CheatCluster bc = new CheatCluster(trmc);
+ bc.addCluster(c);
+ genfsclmap.put(trmc, bc);
}
}
}
- for(MCParticle p:rm)
- {
- reconfs.remove(p);
- }
+//
+// Have the map, make the list of clusters
+//
+ fsclusters.addAll(fsclmap.values());
+//
+// Output the lists to event
+//
event.put(outFSname,reconfs,MCParticle.class,0);
+ event.put(outTrFSname,trreconfs,MCParticle.class,0);
int flag = 1<<LCIOConstants.TRBIT_HITS;
- event.put(outTname,fstracks,BaseTrackMC.class,flag);
+ event.put(outTname,tracks,BaseTrackMC.class,flag);
+ event.put(outFSTname,fstracks,BaseTrackMC.class,flag);
flag = 1<<LCIOConstants.CLBIT_HITS;
- event.put(outCname,new ArrayList(fsclmap.values()),CheatCluster.class,flag);
+ event.put(outCname,fsclusters,CheatCluster.class,flag);
event.put("NonFSReconClusters",nonfsclusters,CheatCluster.class,flag);
- event.put("UntrackedChargedClusters",untrackedchargedclusters,CheatCluster.class,flag);
- List<ReconstructedParticle> rrl = new ArrayList<ReconstructedParticle>();
for(MCParticle p:reconfs)
{
+//
+// Eliminate the neutrinos
+//
+ int pid = Math.abs(p.getPDGID());
+ if( (pid == 12)||(pid == 14)||(pid == 16) )continue;
MCReconstructedParticle rp = new MCReconstructedParticle(p);
if(tmap.containsKey(p))rp.addTrack(tmap.get(p));
if(fsclmap.containsKey(p))rp.addCluster(fsclmap.get(p));
- rrl.add(rp);
+ perfect.add(rp);
+ if(tmap.containsKey(p)||fsclmap.containsKey(p))visperfect.add(rp);
+ }
+ event.put(outRname,perfect);
+ event.put(outVRname,visperfect);
+ for(MCParticle p:initfs)
+ {
+//
+// Eliminate the neutrinos
+//
+ int pid = Math.abs(p.getPDGID());
+ if( (pid == 12)||(pid == 14)||(pid == 16) )continue;
+ MCReconstructedParticle rp = new MCReconstructedParticle(p);
+ if(tmap.containsKey(p))rp.addTrack(tmap.get(p));
+ List<MCParticle> dlist = new ArrayList<MCParticle>(p.getDaughters());
+ while(dlist.size() > 0)
+ {
+ MCParticle d = dlist.get(0);
+ dlist.remove(d);
+ if(tmap.containsKey(d))rp.addTrack(tmap.get(d));
+ dlist.addAll(d.getDaughters());
+ }
+ if(genfsclmap.containsKey(p))rp.addCluster(genfsclmap.get(p));
+ genperfect.add(rp);
+ if(tmap.containsKey(p)||genfsclmap.containsKey(p))genvisperfect.add(rp);
}
- event.put(outRname,rrl);
+ event.put(outGenRname,genperfect);
+ event.put(outGenVRname,genvisperfect);
}
private boolean isTrackable(MCParticle p,Track t)
{
@@ -477,19 +405,6 @@
}
return tr;
}
- private MCParticle getTrackableParent(MCParticle p)
- {
- MCParticle ret = null;
- boolean htp = false;
- MCParticle par = p;
- while(!htp&&(par.getParents().size() == 1))
- {
- par = par.getParents().get(0);
- if(tmap.containsKey(par))htp = true;
- }
- if(htp)ret = par;
- return ret;
- }
private boolean hasTrackableDaughters(MCParticle p)
{
boolean htd = false;
@@ -517,70 +432,33 @@
}
return htd;
}
- private boolean tooClose(Hep3Vector v)
- {
- boolean tc = false;
- double vr = Math.sqrt(v.x()*v.x() + v.y()*v.y());
- if( (vr > emR-tccut)||(Math.abs(v.z()) > emZ-tccut) )
- {
- tc = true;
- }
- return tc;
- }
- private MCParticle getLongestNeutralParent(MCParticle p)
+ protected void detectorChanged(Detector detector)
{
- MCParticle ret = null;
- double length = minneutlen;
- MCParticle par = p;
- while(par.getParents().size() == 1)
+ super.detectorChanged(detector);
+
+ Map<String,Subdetector> subdets = detector.getSubdetectors();
+ for (Subdetector subdet : subdets.values())
{
- par = par.getParents().get(0);
- if(par.getCharge() != 0)
+ if (subdet.isCalorimeter())
{
- Hep3Vector v = par.getOrigin();
- Hep3Vector e = par.getEndPoint();
- if(tooClose(e))
- {
- double l = Math.sqrt( (e.x()-v.x())*(e.x()-v.x()) +
- (e.y()-v.y())*(e.y()-v.y()) +
- (e.z()-v.z())*(e.z()-v.z()) );
- if(l > length)
- {
- length = l;
- ret = par;
- }
+ if (subdet.getName().equals("EMBarrel"))
+ {
+ Tube tube = (Tube)subdet.getDetectorElement().getGeometry().getLogicalVolume().getSolid();
+ emR = tube.getInnerRadius();
+ double outerRadius = tube.getOuterRadius();
+ double zHalfLength = tube.getZHalfLength();
+ }
+ if (subdet.getName().equals("EMEndcap"))
+ {
+ Tube tube = (Tube)subdet.getDetectorElement().getChildren().get(0).getGeometry().getLogicalVolume().getSolid();
+ double zHalfLength = tube.getZHalfLength();
+ Hep3Vector v = subdet.getDetectorElement().getChildren().get(0).getGeometry().getPosition();
+ emZ = Math.abs(v.z()) - zHalfLength;
}
}
}
- return ret;
- }
- protected void detectorChanged(Detector detector)
- {
- super.detectorChanged(detector);
-
- Map<String,Subdetector> subdets = detector.getSubdetectors();
- for (Subdetector subdet : subdets.values())
- {
- if (subdet.isCalorimeter())
- {
- if (subdet.getName().equals("EMBarrel"))
- {
- Tube tube = (Tube)subdet.getDetectorElement().getGeometry().getLogicalVolume().getSolid();
- emR = tube.getInnerRadius();
- double outerRadius = tube.getOuterRadius();
- double zHalfLength = tube.getZHalfLength();
- }
- if (subdet.getName().equals("EMEndcap"))
- {
- Tube tube = (Tube)subdet.getDetectorElement().getChildren().get(0).getGeometry().getLogicalVolume().getSolid();
- double zHalfLength = tube.getZHalfLength();
- Hep3Vector v = subdet.getDetectorElement().getChildren().get(0).getGeometry().getPosition();
- emZ = Math.abs(v.z()) - zHalfLength;
- }
- }
- }
System.out.println("Detector change: innerZ = "+emZ+" innerR = "+emR);
-
- }
+
+ }
}
