lcsim-contrib/src/main/java/org/lcsim/contrib/LGilbert
diff -N QqbarAnalysisDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ QqbarAnalysisDriver.java 29 Jul 2010 00:43:35 -0000 1.1
@@ -0,0 +1,335 @@
+package org.lcsim.contrib.LGilbert;
+
+import org.lcsim.event.ReconstructedParticle;
+import hep.physics.vec.*;
+import hep.aida.*;
+import java.util.*;
+import org.lcsim.recon.analysis.RMS90Calculator;
+
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.util.*;
+import org.lcsim.recon.util.CalorimeterInformation;
+import org.lcsim.geometry.Calorimeter.CalorimeterType;
+import org.lcsim.recon.cluster.util.BasicCluster;
+import org.lcsim.recon.cluster.util.ClusterEnergyCalculator;
+import org.lcsim.util.aida.*;
+
+
+public class QqbarAnalysisDriver extends Driver
+{
+ private AIDA aida = AIDA.defaultInstance();
+ int ievt;
+ int nmax = 1000000;
+ double[] phil = {3.333,6.667,10.,13.333,16.667,20.,23.333,26.667,30.};
+ double[] phiv;
+ double phidf;
+ double ctcut = 0.8;
+ double ZE;
+ RMS90Calculator calc;
+ IAnalysisFactory analysisFactory;
+ ITreeFactory treeFactory;
+ ITree tree;
+ IDataPointSetFactory dpsf ;
+ String pre;
+ String pre2;
+
+ double phi;
+ boolean first;
+ double numberBarrelHits;
+ double numberEndcapHits;
+ double totalHits;
+ double meanEvents;
+ ArrayList<Double> numberBarrelEvents;
+
+ CalorimeterInformation ci;
+ ClusterEnergyCalculator cecnh;
+ ClusterEnergyCalculator cecph;
+ IProfile1D pp1;
+ IProfile1D pp2;
+ IProfile1D pp3;
+ IProfile1D pp4;
+//
+// Keep track of the different cm Energies for plots
+//
+ int nE;
+ int[] Eevt;
+
+ public QqbarAnalysisDriver()
+ {
+ ievt = 0;
+ analysisFactory = IAnalysisFactory.create();
+ treeFactory = analysisFactory.createTreeFactory();
+ tree = treeFactory.create();
+ dpsf = analysisFactory.createDataPointSetFactory(tree);
+ calc = new RMS90Calculator();
+ first = true;
+ nE = 0;
+ Eevt = new int[10];
+ phiv = new double[phil.length];
+ double lv = 0.;
+
+ for(int i=0;i<phil.length;i++)
+ {
+ phiv[i] = lv + (phil[i] - lv)/2.;
+ lv = phil[i];
+ }
+ }
+ protected void process(EventHeader event)
+ {
+ List<MCParticle> mcl = event.get(MCParticle.class,"MCParticle");
+ double Zmass = 0.;
+ ZE = 0.;
+ double ct = 0.;
+
+ if(first)
+ {
+ ci = CalorimeterInformation.instance();
+ int El = (int) (ZE + .5);
+ Eevt[0] = El;
+ nE = 1;
+ pre = "cmE="+Eevt[0]+"/";
+ System.out.println("New CM Energy = "+El+"GeV at event "+ievt);
+ first = false;
+ }
+
+ for(MCParticle p:mcl)
+ {
+ int id = Math.abs(p.getPDGID());
+ if( (id == 1 )||(id == 2 )||(id == 3 ) )
+ {
+ List<MCParticle> parents = p.getParents();
+ for(MCParticle parent:parents)
+ {
+ if(parent.getPDGID() == 23 )
+ {
+ Zmass = parent.getMass();
+ ZE = parent.getEnergy();
+ Hep3Vector pp = p.getMomentum();
+ ct = Math.abs(pp.z()/pp.magnitude());
+ double phi = Math.atan2(pp.y(),pp.x());
+ if(phi < 0.)phi += 2.*Math.PI;
+ double phideg = phi*180./Math.PI;
+ int n12 = (int) (phideg/30.);
+ phidf = phideg - 30.*n12;
+ }
+ }
+ }
+ }
+
+ System.out.println("Processing event "+ievt);
+ if(ct > ctcut)
+ {
+ ievt++;
+ return;
+ }
+
+ super.process(event);
+
+ int phib = 0;
+
+ for(int i=0;i<phil.length;i++)
+ {
+ if(phidf < phil[i])break;
+ phib++;
+ }
+ pre2 = "phibin "+phib+"/";
+//
+// See if cmE changed
+//
+ int Ec = (int) (ZE + .5);
+ boolean match = false;
+ for(int i=0;i<nE;i++)
+ {
+ if(Ec == Eevt[i])
+ {
+ match = true;
+ break;
+ }
+ }
+ if(!match)
+ {
+ System.out.println("New CM Energy = "+Ec+"GeV at event "+ievt);
+ Eevt[nE] = Ec;
+// tree.mkdirs("cmE="+Eevt[nE]+"/");
+ pre = "cmE="+Eevt[nE]+"/";
+ nE++;
+ }
+
+ List<ReconstructedParticle> rl = event.get(ReconstructedParticle.class,"ReconstructedParticles");
+ double Etot = 0.;
+ double px = 0.;
+ double py = 0.;
+ double pz = 0.;
+ for(ReconstructedParticle p:rl)
+ {
+ Etot += p.getEnergy();
+ Hep3Vector mom = p.getMomentum();
+ px += mom.x();
+ py += mom.y();
+ pz += mom.z();
+ }
+
+ double innerRadius = ci.getRMin(CalorimeterType.EM_BARREL);
+ double z = ci.getZMin(CalorimeterType.EM_ENDCAP);
+
+ List<MCParticle> singleParticle = event.get(MCParticle.class, "MCParticle");
+ int numberMC = singleParticle.size();
+ MCParticle inputParticle = singleParticle.get(numberMC - 1);
+ double pdgID = inputParticle.getPDGID();
+
+ Hep3Vector end = inputParticle.getEndPoint();
+ double magnitudeEndRadius = Math.sqrt(end.x() * end.x() + end.y() * end.y());
+ double magnitudeEndZ = Math.abs(end.z());
+ double phi = Math.atan2(py, px);
+
+ if(phi < 0.)
+ {
+ phi += 2.*Math.PI;
+ }
+
+ double phideg = phi*180./Math.PI;
+ int n12 = (int) (phideg/30.);
+ double phidf = phideg - 30.*n12;
+ int phibin = (int) (phidf/2.);
+
+ BasicCluster allcal = new BasicCluster();
+ BasicCluster allcalEM = new BasicCluster();
+ BasicCluster allcalHAD = new BasicCluster();
+
+ if (magnitudeEndRadius > innerRadius || magnitudeEndZ > z)
+ {
+ List<SimCalorimeterHit> barrelHits = event.get(SimCalorimeterHit.class, ci.getDigiCollectionName(CalorimeterType.HAD_BARREL));
+ List<SimCalorimeterHit> endcapHits = event.get(SimCalorimeterHit.class, ci.getDigiCollectionName(CalorimeterType.HAD_ENDCAP));
+
+ numberBarrelHits = barrelHits.size();
+ numberEndcapHits = endcapHits.size();
+ totalHits = numberBarrelHits + numberEndcapHits;
+ numberBarrelEvents.add(numberBarrelHits);
+
+ double totCalEsf = 0.;
+ double totCalEEMsf = 0.;
+ double totCalEHADsf = 0.;
+ int totHADhits = 0;
+ int totEMhits = 0;
+ String[] collnames = {ci.getDigiCollectionName("EM_BARREL"), ci.getDigiCollectionName("EM_ENDCAP"), ci.getDigiCollectionName("HAD_BARREL"), ci.getDigiCollectionName("HAD_ENDCAP")};
+
+ for(int i=0;i<collnames.length;i++)
+ {
+ for(CalorimeterHit h:event.get(CalorimeterHit.class,collnames[i]))
+ {
+ allcal.addHit(h);
+ totCalEsf += h.getCorrectedEnergy();
+
+ if(i < 2)
+ {
+ allcalEM.addHit(h);
+ totCalEEMsf += h.getCorrectedEnergy();
+ totEMhits++;
+ }
+ else
+ {
+ allcalHAD.addHit(h);
+ totCalEHADsf += h.getCorrectedEnergy();
+ totHADhits++;
+ }
+ }
+
+ ClusterEnergyCalculator cec = cecnh;
+
+
+ aida.cloud1D("Number of HCAL Barrel Hits").fill(numberBarrelHits);
+ aida.cloud1D("Number of HCAL Endcap Hits").fill(numberEndcapHits);
+ aida.cloud1D("Total Number of HCAL Hits").fill(totalHits);
+ aida.cloud1D("Energy").fill(Etot);
+ aida.cloud1D("phi").fill(phi);
+
+ pp2 = aida.profile1D("HCal only Cal energy vs phi",15,0.,30.);
+ pp2.fill(phidf,cec.getEnergy(allcal));
+ pp1 = aida.profile1D("HCal only HADCal hits vs phi",15,0.,30.);
+ pp1.fill(phidf,totalHits);
+
+ }
+
+ double M = Math.sqrt(Etot*Etot - px*px - py*py - pz*pz);
+ double dE = Etot - ZE;
+ double dM = M - Zmass;
+ aida.cloud1D(pre+"Event Energy",nmax).fill(Etot);
+ aida.cloud1D(pre+"Event Mass",nmax).fill(M);
+ aida.cloud1D(pre+"delta Event Energy",nmax).fill(dE);
+ aida.cloud1D(pre+"delta Event Mass",nmax).fill(dM);
+ aida.cloud1D(pre+pre2+"Event Energy",nmax).fill(Etot);
+ aida.cloud1D(pre+pre2+"Event Mass",nmax).fill(M);
+ aida.cloud1D(pre+pre2+"delta Event Energy",nmax).fill(dE);
+ aida.cloud1D(pre+pre2+"delta Event Mass",nmax).fill(dM);
+ ievt++;
+ }
+ }
+
+
+ protected void suspend()
+ {
+ IDataPointSet dpsE;
+ IDataPointSet[] dpsE2 = new IDataPointSet[nE];
+ IDataPointSet dpsM;
+ IDataPointSet[] dpsM2 = new IDataPointSet[nE];
+ dpsE = dpsf.create("dE over E vs cmE","dEoE vs E",2);
+ dpsM = dpsf.create("mean90 vs cmE","m90 vs E",2);
+ for(int iE=0;iE<nE;iE++)
+ {
+ dpsE2[iE] = dpsf.create("E="+Eevt[iE]+":dE over E vs phibin","Ecm="+Eevt[iE]+":dEoE vs phi",2);
+ dpsM2[iE] = dpsf.create("E="+Eevt[iE]+":mean90 vs phibin","Ecm="+Eevt[iE]+":m90 vs phi",2);
+ ZE = Eevt[iE];
+ String rname = "cmE="+Eevt[iE]+"/Event Energy";
+ ICloud1D thisc = aida.cloud1D(rname);
+ double r90 = calc.calculateRMS90(thisc);
+ double m90 = calc.getMEAN90();
+ int ent = thisc.entries();
+ double m = thisc.mean();
+ double r = thisc.rms();
+ double dEoE = r90/(m90);
+ double err2 = (1.1*r90/Math.sqrt(.9*ent));
+ double err = (1.1*r90/Math.sqrt(1.8*ent))/(m90);
+ dpsE.addPoint();
+ dpsM.addPoint();
+ IDataPoint dp = dpsE.point(iE);
+ IDataPoint dp2 = dpsM.point(iE);
+ dp.coordinate(0).setValue(ZE);
+ dp2.coordinate(0).setValue(ZE);
+ dp.coordinate(1).setValue(dEoE);
+ dp2.coordinate(1).setValue(m90);
+ dp.coordinate(1).setErrorPlus(err);
+ dp.coordinate(1).setErrorMinus(err);
+ dp2.coordinate(1).setErrorPlus(err2);
+ dp2.coordinate(1).setErrorMinus(err2);
+ for(int i=0;i<phil.length;i++)
+ {
+ rname = "cmE="+Eevt[iE]+"/phibin "+i+"/Event Energy";
+ thisc = aida.cloud1D(rname);
+ r90 = calc.calculateRMS90(thisc);
+ m90 = calc.getMEAN90();
+ ent = thisc.entries();
+ m = thisc.mean();
+ r = thisc.rms();
+ dEoE = r90/(m90);
+ err2 = (1.1*r90/Math.sqrt(.9*ent));
+ err = (1.1*r90/Math.sqrt(1.8*ent))/(m90);
+ dpsE2[iE].addPoint();
+ dpsM2[iE].addPoint();
+ dp = dpsE2[iE].point(i);
+ dp2 = dpsM2[iE].point(i);
+ dp.coordinate(0).setValue(phiv[i]);
+ dp2.coordinate(0).setValue(phiv[i]);
+ dp.coordinate(1).setValue(dEoE);
+ dp2.coordinate(1).setValue(m90);
+ dp.coordinate(1).setErrorPlus(err);
+ dp.coordinate(1).setErrorMinus(err);
+ dp2.coordinate(1).setErrorPlus(err2);
+ dp2.coordinate(1).setErrorMinus(err2);
+ }
+ }
+ super.suspend();
+ }
+}
\ No newline at end of file