slic/analysis
diff -N Slic_DiagnosticHists.java
--- Slic_DiagnosticHists.java 18 Mar 2005 01:49:06 -0000 1.6
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,699 +0,0 @@
-import hep.aida.*;
-import hep.lcio.event.*;
-import org.freehep.record.loop.event.*;
-import java.util.Vector;
-
-public class Slic_DiagnosticHists extends RecordAdapter
-{
- public Slic_DiagnosticHists()
- {
- IAnalysisFactory analysisFactory = IAnalysisFactory.create();
- ITree tree = analysisFactory.createTreeFactory().create();
-
- // Create histogram here
- IHistogramFactory histogramFactory = analysisFactory.createHistogramFactory(tree);
-
- tree.mkdir("/Temp");
- tree.cd("/Temp");
-
- hVEdistA = histogramFactory.createCloud1D("Vtx to EP dist: KE > 1 MeV");
- hVEdistB = histogramFactory.createCloud1D("Vtx to EP dist: KE < 1 MeV");
- hVE0PDG2 = histogramFactory.createCloud1D("Vtx to EP dist = 0: Predecay PDG");
- hVE0PDG5 = histogramFactory.createCloud1D("Vtx to EP dist = 0: Sim decay PDG");
- hVE0ci = histogramFactory.createHistogram1D("Vtx to EP dist = 0: creation index",6,-0.5,5.5);
- hVE0di = histogramFactory.createHistogram1D("Vtx to EP dist = 0: demise index",6,-0.5,5.5);
-
- tree.mkdir("/Event");
- tree.cd("/Event");
-
- hnMCP = histogramFactory.createCloud1D("# MCParticles");
- hFSE = histogramFactory.createCloud1D("Gen Final State Energy");
- hSimFSE = histogramFactory.createCloud1D("Sim Final State Energy");
-
- tree.mkdir("/Particles");
- tree.cd("/Particles");
-
- hsimstat = histogramFactory.createHistogram1D("Simulator Status",128,-.5,127.5);
- hgenstat = histogramFactory.createHistogram1D("Generator Status",4,-.5,3.5);
- pgenvssim = histogramFactory.createHistogram2D("Status:gen vs sim",4,-.5,3.5,128,-.5,127.5);
- hstatbits = histogramFactory.createHistogram1D("Simulator Status bits",7,.5,7.5);
- hgFSE = histogramFactory.createCloud1D("Energy GenFS particles");
- hgFSPDG = histogramFactory.createCloud1D("PDG GenFS particles");
- hsFSE = histogramFactory.createCloud1D("Energy SimFS particles");
- hsFSPDG = histogramFactory.createCloud1D("PDG SimFS particles");
- int nccat = 6;
- int necat = 5;
- hEce = new ICloud1D[nccat][necat];
- hPDGce = new ICloud1D[nccat][necat];
- hVX = new ICloud2D[nccat][necat];
- hEP = new ICloud2D[nccat][necat];
- hEc = new ICloud1D[nccat];
- hPDGc = new ICloud1D[nccat];
- String[] cname = {"Document","Initial","PreDecay","Backscatter","Continuous","Decay"};
- String[] ename = {"None","DecayT","DecayC","Left","Stopped"};
- for(int i=0;i<nccat;i++)
- {
- hEc[i] = histogramFactory.createCloud1D("Type "+cname[i] + ":: Energy");
- hPDGc[i] = histogramFactory.createCloud1D("Type "+cname[i] + ":: PDG");
- for(int j=0;j<necat;j++)
- {
- hEce[i][j] = histogramFactory.createCloud1D("Type "+cname[i] + ": End "+ename[j]+":: Energy");
- hPDGce[i][j] = histogramFactory.createCloud1D("Type "+cname[i] + ": End "+ename[j]+ ":: PDG");
- hVX[i][j] = histogramFactory.createCloud2D("Type "+cname[i] + ": End "+ename[j]+":: Vtx");
- hEP[i][j] = histogramFactory.createCloud2D("Type "+cname[i] + ": End "+ename[j]+":: Endpoint");
- }
- }
- hbadstatGB = histogramFactory.createCloud1D("bad status: Gen particle with Backscatter");
- hbadstatGVNE = histogramFactory.createCloud1D("bad status: Gen particle with VtxnotEPofParent");
- hbadstatGFnoE = histogramFactory.createCloud1D("bad status: Gen particle in Sim with no End status");
- hbadstatSnotS = histogramFactory.createCloud1D("bad status: Sim particle without Sim set");
- hbadstatSnoE = histogramFactory.createCloud1D("bad status: Sim particle with no end status");
- hbadstatME12 = histogramFactory.createCloud1D("mult end status: Dec in T + Decay in C");
- hbadstatME13 = histogramFactory.createCloud1D("mult end status: Dec in T + Left");
- hbadstatME14 = histogramFactory.createCloud1D("mult end status: Dec in T + Stopped");
- hbadstatME23 = histogramFactory.createCloud1D("mult end status: Dec in C + Left");
- hbadstatME24 = histogramFactory.createCloud1D("mult end status: Dec in C + Stopped");
- hbadstatME34 = histogramFactory.createCloud1D("mult end status: Left + Stopped");
-
- hnCal = new ICloud1D[ncalcolls][3];
- hECal = new ICloud1D[ncalcolls][3];
- RZCal = new ICloud2D[ncalcolls][3];
- hlayCal = new IHistogram1D[ncalcolls][3];
- hphibinCal = new ICloud1D[ncalcolls][3];
- hthetabinCal = new ICloud1D[ncalcolls][3];
- hncontperhit = new IHistogram1D[ncalcolls][3];
- hphivsbinCal = new ICloud2D[ncalcolls][3];
- hthetavsbinCal = new ICloud2D[ncalcolls][3];
- hcontPDG = new ICloud1D[ncalcolls][3];
- hcontPDGsingle = new ICloud1D[ncalcolls][3];
- nCalhits = new int[ncalcolls][3];
- for(int i=0;i<ncalcolls;i++)
- {
- tree.mkdir("/"+calcoll[i]);
- for(int j=0;j<3;j++)
- {
- tree.mkdir("/"+calcoll[i]+"/"+which[j]);
- tree.cd("/"+calcoll[i]+"/"+which[j]);
- nCalhits[i][j] = 0;
- hnCal[i][j] = histogramFactory.createCloud1D(" #hits");
- hECal[i][j] = histogramFactory.createCloud1D(" Energy per hit");
- RZCal[i][j] = histogramFactory.createCloud2D(" R vs Z");
- hlayCal[i][j] = histogramFactory.createHistogram1D(" layer",50,-.5,49.5);
- hphibinCal[i][j] = histogramFactory.createCloud1D(" phi bin");
- hthetabinCal[i][j] = histogramFactory.createCloud1D(" theta bin");
- hphivsbinCal[i][j] = histogramFactory.createCloud2D(" phibin vs phi");
- hthetavsbinCal[i][j] = histogramFactory.createCloud2D(" thetabin vs theta");
- hncontperhit[i][j] = histogramFactory.createHistogram1D("#parts per hit",20,0.5,20.5);
- hcontPDG[i][j] = histogramFactory.createCloud1D("PDG of contributors");
- hcontPDGsingle[i][j] = histogramFactory.createCloud1D("single contibutor PDG");
- }
- }
- hnTrk = new ICloud1D[ntrkcolls][3];
- hETrk = new ICloud1D[ntrkcolls][3];
- RZTrk = new ICloud2D[ntrkcolls][3];
- hlayTrk = new IHistogram1D[ntrkcolls][3];
- nTrkhits = new int[ntrkcolls][3];
- for(int i=0;i<ntrkcolls;i++)
- {
- tree.mkdir("/"+trkcoll[i]);
- for(int j=0;j<3;j++)
- {
- tree.mkdir("/"+trkcoll[i]+"/"+which[j]);
- tree.cd("/"+trkcoll[i]+"/"+which[j]);
- nTrkhits[i][j] = 0;
- hnTrk[i][j] = histogramFactory.createCloud1D(" #hits");
- hETrk[i][j] = histogramFactory.createCloud1D(" Energy per hit");
- RZTrk[i][j] = histogramFactory.createCloud2D(" R vs Z");
- hlayTrk[i][j] = histogramFactory.createHistogram1D(" layer",50,-.5,49.5);
- }
- }
-
-
-
- nmaxhits = 99000;
- ievt = 0;
-
- }
- public void recordSupplied(RecordSuppliedEvent recordSuppliedEvent)
- {
- LCEvent event = (LCEvent) recordSuppliedEvent.getRecord();
-
- LCCollection collection = event.getCollection("MCParticle");
- int nMc = collection.getNumberOfElements();
-// System.out.println(" # MCParticles = " + nMc);
- MCParticle testmc = (MCParticle) collection.getElementAt(nMc-1);
-// if(testmc.getEnergy() < 9.)return;
- boolean done = true;
-/*
- if( ievt == 4708 )
- {
- for(int i = 0; i < nMc; i++)
- {
- dumpMCParticle( (MCParticle) collection.getElementAt(i),i);
- }
- }
- if(done)
- {
- ievt++;
- return;
- }
-*/
- MCParticle[] gFS = new MCParticle[nMc];
- int ngFS =0;
- hnMCP.fill(nMc);
- double totFSE = 0.;
- double totSimFSE = 0.;
- for(int i = 0; i < nMc; i++)
- {
- MCParticle mcparticle = (MCParticle) collection.getElementAt(i);
-// dumpMCParticle(mcparticle,i);
- if(mcparticle.isCreatedInSimulation())hstatbits.fill(1);
- if(mcparticle.isBackscatter())hstatbits.fill(2);
- if(mcparticle.vertexIsNotEndpointOfParent())hstatbits.fill(3);
- if(mcparticle.isDecayedInTracker())hstatbits.fill(4);
- if(mcparticle.isDecayedInCalorimeter())hstatbits.fill(5);
- if(mcparticle.hasLeftDetector())hstatbits.fill(6);
- if(mcparticle.isStopped())hstatbits.fill(7);
- double[] ep = mcparticle.getEndpoint();
- double R = Math.sqrt(ep[0]*ep[0]+ep[1]*ep[1]);
- double Z = ep[2];
- double[] vx = mcparticle.getVertex();
- double vR = Math.sqrt(vx[0]*vx[0]+vx[1]*vx[1]);
- double vZ = vx[2];
- double E = mcparticle.getEnergy();
- int PDG = mcparticle.getPDG();
- int gstat = mcparticle.getGeneratorStatus();
- int sstat = mcparticle.getSimulatorStatus();
- sstat &= 0x7fffffff;
- sstat = sstat >> 24;
- hsimstat.fill(sstat);
- hgenstat.fill(gstat);
- pgenvssim.fill(gstat,sstat);
-//
-// Check for inconsistencies
-//
- boolean badstat = checkStatusConsistency(mcparticle);
- if(badstat)
- {
- dumpMCParticle(mcparticle,i);
- continue;
- }
-
-//
-// Catagorize the particles
-//
- boolean fromGenerator = gstat > 0;
- boolean createdInSim = mcparticle.isCreatedInSimulation();
- int nendd = 0;
- int nd = mcparticle.getNumberOfDaughters();
- for(int k=0;k<nd;k++)
- {
- MCParticle d = mcparticle.getDaughter(k);
- if(d.vertexIsNotEndpointOfParent())continue;
- if(d.isBackscatter())continue;
- nendd++;
- }
- boolean hasEPDaughters = nendd > 0;
- int creationIndex = 0;
- int demiseIndex = 0;
- if(mcparticle.isDecayedInTracker())demiseIndex = 1;
- else if(mcparticle.isDecayedInCalorimeter())demiseIndex = 2;
- else if(mcparticle.hasLeftDetector())demiseIndex = 3;
- else if(mcparticle.isStopped())demiseIndex = 4;
- if(fromGenerator)
- {
- int np = mcparticle.getNumberOfParents();
- boolean parentinsim = false;
- if(np == 1)
- {
- MCParticle p = mcparticle.getParent(0);
- if(p.isDecayedInTracker() || p.isDecayedInCalorimeter() ||
- p.hasLeftDetector() || p.isStopped() )
- {
- parentinsim = true;
- }
- }
- if( demiseIndex == 0)
- {
-//
-// Simulator never saw this
-//
- creationIndex = 0;
- }
- else if(!createdInSim && !parentinsim)
- {
-//
-// Initial generator particle passed to simulator
-//
- creationIndex = 1;
- }
- else
- {
-//
-// Predecay
-//
- creationIndex = 2;
- }
- }
- else
- {
-//
-// Simulator created it on its own
-//
- if(mcparticle.isBackscatter())
- {
-//
-// backscatter
-//
- creationIndex = 3;
- }
- else if(mcparticle.vertexIsNotEndpointOfParent())
- {
-//
-// created without killing parent
-//
- creationIndex = 4;
- }
- else
- {
-//
-// result of parent decay or interaction
-//
- creationIndex = 5;
- }
- }
- hEce[creationIndex][demiseIndex].fill(E);
- hPDGce[creationIndex][demiseIndex].fill(PDG);
- hVX[creationIndex][demiseIndex].fill(vZ,vR);
- hEP[creationIndex][demiseIndex].fill(Z,R);
- hEc[creationIndex].fill(E);
- hPDGc[creationIndex].fill(PDG);
- if(creationIndex > 0)
- {
- double KE = E - mcparticle.getMass();
- double VEdist = Math.sqrt( (vx[0] - ep[0])*(vx[0] - ep[0]) +
- (vx[1] - ep[1])*(vx[1] - ep[1]) +
- (vx[2] - ep[2])*(vx[2] - ep[2]) );
- if(KE < .001)hVEdistB.fill(VEdist);
- else hVEdistA.fill(VEdist);
-// if( (KE >= .001) && (VEdist < .001) )dumpMCParticle(mcparticle,i);
-// dumpMCParticle(mcparticle,i);
- if(VEdist == 0.)
- {
- if(creationIndex == 2)hVE0PDG2.fill(PDG);
- else if(creationIndex == 5)hVE0PDG5.fill(PDG);
- hVE0ci.fill(creationIndex);
- hVE0di.fill(demiseIndex);
- if(PDG == 3122)
- {
- for(int j = 0; j < nMc; j++)
- {
- dumpMCParticle( (MCParticle) collection.getElementAt(j),j);
- }
- }
- }
- }
- if(creationIndex == 1)
- {
-//
-// Generator "Final State"
-//
- hgFSE.fill(E);
- hgFSPDG.fill(PDG);
- totFSE += E;
- gFS[ngFS] = mcparticle;
- ngFS++;
- }
- }
- hFSE.fill(totFSE);
- for(int i=0;i<ngFS;i++)
- {
- MCParticle[] sFS = getSimFSparticles(gFS[i]);
- for(int j=0;j<sFS.length;j++)
- {
- hsFSE.fill(sFS[j].getEnergy());
- hsFSPDG.fill(sFS[j].getPDG());
- totSimFSE += sFS[j].getEnergy();
- }
- }
- hSimFSE.fill(totSimFSE);
-
- for(int i=0;i<ncalcolls;i++)
- {
- collection = event.getCollection(calcoll[i]);
- int nhits = collection.getNumberOfElements();
- hnCal[i][0].fill(nhits);
- int n1 = 0;
- int n2 = 0;
- for(int j = 0; j < nhits; j++)
- {
- SimCalorimeterHit hit = (SimCalorimeterHit) collection.getElementAt(j);
- float[] pos = hit.getPosition();
- double R = Math.sqrt(pos[0]*pos[0] + pos[1]*pos[1]);
- int id0 = hit.getCellID0();
- int id1 = hit.getCellID1();
- int lay = id0 & 0x7F;
- int tbin = id1 & 0x7FF;
- int pbin = (id1 >> 11) & 0x7FF;
- double phi = Math.atan2(pos[1],pos[0]);
- if(phi < 0)phi += 2.*Math.PI;
- double theta = Math.atan2(pos[2],R);
- if(nCalhits[i][0] < nmaxhits)
- {
- nCalhits[i][0]++;
- RZCal[i][0].fill(pos[2],R);
- hECal[i][0].fill(hit.getEnergy());
- hthetabinCal[i][0].fill(tbin);
- hphibinCal[i][0].fill(pbin);
- hphivsbinCal[i][0].fill(phi,pbin);
- hthetavsbinCal[i][0].fill(theta,tbin);
- for(int m=0;m<hit.getNMCParticles();m++)
- {
- hcontPDG[i][0].fill(hit.getPDGCont(m));
- if(hit.getNMCParticles() == 1)
- {
- hcontPDGsingle[i][0].fill(hit.getPDGCont(m));
- }
- }
- }
- hlayCal[i][0].fill(lay);
- hncontperhit[i][0].fill(hit.getNMCParticles());
- int k = 1;
- if(hit.getEnergy() > calEcut[i])
- {
- k = 2;
- n2++;
- }
- else n1++;
- if(nCalhits[i][k] < nmaxhits)
- {
- nCalhits[i][k]++;
- RZCal[i][k].fill(pos[2],R);
- hECal[i][k].fill(hit.getEnergy());
- hthetabinCal[i][k].fill(tbin);
- hphibinCal[i][k].fill(pbin);
- hphivsbinCal[i][k].fill(phi,pbin);
- hthetavsbinCal[i][k].fill(theta,tbin);
- for(int m=0;m<hit.getNMCParticles();m++)
- {
- hcontPDG[i][k].fill(hit.getPDGCont(m));
- if(hit.getNMCParticles() == 1)
- {
- hcontPDGsingle[i][k].fill(hit.getPDGCont(m));
- }
- }
- }
- hlayCal[i][k].fill(lay);
- hncontperhit[i][k].fill(hit.getNMCParticles());
- }
- hnCal[i][1].fill(n1);
- hnCal[i][2].fill(n2);
- }
- for(int i=0;i<ntrkcolls;i++)
- {
- collection = event.getCollection(trkcoll[i]);
- int nhits = collection.getNumberOfElements();
- hnTrk[i][0].fill(nhits);
- int n1 = 0;
- int n2 = 0;
- for(int j = 0; j < nhits; j++)
- {
- SimTrackerHit hit = (SimTrackerHit) collection.getElementAt(j);
- double[] pos = hit.getPosition();
- double R = Math.sqrt(pos[0]*pos[0] + pos[1]*pos[1]);
- if(nTrkhits[i][0] < nmaxhits)
- {
- nTrkhits[i][0]++;
- RZTrk[i][0].fill(pos[2],R);
- hETrk[i][0].fill(hit.getdEdx());
- }
- hlayTrk[i][0].fill(hit.getCellID() & 0x7F);
- int k = 1;
- if(hit.getdEdx() > trkEcut[i])
- {
- k = 2;
- n2++;
- }
- else n1++;
- if(nTrkhits[i][k] < nmaxhits)
- {
- nTrkhits[i][k]++;
- RZTrk[i][k].fill(pos[2],R);
- hETrk[i][k].fill(hit.getdEdx());
- }
- hlayTrk[i][k].fill(hit.getCellID() & 0x7F);
- }
- hnTrk[i][1].fill(n1);
- hnTrk[i][2].fill(n2);
- }
- ievt++;
- }
-
- private MCParticle[] getSimFSparticles(MCParticle gpart)
- {
- Vector fsp = new Vector();
- traceit(gpart, fsp);
- MCParticle[] returnp = new MCParticle[fsp.size()];
- for(int i=0;i<fsp.size();i++)
- {
- returnp[i] = (MCParticle) fsp.get(i);
- }
- return returnp;
- }
- private void traceit(MCParticle p, Vector v)
- {
- int nepd = 0;
- int nd = p.getNumberOfDaughters();
- for(int k=0;k<nd;k++)
- {
- MCParticle d = p.getDaughter(k);
- if(d.vertexIsNotEndpointOfParent())continue;
- if(d.isBackscatter())continue;
- nepd++;
- }
- if(nepd > 0)
- {
- for(int i=0;i<nd;i++)
- {
- MCParticle dp = p.getDaughter(i);
- if(!dp.isBackscatter())traceit(dp,v);
- }
- }
- else
- {
- v.add(p);
- }
- }
-
- private boolean checkStatusConsistency(MCParticle mcparticle)
- {
- boolean badstat = false;
- int gstat = mcparticle.getGeneratorStatus();
- if(gstat > 0)
- {
- if(mcparticle.isBackscatter())
- {
- hbadstatGB.fill(mcparticle.getPDG());
- badstat = true;
- }
- if(mcparticle.vertexIsNotEndpointOfParent())
- {
- hbadstatGVNE.fill(mcparticle.getPDG());
- badstat = true;
- }
- if((gstat == 1) || (mcparticle.isCreatedInSimulation()))
- {
- if(!(mcparticle.isDecayedInTracker()||mcparticle.isDecayedInCalorimeter()||
- mcparticle.hasLeftDetector()||mcparticle.isStopped()))
- {
- hbadstatGFnoE.fill(mcparticle.getPDG());
-// badstat = true;
- }
- else
- {
- if(mcparticle.isDecayedInTracker())
- {
- if(mcparticle.isDecayedInCalorimeter())
- {
- hbadstatME12.fill(mcparticle.getPDG());
- badstat = true;
- }
- if(mcparticle.hasLeftDetector())
- {
- hbadstatME13.fill(mcparticle.getPDG());
- badstat = true;
- }
- if(mcparticle.isStopped())
- {
- hbadstatME14.fill(mcparticle.getPDG());
- badstat = true;
- }
- }
- else if(mcparticle.isDecayedInCalorimeter())
- {
- if(mcparticle.hasLeftDetector())
- {
- hbadstatME23.fill(mcparticle.getPDG());
- badstat = true;
- }
- if(mcparticle.isStopped())
- {
- hbadstatME24.fill(mcparticle.getPDG());
- badstat = true;
- }
- }
- else if(mcparticle.hasLeftDetector())
- {
- if(mcparticle.isStopped())
- {
- hbadstatME34.fill(mcparticle.getPDG());
- badstat = true;
- }
- }
- }
- }
- }
- else
- {
- if(!mcparticle.isCreatedInSimulation())
- {
- hbadstatSnotS.fill(mcparticle.getPDG());
- badstat = true;
- }
- else
- {
- if(!(mcparticle.isDecayedInTracker()||mcparticle.isDecayedInCalorimeter()||
- mcparticle.hasLeftDetector()||mcparticle.isStopped()))
- {
- hbadstatSnoE.fill(mcparticle.getPDG());
-// badstat = true;
- }
- else
- {
- if(mcparticle.isDecayedInTracker())
- {
- if(mcparticle.isDecayedInCalorimeter())
- {
- hbadstatME12.fill(mcparticle.getPDG());
- badstat = true;
- }
- if(mcparticle.hasLeftDetector())
- {
- hbadstatME13.fill(mcparticle.getPDG());
- badstat = true;
- }
- if(mcparticle.isStopped())
- {
- hbadstatME14.fill(mcparticle.getPDG());
- badstat = true;
- }
- }
- else if(mcparticle.isDecayedInCalorimeter())
- {
- if(mcparticle.hasLeftDetector())
- {
- hbadstatME23.fill(mcparticle.getPDG());
- badstat = true;
- }
- if(mcparticle.isStopped())
- {
- hbadstatME24.fill(mcparticle.getPDG());
- badstat = true;
- }
- }
- else if(mcparticle.hasLeftDetector())
- {
- if(mcparticle.isStopped())
- {
- hbadstatME34.fill(mcparticle.getPDG());
- badstat = true;
- }
- }
- }
- }
- }
-
- return badstat;
- }
- private void dumpMCParticle(MCParticle p, int i)
- {
- System.out.println("Event "+ievt+" particle "+i+" PDG = "+p.getPDG()+" E = "+p.getEnergy());
- System.out.println("Sim status: Created in Sim = "+p.isCreatedInSimulation()+" Backscatter = "+p.isBackscatter()+
- " Continuous process creation = "+p.vertexIsNotEndpointOfParent()+" Int/decay in tracker = "+
- p.isDecayedInTracker()+" Int/decay in cal = "+p.isDecayedInCalorimeter()+" Stopped = "+
- p.isStopped()+" Left det = "+p.hasLeftDetector());
- }
-
- private ICloud1D hVEdistA;
- private ICloud1D hVEdistB;
- private ICloud1D hVE0PDG2;
- private ICloud1D hVE0PDG5;
- private IHistogram1D hVE0ci;
- private IHistogram1D hVE0di;
-
- private IHistogram1D hsimstat;
- private IHistogram1D hgenstat;
- private IHistogram1D hstatbits;
- private ICloud1D hnMCP;
- private ICloud1D hFSE;
- private ICloud1D hSimFSE;
- private IHistogram2D pgenvssim;
- private ICloud1D hgFSE;
- private ICloud1D hgFSPDG;
- private ICloud1D hsFSE;
- private ICloud1D hsFSPDG;
- private ICloud1D[][] hEce;
- private ICloud1D[][] hPDGce;
- private ICloud2D[][] hVX;
- private ICloud2D[][] hEP;
- private ICloud1D[] hEc;
- private ICloud1D[] hPDGc;
- private ICloud1D hbadstatGB;
- private ICloud1D hbadstatGVNE;
- private ICloud1D hbadstatGFnoE;
- private ICloud1D hbadstatSnotS;
- private ICloud1D hbadstatSnoE;
- private ICloud1D hbadstatME12;
- private ICloud1D hbadstatME13;
- private ICloud1D hbadstatME14;
- private ICloud1D hbadstatME23;
- private ICloud1D hbadstatME24;
- private ICloud1D hbadstatME34;
-
- private ICloud1D[][] hnCal;
- private ICloud1D[][] hECal;
- private ICloud2D[][] RZCal;
- private IHistogram1D[][] hlayCal;
- private ICloud1D[][] hphibinCal;
- private ICloud1D[][] hthetabinCal;
- private ICloud2D[][] hphivsbinCal;
- private ICloud2D[][] hthetavsbinCal;
- private IHistogram1D[][] hncontperhit;
- private ICloud1D[][] hcontPDG;
- private ICloud1D[][] hcontPDGsingle;
- private ICloud1D[][] hnTrk;
- private ICloud1D[][] hETrk;
- private ICloud2D[][] RZTrk;
- private IHistogram1D[][] hlayTrk;
-
- private int[][] nCalhits;
- private int[][] nTrkhits;
-
- private int nmaxhits;
- private int nEMhits;
- private int nHADhits;
- private int nLUMhits;
- private int nMUONhits;
- private int nCENhits;
- private int nVXDhits;
- private String[] statflags = {"Stopped","Left","CalDecay","TrkDecay","Cont","Backscatter","SimCreated"};
- private String[] calcoll = {"EcalBarrHits","EcalEndcapHits","HcalBarrHits","HcalEndcapHits",
- "MuonBarrHits","MuonEndcapHits"};
- private String[] trkcoll = {"VtxBarrHits","TrkBarrHits","TrkEndcapHits"};
- private String[] which = {"All: ","E<cut: ","E>cut: "};
- private int ncalcolls = 6;
- private int ntrkcolls = 3;
- private double calEcut[] = {.00006,.00006,.001,.001,.000001,.000001};
- private double trkEcut[] = {.000015,.00003,.00003};
- private int ievt;
-}
slic/analysis
diff -N Slic_HitAssociation.java
--- Slic_HitAssociation.java 14 Mar 2005 18:42:50 -0000 1.2
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,398 +0,0 @@
-import hep.aida.*;
-import hep.lcio.event.*;
-import org.freehep.record.loop.event.*;
-import java.util.Vector;
-
-public class Slic_HitAssociation extends RecordAdapter
-{
- public Slic_HitAssociation()
- {
- IAnalysisFactory analysisFactory = IAnalysisFactory.create();
- ITree tree = analysisFactory.createTreeFactory().create();
-
- // Create histogram here
- IHistogramFactory histogramFactory = analysisFactory.createHistogramFactory(tree);
-
- tree.mkdir("/hists");
- tree.cd("/hists");
- hPDGall = histogramFactory.createHistogram1D("PDG all visible particles",5555,-.5,5554.5);
- hdIndex = histogramFactory.createHistogram1D("demise Index all visible particles",5,-0.5,4.5);
- hcIndex = histogramFactory.createHistogram1D("creation Index all visible particles",6,-0.5,5.5);
- hccIndex = histogramFactory.createHistogram1D("creation Index particles with cal hits",6,-0.5,5.5);
-
- hnthits = new ICloud1D[nparticledivisions][ncutdivisions];
- hnchits = new ICloud1D[nparticledivisions][ncutdivisions];
- hE = new ICloud1D[nparticledivisions][ncutdivisions];
- hEME = new ICloud1D[nparticledivisions][ncutdivisions];
- hnEMhits = new ICloud1D[nparticledivisions][ncutdivisions];
- hnHADhits = new ICloud1D[nparticledivisions][ncutdivisions];
- hEvsEME = new ICloud2D[nparticledivisions][ncutdivisions];
- hEvsHADE = new ICloud2D[nparticledivisions][ncutdivisions];
- hEvsCALE = new ICloud2D[nparticledivisions][ncutdivisions];
- entries = new int[nparticledivisions][ncutdivisions];
- for(int i=0;i<nparticledivisions;i++)
- {
- tree.mkdir("/"+pdnames[i]);
- for(int j=0;j<ncutdivisions;j++)
- {
- tree.mkdir("/"+pdnames[i]+"/"+cutnames[j]);
- tree.cd("/"+pdnames[i]+"/"+cutnames[j]);
- entries[i][j] = 0;
- hnthits[i][j] = histogramFactory.createCloud1D("# tracker hits");
- hnchits[i][j] = histogramFactory.createCloud1D("# calorimeter hits");
- hE[i][j] = histogramFactory.createCloud1D("creation Energy");
- hEME[i][j] = histogramFactory.createCloud1D("Energy in EM");
- hnEMhits[i][j] = histogramFactory.createCloud1D("#hits in EM");
- hEvsEME[i][j] = histogramFactory.createCloud2D("Energy vs EM energy");
- hnHADhits[i][j] = histogramFactory.createCloud1D("#hits in HAD");
- hEvsHADE[i][j] = histogramFactory.createCloud2D("Energy vs HAD energy");
- hEvsCALE[i][j] = histogramFactory.createCloud2D("Energy vs CAL energy");
- }
- }
- ng = 0;
- np = 0;
- ievt = 0;
- sfEM = .012;
- sfHAD = .06;
- }
- public void recordSupplied(RecordSuppliedEvent recordSuppliedEvent)
- {
- LCEvent event = (LCEvent) recordSuppliedEvent.getRecord();
- Vector MCPwithhits = new Vector();
- Vector ThitColl = new Vector();
- Vector ChitColl = new Vector();
- int nMCwithhits = 0;
-
- for(int ic = 0;ic<nTcollections;ic++)
- {
- LCCollection collection = event.getCollection(TCname[ic]);
- int n = collection.getNumberOfElements();
- for(int i = 0; i < n; i++)
- {
- SimTrackerHit hit = (SimTrackerHit) collection.getElementAt(i);
- MCParticle mcp = hit.getMCParticle();
- TCluster tc = null;
- for(int ip=0;ip<nMCwithhits;ip++)
- {
- MCParticle mcs = (MCParticle) MCPwithhits.get(ip);
- if(mcs == mcp)
- {
- tc = (TCluster) ThitColl.get(ip);
- ip = nMCwithhits;
- }
- }
- if(tc == null)
- {
- tc = new TCluster();
- CCluster cc = new CCluster();
- MCPwithhits.add(mcp);
- ThitColl.add(tc);
- ChitColl.add(cc);
- nMCwithhits++;
- }
- tc.addhit(hit);
- }
- }
-
- for(int ic = 0;ic<nCcollections;ic++)
- {
- LCCollection collection = event.getCollection(CCname[ic]);
- int n = collection.getNumberOfElements();
- for(int i = 0; i < n; i++)
- {
- SimCalorimeterHit hit = (SimCalorimeterHit) collection.getElementAt(i);
- int ncont = hit.getNMCParticles();
- double emax = 0.;
- MCParticle mcp = null;
- for(int j=0;j<ncont;j++)
- {
- if(hit.getEnergyCont(j) > emax)
- {
- mcp = hit.getParticleCont(j);
- emax = hit.getEnergyCont(j);
- }
- }
- CCluster cc = null;
- for(int ip=0;ip<nMCwithhits;ip++)
- {
- MCParticle mcs = (MCParticle) MCPwithhits.get(ip);
- if(mcs == mcp)
- {
- cc = (CCluster) ChitColl.get(ip);
- ip = nMCwithhits;
- }
- }
- if(cc == null)
- {
- cc = new CCluster();
- TCluster tc = new TCluster();
- MCPwithhits.add(mcp);
- ThitColl.add(tc);
- ChitColl.add(cc);
- nMCwithhits++;
- }
- cc.addhit(hit);
- }
- }
- for(int indx = 0;indx<nMCwithhits;indx++)
- {
- MCParticle mc = (MCParticle) MCPwithhits.get(indx);
- int pdg = Math.abs(mc.getPDG());
- hPDGall.fill(pdg);
- int demiseIndex = 0;
- if(mc.isDecayedInTracker())demiseIndex = 1;
- else if(mc.isDecayedInCalorimeter())demiseIndex = 2;
- else if(mc.hasLeftDetector())demiseIndex = 3;
- else if(mc.isStopped())demiseIndex = 4;
- int creationIndex = getCreationIndex(mc,demiseIndex);
- hcIndex.fill(creationIndex);
- hdIndex.fill(demiseIndex);
- CCluster cc = (CCluster) ChitColl.get(indx);
- TCluster tc = (TCluster) ThitColl.get(indx);
- int nthits = tc.getNhits();
- int nchits = cc.getNhits();
- if(nchits > 0)
- {
- hccIndex.fill(creationIndex);
- }
- double E = mc.getEnergy() - mc.getMass();
- if(creationIndex == 0)
- {
- System.out.println(" Event "+ievt);
- System.out.println(" Doc particle with hits!?!?!? PDG = "+pdg+" KE = "+E+" nThits = "+nthits+" nChits = "+
- nchits);
- }
- SimCalorimeterHit[] chits = cc.getHits();
- SimTrackerHit[] thits = tc.getHits();
- double EMenergy = 0;
- int nEMhits = 0;
- double HADenergy = 0;
- int nHADhits = 0;
- for(int i=0;i<nchits;i++)
- {
- int cellid = chits[i].getCellID0();
- int sys = (cellid >> 7) & 7;
- if(sys == 4)
- {
- EMenergy += chits[i].getEnergy();
- nEMhits ++;
- }
- if(sys == 5)
- {
- HADenergy += chits[i].getEnergy();
- nHADhits ++;
- }
- }
- double CALenergy = EMenergy/sfEM + HADenergy/sfHAD;
- int ind = -1;
- for(int k=0;k<nparticledivisions-2;k++)
- {
- if(pdg == pdpdg[k])
- {
- ind = k;
- continue;
- }
- }
- if(ind < 0)
- {
- ind = nparticledivisions - 1;
- if(mc.getCharge() != 0)ind --;
- }
- if(entries[ind][0] < maxentries)
- {
- entries[ind][0]++;
- hnthits[ind][0].fill(nthits);
- hnchits[ind][0].fill(nchits);
- hE[ind][0].fill(E);
- hEME[ind][0].fill(EMenergy);
- hnEMhits[ind][0].fill(nEMhits);
- hEvsEME[ind][0].fill(E,EMenergy);
- hnHADhits[ind][0].fill(nHADhits);
- hEvsHADE[ind][0].fill(E,HADenergy);
- hEvsCALE[ind][0].fill(E,CALenergy);
- }
- if(mc.getGeneratorStatus() > 0)
- {
- if(entries[ind][1] < maxentries)
- {
- entries[ind][1]++;
- hnthits[ind][1].fill(nthits);
- hnchits[ind][1].fill(nchits);
- hE[ind][1].fill(E);
- hEME[ind][1].fill(EMenergy);
- hnEMhits[ind][1].fill(nEMhits);
- hEvsEME[ind][1].fill(E,EMenergy);
- hnHADhits[ind][1].fill(nHADhits);
- hEvsHADE[ind][1].fill(E,HADenergy);
- hEvsCALE[ind][1].fill(E,CALenergy);
- }
- if(nEMhits+nHADhits > 10)
- {
- if(entries[ind][2] < maxentries)
- {
- entries[ind][2]++;
- hnthits[ind][2].fill(nthits);
- hnchits[ind][2].fill(nchits);
- hE[ind][2].fill(E);
- hEME[ind][2].fill(EMenergy);
- hnEMhits[ind][2].fill(nEMhits);
- hEvsEME[ind][2].fill(E,EMenergy);
- hnHADhits[ind][2].fill(nHADhits);
- hEvsHADE[ind][2].fill(E,HADenergy);
- hEvsCALE[ind][2].fill(E,CALenergy);
- }
- }
- }
- }
- ievt++;
-
- }
- private int getCreationIndex(MCParticle mcparticle, int demiseIndex)
- {
- int creationIndex = 0;
- boolean fromGenerator = mcparticle.getGeneratorStatus() > 0;
- boolean createdInSim = mcparticle.isCreatedInSimulation();
- if(fromGenerator)
- {
- int np = mcparticle.getNumberOfParents();
- boolean parentinsim = false;
- if(np == 1)
- {
- MCParticle p = mcparticle.getParent(0);
- if(p.isDecayedInTracker() || p.isDecayedInCalorimeter() ||
- p.hasLeftDetector() || p.isStopped() )
- {
- parentinsim = true;
- }
- }
- if( (demiseIndex == 0)&&(mcparticle.getGeneratorStatus() != 1) )
- {
-//
-// Simulator never saw this
-//
- creationIndex = 0;
- }
- else if(!createdInSim && !parentinsim)
- {
-//
-// Initial generator particle passed to simulator
-//
- creationIndex = 1;
- }
- else
- {
-//
-// Predecay
-//
- creationIndex = 2;
- }
- }
- else
- {
-//
-// Simulator created it on its own
-//
- if(mcparticle.isBackscatter())
- {
-//
-// backscatter
-//
- creationIndex = 3;
- }
- else if(mcparticle.vertexIsNotEndpointOfParent())
- {
-//
-// created without killing parent
-//
- creationIndex = 4;
- }
- else
- {
-//
-// result of parent decay or interaction
-//
- creationIndex = 5;
- }
- }
- return creationIndex;
- }
- int nparticledivisions = 12;
- int ncutdivisions = 3;
- String[] pdnames = {"gamma","electron","muon","pion","K charged","K0S","K0L","proton","neutron","ions","other+-","other0"};
- int[] pdpdg = {22,11,13,211,321,310,130,2212,2112,0,-1,-1};
- String[] cutnames = {"all","Genparts","nCalHits>10"};
- int[][] entries;
- int maxentries = 99999;
- IHistogram1D hPDGall;
- IHistogram1D hdIndex;
- IHistogram1D hcIndex;
- IHistogram1D hccIndex;
- ICloud1D[][] hnthits;
- ICloud1D[][] hnchits;
- ICloud1D[][] hE;
- ICloud1D[][] hEME;
- ICloud1D[][] hnEMhits;
- ICloud1D[][] hnHADhits;
- ICloud2D[][] hEvsEME;
- ICloud2D[][] hEvsHADE;
- ICloud2D[][] hEvsCALE;
- private int nTcollections = 3;
- private int nCcollections = 6;
- private int ng;
- private int np;
- private int ievt;
- private double sfEM;
- private double sfHAD;
-// private String[] TCname = {"VXDtCollection","CENtCollection"};
-// private String[] CCname = {"LUMcalCollection","EMcalCollection","HADcalCollection","MUONcalCollection"};
-// private String[] TCname = {"vxd","trk"};
-// private String[] CCname = {"lum","ecal","hcal","muon"};
- private String[] CCname = {"EcalBarrHits","EcalEndcapHits","HcalBarrHits","HcalEndcapHits",
- "MuonBarrHits","MuonEndcapHits"};
- private String[] TCname = {"VtxBarrHits","TrkBarrHits","TrkEndcapHits"};
-
-}
-class CCluster
-{
- int nhits;
- Vector hits;
- public CCluster()
- {
- nhits = 0;
- hits = new Vector();
- }
- public void addhit(SimCalorimeterHit hit)
- {
- hits.add(hit);
- nhits ++;
- }
- public int getNhits(){return nhits;}
- public SimCalorimeterHit[] getHits()
- {
- SimCalorimeterHit[] thits = new SimCalorimeterHit[nhits];
- for(int i=0;i<nhits;i++){thits[i] = (SimCalorimeterHit) hits.get(i);}
- return thits;
- }
-}
-class TCluster
-{
- int nhits;
- Vector hits;
- public TCluster()
- {
- nhits = 0;
- hits = new Vector();
- }
- public void addhit(SimTrackerHit hit)
- {
- hits.add(hit);
- nhits ++;
- }
- public int getNhits(){return nhits;}
- public SimTrackerHit[] getHits()
- {
- SimTrackerHit[] thits = new SimTrackerHit[nhits];
- for(int i=0;i<nhits;i++){thits[i] = (SimTrackerHit) hits.get(i);}
- return thits;
- }
-}
\ No newline at end of file