lcsim-contrib/src/main/java/org/lcsim/contrib/LGilbert
diff -N singlePlotsid03.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ singlePlotsid03.java 13 Aug 2010 23:42:37 -0000 1.1
@@ -0,0 +1,406 @@
+package org.lcsim.contrib.LGilbert;
+
+import hep.aida.*;
+import hep.physics.vec.Hep3Vector;
+
+import java.util.*;
+
+import org.lcsim.digisim.DigiPackageDriver;
+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.recon.cluster.util.QNeutralHadronClusterEnergyCalculator;
+import org.lcsim.recon.cluster.util.QPhotonClusterEnergyCalculator;
+import org.lcsim.util.aida.*;
+import org.lcsim.recon.util.CalInfoDriver;
+
+/* There is not currently a cut on pions, as I do not recall how to determine cos(theta) */
+
+public class singlePlotsid03 extends Driver
+{
+ private AIDA aida = AIDA.defaultInstance();
+
+ double phi;
+ boolean first;
+ double numberBarrelHits;
+ double numberEndcapHits;
+ double totalHits;
+ double meanEvents;
+ ArrayList<Double> numberBarrelEvents;
+ ArrayList<String> prefix;
+ ArrayList<String> id;
+
+ CalorimeterInformation ci;
+ ClusterEnergyCalculator cecnh;
+ ClusterEnergyCalculator cecph;
+ IAnalysisFactory analysisFactory;
+ ITreeFactory treeFactory;
+ ITree tree;
+ IDataPointSetFactory dpsf ;
+ IProfile1D pp1;
+ IProfile1D pp2;
+ IProfile1D pp3;
+ IProfile1D pp4;
+
+ IProfile1D nHenergy;
+ IProfile1D nbarHenergy;
+ IProfile1D piHenergy;
+ IProfile1D muHenergy;
+ IProfile1D k0lHenergy;
+
+ IProfile1D nHhits;
+ IProfile1D nbarHhits;
+ IProfile1D piHhits;
+ IProfile1D muHhits;
+ IProfile1D k0lHhits;
+
+ IProfile1D nAllenergy;
+ IProfile1D nbarAllenergy;
+ IProfile1D piAllenergy;
+ IProfile1D muAllenergy;
+ IProfile1D k0lAllenergy;
+
+ IProfile1D nAllhits;
+ IProfile1D nbarAllhits;
+ IProfile1D piAllhits;
+ IProfile1D muAllhits;
+ IProfile1D k0lAllhits;
+
+ public singlePlotsid03()
+ {
+ analysisFactory = IAnalysisFactory.create();
+ treeFactory = analysisFactory.createTreeFactory();
+ tree = treeFactory.create();
+ dpsf = analysisFactory.createDataPointSetFactory(tree);
+
+ add(new CalInfoDriver());
+ add(new DigiPackageDriver());
+
+ first = true;
+ numberBarrelEvents = new ArrayList();
+ prefix = new ArrayList();
+ id = new ArrayList();
+ phi = 0;
+ numberBarrelHits = 0;
+ numberEndcapHits = 0;
+ totalHits = 0;
+ meanEvents = 0;
+ cecnh = new QNeutralHadronClusterEnergyCalculator();
+ cecph = new QPhotonClusterEnergyCalculator();
+ }
+ protected void process(EventHeader event)
+ {
+ super.process(event);
+
+ System.out.println("Processing event" + event.getEventNumber());
+
+ prefix.add("EM & HAD Calorimeters/");
+ prefix.add("HAD Calorimeter Only/");
+ id.add("Other/");
+ id.add("KOLong/");
+ id.add("Mu/");
+ id.add("Neutron/");
+ id.add("Anti-neutron/");
+ id.add("Pion/");
+
+ if(first)
+ {
+ ci = CalorimeterInformation.instance();
+ first = false;
+ }
+
+ List<MCParticle> singleParticle = event.get(MCParticle.class, "MCParticle");
+ int numberMC = singleParticle.size();
+ MCParticle inputParticle = singleParticle.get(numberMC - 1);
+ double pdgID = inputParticle.getPDGID();
+ int idnumber = 0;
+
+ if (pdgID == 130)
+ {
+ idnumber = 1;
+ }
+ else if (pdgID == 13)
+ {
+ idnumber = 2;
+ }
+ else if (pdgID == 2112)
+ {
+ idnumber = 3;
+ }
+ else if (pdgID == -2112)
+ {
+ idnumber = 4;
+ }
+ else if (pdgID == 211 || pdgID == -211)
+ {
+ idnumber = 5;
+ }
+
+ double innerRadius = ci.getRMin(CalorimeterType.EM_BARREL);
+ double z = ci.getZMin(CalorimeterType.EM_ENDCAP);
+
+ Hep3Vector end = inputParticle.getEndPoint();
+ double magnitudeEndRadius = Math.sqrt(end.x() * end.x() + end.y() * end.y());
+ double magnitudeEndZ = Math.abs(end.z());
+ double px = inputParticle.getPX();
+ double py = inputParticle.getPY();
+ double energy = inputParticle.getEnergy();
+ 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.);
+ double mag = Math.sqrt(Math.pow(inputParticle.getPX(), 2) + Math.pow(inputParticle.getPY(), 2) + Math.pow(inputParticle.getPZ(), 2));
+ double cosTheta = Math.abs(inputParticle.getPZ() / mag);
+
+ BasicCluster allcal = new BasicCluster();
+ BasicCluster allcalEM = new BasicCluster();
+ BasicCluster allcalHAD = new BasicCluster();
+
+ if (magnitudeEndRadius > innerRadius || magnitudeEndZ > z || cosTheta > 0.7)
+ {
+ 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(prefix.get(0) + "Number of Barrel Hits").fill(numberBarrelHits);
+ aida.cloud1D(prefix.get(0) + "Number of Endcap Hits").fill(numberEndcapHits);
+ aida.cloud1D(prefix.get(0) + "Total Number of HCAL Hits").fill(totalHits);
+ aida.cloud1D(prefix.get(0) + "Energy").fill(energy);
+ aida.cloud1D(prefix.get(0) + "phi").fill(phi);
+
+ pp2 = aida.profile1D(prefix.get(0) + "HCal only Cal energy from calibration vs phi",15,0.,30.);
+ pp2.fill(phidf,cec.getEnergy(allcal));
+ pp1 = aida.profile1D(prefix.get(0) + "HCal only HADCal hits vs phi",15,0.,30.);
+ pp1.fill(phidf,totalHits);
+
+ aida.cloud1D(prefix.get(0) + id.get(idnumber) + "Number of Barrel Hits").fill(numberBarrelHits);
+ aida.cloud1D(prefix.get(0) + id.get(idnumber) + "Number of Endcap Hits").fill(numberEndcapHits);
+ aida.cloud1D(prefix.get(0) + id.get(idnumber) + "Total Number of HCAL Hits").fill(totalHits);
+ aida.cloud1D(prefix.get(0) + id.get(idnumber) + "Energy").fill(energy);
+ aida.cloud1D(prefix.get(0) + id.get(idnumber) + "phi").fill(phi);
+
+ if (idnumber == 1)
+ {
+ k0lAllenergy = aida.profile1D(prefix.get(1) + id.get(idnumber) + "K0L HCal only Cal energy from calibration vs phi",15,0.,30.);
+ k0lAllenergy.fill(phidf,cec.getEnergy(allcal));
+ k0lAllhits = aida.profile1D(prefix.get(1) + id.get(idnumber) + "K0L HCal only HADCal hits vs phi",15,0.,30.);
+ k0lAllhits.fill(phidf,totalHits);
+ }
+ else if (idnumber == 2)
+ {
+ muAllenergy = aida.profile1D(prefix.get(1) + id.get(idnumber) + "mu HCal only Cal energy from calibration vs phi",15,0.,30.);
+ muAllenergy.fill(phidf,cec.getEnergy(allcal));
+ muAllhits = aida.profile1D(prefix.get(1) + id.get(idnumber) + "mu HCal only HADCal hits vs phi",15,0.,30.);
+ muAllhits.fill(phidf,totalHits);
+ }
+ else if (idnumber == 3)
+ {
+ nAllenergy = aida.profile1D(prefix.get(1) + id.get(idnumber) + "n HCal only Cal energy from calibration vs phi",15,0.,30.);
+ nAllenergy.fill(phidf,cec.getEnergy(allcal));
+ nAllhits = aida.profile1D(prefix.get(1) + id.get(idnumber) + "n HCal only HADCal hits vs phi",15,0.,30.);
+ nAllhits.fill(phidf,totalHits);
+ }
+ else if (idnumber == 4)
+ {
+ nbarAllenergy = aida.profile1D(prefix.get(1) + id.get(idnumber) + "n bar HCal only Cal energy from calibration vs phi",15,0.,30.);
+ nbarAllenergy.fill(phidf,cec.getEnergy(allcal));
+ nbarAllhits = aida.profile1D(prefix.get(1) + id.get(idnumber) + "n bar HCal only HADCal hits vs phi",15,0.,30.);
+ nbarAllhits.fill(phidf,totalHits);
+ }
+ else if (idnumber == 5)
+ {
+ piAllenergy = aida.profile1D(prefix.get(1) + id.get(idnumber) + "pi HCal only Cal energy from calibration vs phi",15,0.,30.);
+ piAllenergy.fill(phidf,cec.getEnergy(allcal));
+ piAllhits = aida.profile1D(prefix.get(1) + id.get(idnumber) + "pi HCal only HADCal hits vs phi",15,0.,30.);
+ piAllhits.fill(phidf,totalHits);
+
+ // includes both pi + and pi -. Should I split those further? n & n bar produce slightly different results.
+ }
+ }
+ }
+
+ double hadRadius = ci.getRMin(CalorimeterType.HAD_BARREL);
+ double hadZ = ci.getZMin(CalorimeterType.HAD_ENDCAP);
+
+ if (magnitudeEndRadius > hadRadius || magnitudeEndZ > hadZ)
+ {
+ 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(prefix.get(1) + "Number of Barrel Hits").fill(numberBarrelHits);
+ aida.cloud1D(prefix.get(1) + "Number of Endcap Hits").fill(numberEndcapHits);
+ aida.cloud1D(prefix.get(1) + "Total Number of HCAL Hits").fill(totalHits);
+ aida.cloud1D(prefix.get(1) + "Energy").fill(energy);
+ aida.cloud1D(prefix.get(1) + "phi").fill(phi);
+
+ pp3 = aida.profile1D(prefix.get(0) + "HCal Only Cal energy from calibration vs phi",15,0.,30.);
+ pp3.fill(phidf,cec.getEnergy(allcal));
+ pp4 = aida.profile1D(prefix.get(0) + "HCal Only HADCal hits vs phi",15,0.,30.);
+ pp4.fill(phidf,totalHits);
+
+ aida.cloud1D(prefix.get(1) + id.get(idnumber) + "Number of Barrel Hits").fill(numberBarrelHits);
+ aida.cloud1D(prefix.get(1) + id.get(idnumber) + "Number of Endcap Hits").fill(numberEndcapHits);
+ aida.cloud1D(prefix.get(1) + id.get(idnumber) + "Total Number of HCAL Hits").fill(totalHits);
+ aida.cloud1D(prefix.get(1) + id.get(idnumber) + "Energy").fill(energy);
+ aida.cloud1D(prefix.get(1) + id.get(idnumber) + "phi").fill(phi);
+
+ if (idnumber == 1)
+ {
+ k0lHenergy = aida.profile1D(prefix.get(1) + id.get(idnumber) + "K0L Total Cal energy from calibration vs phi",15,0.,30.);
+ k0lHenergy.fill(phidf,cec.getEnergy(allcal));
+ k0lHhits = aida.profile1D(prefix.get(1) + id.get(idnumber) + "K0L Total HADCal hits vs phi",15,0.,30.);
+ k0lHhits.fill(phidf,totalHits);
+ }
+ else if (idnumber == 2)
+ {
+ muHenergy = aida.profile1D(prefix.get(1) + id.get(idnumber) + "mu Total Cal energy from calibration vs phi",15,0.,30.);
+ muHenergy.fill(phidf,cec.getEnergy(allcal));
+ muHhits = aida.profile1D(prefix.get(1) + id.get(idnumber) + "mu Total HADCal hits vs phi",15,0.,30.);
+ muHhits.fill(phidf,totalHits);
+ }
+ else if (idnumber == 3)
+ {
+ nHenergy = aida.profile1D(prefix.get(1) + id.get(idnumber) + "n Total Cal energy from calibration vs phi",15,0.,30.);
+ nHenergy.fill(phidf,cec.getEnergy(allcal));
+ nHhits = aida.profile1D(prefix.get(1) + id.get(idnumber) + "n Total HADCal hits vs phi",15,0.,30.);
+ nHhits.fill(phidf,totalHits);
+ }
+ else if (idnumber == 4)
+ {
+ nbarHenergy = aida.profile1D(prefix.get(1) + id.get(idnumber) + "n bar Total Cal energy from calibration vs phi",15,0.,30.);
+ nbarHenergy.fill(phidf,cec.getEnergy(allcal));
+ nbarHhits = aida.profile1D(prefix.get(1) + id.get(idnumber) + "n bar Total HADCal hits vs phi",15,0.,30.);
+ nbarHhits.fill(phidf,totalHits);
+ }
+ else if (idnumber == 5)
+ {
+ piHenergy = aida.profile1D(prefix.get(1) + id.get(idnumber) + "pi Total Cal energy from calibration vs phi",15,0.,30.);
+ piHenergy.fill(phidf,cec.getEnergy(allcal));
+ piHhits = aida.profile1D(prefix.get(1) + id.get(idnumber) + "pi Total HADCal hits vs phi",15,0.,30.);
+ piHhits.fill(phidf,totalHits);
+ }
+ }
+ }
+ }
+ protected void suspend()
+ {
+ IDataPointSet dps1 = dpsf.create("<#HCAL hits20> vs fphi","<hcal hits>20 vs phi",2);
+ IDataPointSet dps2 = dpsf.create("<calE20> vs fphi","<calE>20 vs phi",2);
+ IAxis axis = pp1.axis();
+ int nb = axis.bins();
+ for(int i=0;i<nb;i++)
+ {
+ dps1.addPoint();
+ IDataPoint dp = dps1.point(i);
+ dp.coordinate(0).setValue(axis.binCenter(i));
+ dp.coordinate(1).setValue(pp1.binHeight(i));
+ double err = pp1.binRms(i)/Math.sqrt(pp1.binEntries(i));
+ dp.coordinate(1).setErrorPlus(err);
+ dp.coordinate(1).setErrorMinus(err);
+ }
+ axis = pp2.axis();
+ nb = axis.bins();
+ for(int i=0;i<nb;i++)
+ {
+ dps2.addPoint();
+ IDataPoint dp = dps2.point(i);
+ dp.coordinate(0).setValue(axis.binCenter(i));
+ dp.coordinate(1).setValue(pp2.binHeight(i));
+ double err = pp2.binRms(i)/Math.sqrt(pp2.binEntries(i));
+ dp.coordinate(1).setErrorPlus(err);
+ dp.coordinate(1).setErrorMinus(err);
+ }
+
+ double sumEvents = 0;
+ for (int i = 0; i < numberBarrelEvents.size(); i++)
+ {
+ sumEvents = sumEvents + numberBarrelEvents.get(i);
+ }
+
+ meanEvents = sumEvents / numberBarrelEvents.size();
+ System.out.println("Mean number of events: " + meanEvents);
+ }
+}
lcsim-contrib/src/main/java/org/lcsim/contrib/LGilbert
diff -N singlePlotsid02.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ singlePlotsid02.java 13 Aug 2010 23:42:37 -0000 1.1
@@ -0,0 +1,404 @@
+package org.lcsim.contrib.LGilbert;
+import hep.aida.*;
+import hep.physics.vec.Hep3Vector;
+
+import java.util.*;
+
+import org.lcsim.contrib.Cassell.recon.RemoveHcalModuleNonProjBorderHits;
+import org.lcsim.digisim.DigiPackageDriver;
+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.recon.cluster.util.QNeutralHadronClusterEnergyCalculator;
+import org.lcsim.recon.cluster.util.QPhotonClusterEnergyCalculator;
+import org.lcsim.util.aida.*;
+import org.lcsim.recon.util.CalInfoDriver;
+
+/* There is not currently a cut on pions, as I do not recall how to determine cos(theta) */
+
+public class singlePlotsid02 extends Driver
+{
+ private AIDA aida = AIDA.defaultInstance();
+
+ double phi;
+ boolean first;
+ double numberBarrelHits;
+ double numberEndcapHits;
+ double totalHits;
+ double meanEvents;
+ ArrayList<Double> numberBarrelEvents;
+ ArrayList<String> prefix;
+ ArrayList<String> id;
+
+ CalorimeterInformation ci;
+ ClusterEnergyCalculator cecnh;
+ ClusterEnergyCalculator cecph;
+ IAnalysisFactory analysisFactory;
+ ITreeFactory treeFactory;
+ ITree tree;
+ IDataPointSetFactory dpsf ;
+ IProfile1D pp1;
+ IProfile1D pp2;
+ IProfile1D pp3;
+ IProfile1D pp4;
+
+ IProfile1D nHenergy;
+ IProfile1D nbarHenergy;
+ IProfile1D piHenergy;
+ IProfile1D muHenergy;
+ IProfile1D k0lHenergy;
+
+ IProfile1D nHhits;
+ IProfile1D nbarHhits;
+ IProfile1D piHhits;
+ IProfile1D muHhits;
+ IProfile1D k0lHhits;
+
+ IProfile1D nAllenergy;
+ IProfile1D nbarAllenergy;
+ IProfile1D piAllenergy;
+ IProfile1D muAllenergy;
+ IProfile1D k0lAllenergy;
+
+ IProfile1D nAllhits;
+ IProfile1D nbarAllhits;
+ IProfile1D piAllhits;
+ IProfile1D muAllhits;
+ IProfile1D k0lAllhits;
+
+ public singlePlotsid02()
+ {
+ analysisFactory = IAnalysisFactory.create();
+ treeFactory = analysisFactory.createTreeFactory();
+ tree = treeFactory.create();
+ dpsf = analysisFactory.createDataPointSetFactory(tree);
+
+ add(new CalInfoDriver());
+ add(new DigiPackageDriver());
+ add(new RemoveHcalModuleNonProjBorderHits(10));
+
+ first = true;
+ numberBarrelEvents = new ArrayList();
+ phi = 0;
+ numberBarrelHits = 0;
+ numberEndcapHits = 0;
+ totalHits = 0;
+ meanEvents = 0;
+ cecnh = new QNeutralHadronClusterEnergyCalculator();
+ cecph = new QPhotonClusterEnergyCalculator();
+ }
+ protected void process(EventHeader event)
+ {
+ super.process(event);
+
+
+ prefix.add("EM & HAD Calorimeters/");
+ prefix.add("HAD Calorimeter Only/");
+ id.add("Other/");
+ id.add("KOLong/");
+ id.add("Mu/");
+ id.add("Neutron/");
+ id.add("Anti-neutron/");
+ id.add("Pion/");
+
+ if(first)
+ {
+ ci = CalorimeterInformation.instance();
+ first = false;
+ }
+
+ List<MCParticle> singleParticle = event.get(MCParticle.class, "MCParticle");
+ int numberMC = singleParticle.size();
+ MCParticle inputParticle = singleParticle.get(numberMC - 1);
+ double pdgID = inputParticle.getPDGID();
+ int idnumber = 0;
+
+ if (pdgID == 130)
+ {
+ idnumber = 1;
+ }
+ else if (pdgID == 13)
+ {
+ idnumber = 2;
+ }
+ else if (pdgID == 2112)
+ {
+ idnumber = 3;
+ }
+ else if (pdgID == -2112)
+ {
+ idnumber = 4;
+ }
+ else if (pdgID == 211)
+ {
+ idnumber = 5;
+ }
+ else if (pdgID == -211)
+ {
+ idnumber = 6;
+ }
+
+ double innerRadius = ci.getRMin(CalorimeterType.EM_BARREL);
+ double z = ci.getZMin(CalorimeterType.EM_ENDCAP);
+
+ Hep3Vector end = inputParticle.getEndPoint();
+ double magnitudeEndRadius = Math.sqrt(end.x() * end.x() + end.y() * end.y());
+ double magnitudeEndZ = Math.abs(end.z());
+ double px = inputParticle.getPX();
+ double py = inputParticle.getPY();
+ double energy = inputParticle.getEnergy();
+ 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(prefix.get(1) + "Number of Barrel Hits").fill(numberBarrelHits);
+ aida.cloud1D(prefix.get(1) + "Number of Endcap Hits").fill(numberEndcapHits);
+ aida.cloud1D(prefix.get(1) + "Total Number of HCAL Hits").fill(totalHits);
+ aida.cloud1D(prefix.get(1) + "Energy").fill(energy);
+ aida.cloud1D(prefix.get(1) + "phi").fill(phi);
+
+ pp2 = aida.profile1D(prefix.get(1) + "HCal only Cal energy from calibration vs phi",15,0.,30.);
+ pp2.fill(phidf,cec.getEnergy(allcal));
+ pp1 = aida.profile1D(prefix.get(1) + "HCal only HADCal hits vs phi",15,0.,30.);
+ pp1.fill(phidf,totalHits);
+
+ aida.cloud1D(prefix.get(1) + id.get(idnumber) + "Number of Barrel Hits").fill(numberBarrelHits);
+ aida.cloud1D(prefix.get(1) + id.get(idnumber) + "Number of Endcap Hits").fill(numberEndcapHits);
+ aida.cloud1D(prefix.get(1) + id.get(idnumber) + "Total Number of HCAL Hits").fill(totalHits);
+ aida.cloud1D(prefix.get(1) + id.get(idnumber) + "Energy").fill(energy);
+ aida.cloud1D(prefix.get(1) + id.get(idnumber) + "phi").fill(phi);
+
+ if (idnumber == 1)
+ {
+ k0lAllenergy = aida.profile1D(prefix.get(1) + id.get(idnumber) + "K0L HCal only Cal energy from calibration vs phi",15,0.,30.);
+ k0lAllenergy.fill(phidf,cec.getEnergy(allcal));
+ k0lAllhits = aida.profile1D(prefix.get(1) + id.get(idnumber) + "K0L HCal only HADCal hits vs phi",15,0.,30.);
+ k0lAllhits.fill(phidf,totalHits);
+ }
+ else if (idnumber == 2)
+ {
+ muAllenergy = aida.profile1D(prefix.get(1) + id.get(idnumber) + "mu HCal only Cal energy from calibration vs phi",15,0.,30.);
+ muAllenergy.fill(phidf,cec.getEnergy(allcal));
+ muAllhits = aida.profile1D(prefix.get(1) + id.get(idnumber) + "mu HCal only HADCal hits vs phi",15,0.,30.);
+ muAllhits.fill(phidf,totalHits);
+ }
+ else if (idnumber == 3)
+ {
+ nAllenergy = aida.profile1D(prefix.get(1) + id.get(idnumber) + "n HCal only Cal energy from calibration vs phi",15,0.,30.);
+ nAllenergy.fill(phidf,cec.getEnergy(allcal));
+ nAllhits = aida.profile1D(prefix.get(1) + id.get(idnumber) + "n HCal only HADCal hits vs phi",15,0.,30.);
+ nAllhits.fill(phidf,totalHits);
+ }
+ else if (idnumber == 4)
+ {
+ nbarAllenergy = aida.profile1D(prefix.get(1) + id.get(idnumber) + "n bar HCal only Cal energy from calibration vs phi",15,0.,30.);
+ nbarAllenergy.fill(phidf,cec.getEnergy(allcal));
+ nbarAllhits = aida.profile1D(prefix.get(1) + id.get(idnumber) + "n bar HCal only HADCal hits vs phi",15,0.,30.);
+ nbarAllhits.fill(phidf,totalHits);
+ }
+ else if (idnumber == 5)
+ {
+ piAllenergy = aida.profile1D(prefix.get(1) + id.get(idnumber) + "pi HCal only Cal energy from calibration vs phi",15,0.,30.);
+ piAllenergy.fill(phidf,cec.getEnergy(allcal));
+ piAllhits = aida.profile1D(prefix.get(1) + id.get(idnumber) + "pi HCal only HADCal hits vs phi",15,0.,30.);
+ piAllhits.fill(phidf,totalHits);
+ }
+ }
+ }
+
+ double hadRadius = ci.getRMin(CalorimeterType.HAD_BARREL);
+ double hadZ = ci.getZMin(CalorimeterType.HAD_ENDCAP);
+
+ if (magnitudeEndRadius > hadRadius || magnitudeEndZ > hadZ)
+ {
+ 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(prefix.get(1) + "Number of Barrel Hits").fill(numberBarrelHits);
+ aida.cloud1D(prefix.get(0) + "Number of Endcap Hits").fill(numberEndcapHits);
+ aida.cloud1D(prefix.get(0) + "Total Number of HCAL Hits").fill(totalHits);
+ aida.cloud1D(prefix.get(0) + "Energy").fill(energy);
+ aida.cloud1D(prefix.get(0) + "phi").fill(phi);
+
+ pp3 = aida.profile1D(prefix.get(0) + "HCal Only Cal energy from calibration vs phi",15,0.,30.);
+ pp3.fill(phidf,cec.getEnergy(allcal));
+ pp4 = aida.profile1D(prefix.get(0) + "HCal Only HADCal hits vs phi",15,0.,30.);
+ pp4.fill(phidf,totalHits);
+
+ aida.cloud1D(prefix.get(0) + id.get(idnumber) + "Number of Barrel Hits").fill(numberBarrelHits);
+ aida.cloud1D(prefix.get(0) + id.get(idnumber) + "Number of Endcap Hits").fill(numberEndcapHits);
+ aida.cloud1D(prefix.get(0) + id.get(idnumber) + "Total Number of HCAL Hits").fill(totalHits);
+ aida.cloud1D(prefix.get(0) + id.get(idnumber) + "Energy").fill(energy);
+ aida.cloud1D(prefix.get(0) + id.get(idnumber) + "phi").fill(phi);
+
+ if (idnumber == 1)
+ {
+ k0lHenergy = aida.profile1D(prefix.get(0) + id.get(idnumber) + "K0L Total Cal energy from calibration vs phi",15,0.,30.);
+ k0lHenergy.fill(phidf,cec.getEnergy(allcal));
+ k0lHhits = aida.profile1D(prefix.get(0) + id.get(idnumber) + "K0L Total HADCal hits vs phi",15,0.,30.);
+ k0lHhits.fill(phidf,totalHits);
+ }
+ else if (idnumber == 2)
+ {
+ muHenergy = aida.profile1D(prefix.get(0) + id.get(idnumber) + "mu Total Cal energy from calibration vs phi",15,0.,30.);
+ muHenergy.fill(phidf,cec.getEnergy(allcal));
+ muHhits = aida.profile1D(prefix.get(0) + id.get(idnumber) + "mu Total HADCal hits vs phi",15,0.,30.);
+ muHhits.fill(phidf,totalHits);
+ }
+ else if (idnumber == 3)
+ {
+ nHenergy = aida.profile1D(prefix.get(0) + id.get(idnumber) + "n Total Cal energy from calibration vs phi",15,0.,30.);
+ nHenergy.fill(phidf,cec.getEnergy(allcal));
+ nHhits = aida.profile1D(prefix.get(0) + id.get(idnumber) + "n Total HADCal hits vs phi",15,0.,30.);
+ nHhits.fill(phidf,totalHits);
+ }
+ else if (idnumber == 4)
+ {
+ nbarHenergy = aida.profile1D(prefix.get(0) + id.get(idnumber) + "n bar Total Cal energy from calibration vs phi",15,0.,30.);
+ nbarHenergy.fill(phidf,cec.getEnergy(allcal));
+ nbarHhits = aida.profile1D(prefix.get(0) + id.get(idnumber) + "n bar Total HADCal hits vs phi",15,0.,30.);
+ nbarHhits.fill(phidf,totalHits);
+ }
+ else if (idnumber == 5)
+ {
+ piHenergy = aida.profile1D(prefix.get(0) + id.get(idnumber) + "pi Total Cal energy from calibration vs phi",15,0.,30.);
+ piHenergy.fill(phidf,cec.getEnergy(allcal));
+ piHhits = aida.profile1D(prefix.get(0) + id.get(idnumber) + "pi Total HADCal hits vs phi",15,0.,30.);
+ piHhits.fill(phidf,totalHits);
+ }
+ }
+ }
+ }
+ protected void suspend()
+ {
+ IDataPointSet dps1 = dpsf.create("<#HCAL hits20> vs fphi","<hcal hits>20 vs phi",2);
+ IDataPointSet dps2 = dpsf.create("<calE20> vs fphi","<calE>20 vs phi",2);
+ IAxis axis = pp1.axis();
+ int nb = axis.bins();
+ for(int i=0;i<nb;i++)
+ {
+ dps1.addPoint();
+ IDataPoint dp = dps1.point(i);
+ dp.coordinate(0).setValue(axis.binCenter(i));
+ dp.coordinate(1).setValue(pp1.binHeight(i));
+ double err = pp1.binRms(i)/Math.sqrt(pp1.binEntries(i));
+ dp.coordinate(1).setErrorPlus(err);
+ dp.coordinate(1).setErrorMinus(err);
+ }
+ axis = pp2.axis();
+ nb = axis.bins();
+ for(int i=0;i<nb;i++)
+ {
+ dps2.addPoint();
+ IDataPoint dp = dps2.point(i);
+ dp.coordinate(0).setValue(axis.binCenter(i));
+ dp.coordinate(1).setValue(pp2.binHeight(i));
+ double err = pp2.binRms(i)/Math.sqrt(pp2.binEntries(i));
+ dp.coordinate(1).setErrorPlus(err);
+ dp.coordinate(1).setErrorMinus(err);
+ }
+
+ double sumEvents = 0;
+ for (int i = 0; i < numberBarrelEvents.size(); i++)
+ {
+ sumEvents = sumEvents + numberBarrelEvents.get(i);
+ }
+
+ meanEvents = sumEvents / numberBarrelEvents.size();
+ System.out.println("Mean number of events: " + meanEvents);
+ }
+}
