lcsim/src/org/lcsim/contrib/SteveMagill
diff -u -r1.2 -r1.3
--- TrShowCompDriver.java 29 May 2008 22:37:14 -0000 1.2
+++ TrShowCompDriver.java 20 Aug 2008 20:25:42 -0000 1.3
@@ -15,9 +15,9 @@
import org.lcsim.event.Track;
import org.lcsim.util.hitmap.HitMap;
import org.lcsim.recon.cluster.util.BasicCluster;
-
+import org.lcsim.spacegeom.*;
import org.lcsim.event.ReconstructedParticle;
-import org.lcsim.event.base.MCReconstructedParticle;
+import org.lcsim.event.base.*;
import org.lcsim.event.ParticleID;
import org.lcsim.recon.cluster.cheat.CheatCluster;
import org.lcsim.event.ParticleID;
@@ -28,7 +28,8 @@
import hep.physics.particle.properties.ParticleType;
/**
- * Defines perfect particle flow objects - perfect charged particles, perfect photons and perfect neutral hadrons
+ * Diagnostic routine to calculate purities and efficiencies per particle and event
+ * for charged particles
*/
public class TrShowCompDriver extends Driver
@@ -42,15 +43,44 @@
public void process(EventHeader event)
{
- List<CheatCluster> chclusters = event.get(CheatCluster.class,"PerfectCheatClusters");
-
- // check cores
- int nhitsc = 0; // number of hits in core clusters
- int nchhitsc = 0; // number of charged particle hits in core clusters
+// List<CheatCluster> chclusters = event.get(CheatCluster.class,"PerfectCheatClusters");
+ List<ReconstructedParticle> perfrps = event.get(ReconstructedParticle.class,"PerfectRecoParticles");
+ List<Track> evtracks = event.get(Track.class, "PerfectTracks");
+ Map<Track, Integer> trilmap = (Map<Track, Integer>) event.get("TrackILMap");
+ Map<Track, BasicCluster> trclmap = (Map<Track, BasicCluster>) event.get("TrackClusMap");
+
+ // Check the cheated clusters from Perfect Reco list for eficiency den
+ int nhitsccc = 0;
+ double Eccc = 0.;
+ for (ReconstructedParticle perfrp : perfrps)
+ {
+ double rpcharge = perfrp.getCharge();
+ if (rpcharge!=0)
+ {
+ Hep3Vector pp = perfrp.getMomentum();
+ double pmag = pp.magnitude();
+ List<Cluster> rpclusters = perfrp.getClusters();
+ if (rpclusters.size()>0)
+ {
+ BasicCluster rpclus = (BasicCluster) rpclusters.get(0);
+ nhitsccc += rpclus.getSize();
+// Eccc += perfrp.getEnergy();
+ Eccc += rpclus.getEnergy();
+ aida.cloud1D("Perfect ChCl E over p").fill(rpclus.getEnergy()/pmag);
+ aida.cloud2D("Perfect ChCl E over p vs p").fill(pmag,rpclus.getEnergy()/pmag);
+ }
+ }
+ }
+ // get purities for cores and showers
+ // check cores first
+ int nhitsc = 0; // number of hits found in core clusters
+ int nchhitsc = 0; // number of charged particle hits for found hits in core clusters
+ int nphohitsc = 0; // number of photon hits for found hits in core clusters
+ int nnhhitsc = 0; // number of neutral hadron hits for found hits in core clusters
try
{
List<BasicCluster> tcclusters = event.get(BasicCluster.class,"TrackCoreClusters");
- aida.cloud1D("Number of TrCore Clusters").fill(tcclusters.size());
+// aida.cloud1D("Number of TrCore Clusters").fill(tcclusters.size());
for (BasicCluster tcclus : tcclusters)
{
int ncorehits = tcclus.getSize();
@@ -61,125 +91,99 @@
{
SimCalorimeterHit hit = (SimCalorimeterHit) tchit;
double hitch = Math.abs(hit.getMCParticle(0).getCharge());
+ double hitm = hit.getMCParticle(0).getMass();
if (hitch>0) nchhitsc++;
if (hitch>0) nchcorehits++;
+ if (hitch==0 && hitm==0) nphohitsc++;
+ if (hitch==0 && hitm>0) nnhhitsc++;
}
if (ncorehits>0)
{
double dncorehits = ncorehits;
double dnchcorehits = nchcorehits;
- double hrat = dnchcorehits/dncorehits;
- aida.cloud2D("NClus hits vs NMC hits TrCoreClus").fill(nchcorehits,ncorehits);
- aida.cloud1D("MC Particle Purity per TrCoreClus").fill(hrat);
+ double hrat = dnchcorehits/dncorehits; // purity per core
+ aida.cloud2D("NClus hits vs NMC hits per TrCoreClus").fill(nchcorehits,ncorehits);
+ aida.cloud1D("MC Particle Purity per TrCore").fill(hrat);
}
}
- double dnhitsc = nhitsc;
- double dnchhitsc = nchhitsc;
- double coreTot = dnchhitsc/dnhitsc;
- if (nhitsc>0) aida.cloud1D("Total Particle Purity for TrCores in event").fill(coreTot);
+ if (nhitsc>0)
+ {
+ double dnhitsc = nhitsc;
+ double dnchhitsc = nchhitsc;
+ double dnphohitsc = nphohitsc;
+ double dnnhhitsc = nnhhitsc;
+ aida.cloud1D("MC Particle Purity TrCores per event").fill(dnchhitsc/dnhitsc);
+ aida.cloud1D("Pho conf TrCores per event").fill(dnphohitsc/dnhitsc);
+ aida.cloud1D("NeuH conf TrCores per event").fill(dnnhhitsc/dnhitsc);
+ }
}
catch(java.lang.IllegalArgumentException ex)
{
- System.out.println("requested object not found in event " + _tsclusname);
- }
+ System.out.println("TrackCore Clusters not found in event");
+ }
- int NMipCl = 0;
+ // now same for showers (includes cores)
int NTrCALClus = 0;
- int nhitstm = 0;
- int nchhitstm = 0;
+ int nhitsts = 0;
+ int nchhitsts = 0;
+ int nphohitsts = 0;
+ int nnhhitsts = 0;
double trshclusE = 0.;
double trshclusP = 0.;
- double hratTot = 0.;
-
try
{
List<BasicCluster> tsclusters = event.get(BasicCluster.class,_tsclusname);
- aida.cloud1D("Number of TrCAL Clusters").fill(tsclusters.size());
+// aida.cloud1D("Number of TrCAL Clusters").fill(tsclusters.size());
for (BasicCluster tsclus : tsclusters)
{
- NMipCl++;
NTrCALClus++;
int numhits = tsclus.getSize();
- nhitstm += tsclus.getSize();
+ nhitsts += tsclus.getSize();
double tmclE = tsclus.getEnergy();
trshclusE += tsclus.getEnergy();
- aida.cloud1D("Tr CAL Clus E").fill(tmclE);
- aida.cloud1D("Tr CAL Clus Size").fill(numhits);
- aida.cloud2D("NumHits vs E TrCALClus").fill(tmclE,numhits);
- List<CalorimeterHit> tmcalhits = tsclus.getCalorimeterHits();
+ List<CalorimeterHit> tscalhits = tsclus.getCalorimeterHits();
int numchits = 0;
- for (CalorimeterHit tmhit : tmcalhits)
+ int jhits = 0;
+ for (CalorimeterHit tshit : tscalhits)
{
- SimCalorimeterHit hit = (SimCalorimeterHit) tmhit;
+ SimCalorimeterHit hit = (SimCalorimeterHit) tshit;
double hitch = Math.abs(hit.getMCParticle(0).getCharge());
+ double hitm = hit.getMCParticle(0).getMass();
+ jhits++;
if (hitch>0) numchits++;
- if (hitch>0) nchhitstm++;
+ if (hitch>0) nchhitsts++;
+ if (hitch==0 && hitm==0) nphohitsts++;
+ if (hitch==0 && hitm>0) nnhhitsts++;
}
if (numhits>0)
{
double dnumhits = numhits;
double dnumchits = numchits;
- double hrat = dnumchits/dnumhits;
- aida.cloud2D("NClus hits vs NMC hits TrCALClus").fill(numchits,numhits);
- aida.cloud1D("MC Particle Purity per TrCAL Clus").fill(hrat);
+ double djhits = jhits;
+ aida.cloud2D("NClus hits vs NMC hits per TrCALClus").fill(numchits,numhits);
+ aida.cloud1D("MC Particle Purity per TrCALClus").fill(dnumchits/dnumhits);
+ aida.cloud1D("Test hitsogram").fill(djhits/dnumhits);
}
}
- double dnhitstm = nhitstm;
- double dnchhitstm = nchhitstm;
- hratTot = dnchhitstm/dnhitstm;
- aida.cloud1D("Total Particle Purity for TrCAL in event").fill(hratTot);
- }
- catch(java.lang.IllegalArgumentException ex)
- {
- System.out.println("requested object not found in event " + _tsclusname);
- }
- if (NMipCl == 1)
+ if (nhitsts>0 && nhitsccc>0) // if both denominators are OK
{
- aida.cloud1D("Num Hits TrShow Clus single cluster").fill(nhitstm);
+ double dnhitsts = nhitsts;
+ double dnchhitsts = nchhitsts;
+ double dnphohitsts = nphohitsts;
+ double dnnhhitsts = nnhhitsts;
+ double dnhitsccc = nhitsccc;
+ aida.cloud1D("Purity of TrCALs per event").fill(dnchhitsts/dnhitsts);
+ aida.cloud1D("Pho conf TrCALs per event").fill(dnphohitsts/dnhitsts);
+ aida.cloud1D("NeuH conf TrCALs per event").fill(dnnhhitsts/dnhitsts);
+ aida.cloud1D("Efficiency of TrCALs per event").fill(dnchhitsts/dnhitsccc);
+ aida.cloud2D("Pur vs Eff TrCALs").fill(dnchhitsts/dnhitsccc,dnchhitsts/dnhitsts);
+ aida.cloud1D("Num of chhits missing from TrCALs").fill(dnhitsccc-dnchhitsts);
}
-
- // Check the cheated clusters also
- double chclusE = 0;
- int NCheatCl = 0;
- double PerfChHClE = 0;
- double MCChHE = 0;
- int nhitsingcl = 0;
- double chclchpE = 0;
- double chclchpP = 0;
- aida.cloud1D("Number of Cheated TrCAL Clusters").fill(chclusters.size());
- for (CheatCluster chclus : chclusters)
- {
- chclusE += chclus.getEnergy();
- aida.cloud1D("Num hits in Perfect Cheat Cluster TrCAL").fill(chclus.getSize());
- aida.cloud2D("Num hits vs E PerfCheatCl TrCAL").fill(chclus.getEnergy(),chclus.getSize());
- aida.cloud1D("Num hits over E PerfCheatCl TrCAL").fill(chclus.getSize()/chclus.getEnergy());
- if (chclus.getMCParticle().getCharge() != 0)
- {
- NCheatCl++;
- nhitsingcl += chclus.getSize();
- PerfChHClE += chclus.getEnergy();
- MCChHE += chclus.getMCParticle().getEnergy();
- chclchpE += chclus.getEnergy();
- chclchpP += chclus.getMCParticle().getMomentum().magnitude();
- aida.cloud1D("Cheated Cluster Charge E over P TrCAL").fill(chclchpE/chclchpP);
- }
- }
- if (NCheatCl==1)
- {
- aida.cloud1D("Single CheatClus E over p").fill(chclchpE/chclchpP);
- aida.cloud1D("Single CheatClus E Eff TrCALclus").fill(trshclusE/chclchpE);
- }
- if (NCheatCl > 0)
+ }
+ catch(java.lang.IllegalArgumentException ex)
{
- aida.cloud1D("Cheat Cluster Charge E over P TrCAL").fill(chclchpE/chclchpP);
- aida.cloud1D("Cheated Cluster E single cluster TrCAL").fill(chclusE);
- aida.cloud1D("Num Hits in Cheated Cluster single cluster TrCAL").fill(nhitsingcl);
- double dnhitstm = nhitstm;
- double dnhitsingcl = nhitsingcl;
- aida.cloud1D("NHits Eff TrCALClus CheatClus").fill(dnhitstm/dnhitsingcl);
- aida.cloud1D("E Eff TrCALClus CheatClus").fill(1.176*trshclusE/chclchpE);
- aida.cloud2D("Charged Cluster EnEff vs hitP").fill(hratTot,1.176*trshclusE/chclchpE);
+ System.out.println("requested object not found in event " + _tsclusname);
}
}
