lcsim/src/org/lcsim/contrib/CalAnal
diff -N AnalyzeNeutralsZpole.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ AnalyzeNeutralsZpole.java 27 Nov 2006 17:10:31 -0000 1.1
@@ -0,0 +1,314 @@
+import java.util.*;
+import hep.aida.*;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.EventHeader;
+import org.lcsim.recon.cluster.cheat.*;
+import org.lcsim.recon.cluster.util.*;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.geometry.IDDecoder;
+import hep.physics.vec.*;
+import org.lcsim.event.util.*;
+import org.lcsim.event.util.MCParticleClassifier.MCPClass;
+
+/**
+ * Reconstruction: Clusters
+ */
+public class AnalyzeNeutralsZpole extends Driver
+{
+ private AIDA aida = AIDA.defaultInstance();
+ IDDecoder _decoder;
+ String EMBname = "EcalBarrHitsClusters";
+ String HADBname = "HcalBarrHitsClusters";
+ String EMBhitsname = "EcalBarrHits";
+ String HADBhitsname = "HcalBarrHits";
+ String EMEname = "EcalEndcapHitsClusters";
+ String HADEname = "HcalEndcapHitsClusters";
+ String EMEhitsname = "EcalEndcapHits";
+ String HADEhitsname = "HcalEndcapHits";
+ CheatClusterer _clusterer;
+ int cloudmax = 200000;
+ int ntypes;
+ String[] partnames = {"neutron","k0L","nbar","gamma"};
+ int[] pdg = {2112,130,-2112,22};
+ int ncalculators = 1;
+ double[] emeas;
+ DetailedNeutralHadronClusterEnergyCalculator calc;
+ PhotonClusterEnergyCalculator scalc;
+ MCPClass fs;
+ MCParticleClassifier cl;
+ double emtimecut = 100.;
+ double hadtimecut = 100.;
+ double photonnorm = 1.011531;
+ double nhadnorm;
+
+ Map<MCParticle,CheatCluster> embclmap;
+ Map<MCParticle,CheatCluster> hadbclmap;
+ Map<MCParticle,CheatCluster> emeclmap;
+ Map<MCParticle,CheatCluster> hadeclmap;
+ List<MCParticle> plist;
+ double tanthcut = .25;
+ String hEcut = "Eventcuts/";
+ boolean first;
+ public AnalyzeNeutralsZpole()
+ {
+
+ CreateFinalStateMCParticleList fsDriver = new CreateFinalStateMCParticleList("Gen");
+ add(fsDriver);
+ ntypes = 4;
+ _clusterer = new CheatClusterer();
+ emeas = new double[ncalculators];
+ first = true;
+ }
+ protected void process(EventHeader event)
+ {
+ super.process(event);
+//
+// Cannot instantiate ClusterEnergyCalculators until we have an event
+//
+ if(first)
+ {
+ first = false;
+ scalc = new PhotonClusterEnergyCalculator();
+ calc = new DetailedNeutralHadronClusterEnergyCalculator();
+ }
+//
+// Make event cut
+//
+ double all_all_energy = 0.;
+ double all_gamma_energy = 0.;
+ double all_nh_energy = 0.;
+ double all_neutral_energy = 0.;
+ double cut_all_energy = 0.;
+ double cut_gamma_energy = 0.;
+ double cut_nh_energy = 0.;
+ double cut_neutral_energy = 0.;
+ String fslistname = "GenFinalStateParticles";
+ plist = event.get(MCParticle.class,fslistname);
+ for(MCParticle p:plist)
+ {
+ all_all_energy += p.getEnergy();
+ int pindex = -1;
+ int pdgid = p.getPDGID();
+ for(int i=0;i<ntypes;i++)
+ {
+ if(pdgid == pdg[i])pindex = i;
+ }
+ if(pindex > -1)
+ {
+ all_neutral_energy += p.getEnergy();
+ if(pindex == ntypes-1)
+ {
+ all_gamma_energy += p.getEnergy();
+ }
+ else
+ {
+ all_nh_energy += p.getEnergy();
+ }
+ }
+ Hep3Vector pp = p.getMomentum();
+ double pperp = Math.sqrt(pp.x()*pp.x() + pp.y()*pp.y());
+ double tanth = 999.;
+ if(Math.abs(pp.z()) > 0.001)tanth = pperp/Math.abs(pp.z());
+ if(tanth > tanthcut)
+ {
+ cut_all_energy += p.getEnergy();
+ if(pindex > -1)
+ {
+ cut_neutral_energy += p.getEnergy();
+ if(pindex == ntypes-1)
+ {
+ cut_gamma_energy += p.getEnergy();
+ }
+ else
+ {
+ cut_nh_energy += p.getEnergy();
+ }
+ }
+ }
+ }
+ aida.cloud1D(hEcut+"Sum of all FS energies per event").fill(all_all_energy);
+ aida.cloud1D(hEcut+"Sum of all FS neutral energies per event").fill(all_neutral_energy);
+ aida.cloud1D(hEcut+"Sum of all FS gamma energies per event").fill(all_gamma_energy);
+ aida.cloud1D(hEcut+"Sum of all FS neut hadron energies per event").fill(all_nh_energy);
+ aida.cloud1D(hEcut+"Sum of cut FS energies per event").fill(cut_all_energy);
+ aida.cloud1D(hEcut+"Sum of cut FS neutral energies per event").fill(cut_neutral_energy);
+ aida.cloud1D(hEcut+"Sum of cut FS gamma energies per event").fill(cut_gamma_energy);
+ aida.cloud1D(hEcut+"Sum of cut FS neut hadron energies per event").fill(cut_nh_energy);
+ if(cut_all_energy > 89.)
+ {
+ aida.cloud1D(hEcut+"Cut event Sum of cut FS neutral energies per event").fill(cut_neutral_energy);
+ aida.cloud1D(hEcut+"Cut event Sum of cut FS gamma energies per event").fill(cut_gamma_energy);
+ aida.cloud1D(hEcut+"Cut event Sum of cut FS neut hadron energies per event").fill(cut_nh_energy);
+
+
+
+ double missed_neutral_energy = 0.;
+ double missed_nh_energy = 0.;
+ double missed_gamma_energy = 0.;
+ double[] measured_neutral_energy = new double[ncalculators];
+ double[] measured_nh_energy = new double[ncalculators];
+ double[] measured_gamma_energy = new double[ncalculators];
+ for(int i=0;i<ncalculators;i++)
+ {
+ measured_neutral_energy[i] = 0.;
+ measured_nh_energy[i] = 0.;
+ measured_gamma_energy[i] = 0.;
+ }
+//
+// Use the cluster cheater on the EMBarrel and HADBarrel collections
+//
+ List<List<SimCalorimeterHit>> collections = event.get(SimCalorimeterHit.class);
+ for (List<SimCalorimeterHit> collection : collections)
+ {
+ String name = event.getMetaData(collection).getName();
+ if( (name.compareTo(EMBhitsname) == 0) || (name.compareTo(HADBhitsname) == 0) ||
+ (name.compareTo(EMEhitsname) == 0) || (name.compareTo(HADEhitsname) == 0))
+ {
+ Map<MCParticle,CheatCluster> result = _clusterer.findClusters(collection);
+ if (result.size() > 0) event.put(name+"Clusters",new ArrayList(result.values()));
+ if( (name.compareTo(EMBhitsname) == 0))embclmap = result;
+ if( (name.compareTo(HADBhitsname) == 0))hadbclmap = result;
+ if( (name.compareTo(EMEhitsname) == 0))emeclmap = result;
+ if( (name.compareTo(HADEhitsname) == 0))hadeclmap = result;
+ }
+ }
+ for(MCParticle p:plist)
+ {
+ int pindex = -1;
+ int pdgid = p.getPDGID();
+ for(int i=0;i<ntypes;i++)
+ {
+ if(pdgid == pdg[i])pindex = i;
+ }
+ if(pindex > -1)
+ {
+//
+// Make angle cut
+//
+ Hep3Vector pp = p.getMomentum();
+ double pperp = Math.sqrt(pp.x()*pp.x() + pp.y()*pp.y());
+ double tanth = 999.;
+ if(Math.abs(pp.z()) > 0.001)tanth = pperp/Math.abs(pp.z());
+ if(tanth > tanthcut)
+ {
+//
+// Particle survived cuts, look for hits
+//
+ aida.cloud1D(partnames[pindex]+"/Energy",cloudmax).fill(p.getEnergy());
+ aida.cloud1D("Neutrals/Energy",cloudmax).fill(p.getEnergy());
+ aida.cloud1D(partnames[pindex]+"/wted Energy",cloudmax).fill(p.getEnergy(),p.getEnergy());
+ aida.cloud1D("Neutrals/wted Energy",cloudmax).fill(p.getEnergy(),p.getEnergy());
+ if(pindex < ntypes-1)
+ {
+ aida.cloud1D("NHadron/Energy",cloudmax).fill(p.getEnergy());
+ aida.cloud1D("NHadron/wted Energy",cloudmax).fill(p.getEnergy(),p.getEnergy());
+ }
+ List<MCParticle> plist = new ArrayList<MCParticle>();
+ plist.add(p);
+ int next = 0;
+ while(next < plist.size())
+ {
+ List<MCParticle> dlist = plist.get(next).getDaughters();
+ for(MCParticle d:dlist)
+ {
+ plist.add(d);
+ }
+ next++;
+ }
+ BasicCluster cl = new BasicCluster();
+ for(MCParticle ap:plist)
+ {
+ if(embclmap.containsKey(ap))cl.addCluster(embclmap.get(ap));
+ if(hadbclmap.containsKey(ap))cl.addCluster(hadbclmap.get(ap));
+ if(emeclmap.containsKey(ap))cl.addCluster(emeclmap.get(ap));
+ if(hadeclmap.containsKey(ap))cl.addCluster(hadeclmap.get(ap));
+ }
+ if( cl.getSize() == 0 )
+ {
+ for(int i=0;i<ncalculators;i++)
+ {
+ emeas[i] = 0.;
+ }
+ missed_neutral_energy += p.getEnergy();
+ aida.cloud1D("Neutrals/Missed energy").fill(p.getEnergy());
+ aida.cloud1D(partnames[pindex]+"/Missed energy").fill(p.getEnergy());
+ if(pindex == ntypes-1)missed_gamma_energy += p.getEnergy();
+ else
+ {
+ missed_nh_energy += p.getEnergy();
+ aida.cloud1D("NHadron/Missed energy").fill(p.getEnergy());
+ }
+ }
+ else
+ {
+ if(pindex == ntypes-1)
+ {
+ emeas[0] = scalc.getEnergy(cl);
+ for(int i=0;i<ncalculators;i++)
+ {
+ measured_gamma_energy[i] += emeas[i];
+ measured_neutral_energy[i] += emeas[i];
+ double dE = emeas[i] - p.getEnergy();
+ aida.cloud1D(partnames[pindex]+"/dE"+i,cloudmax).fill(dE);
+ aida.cloud1D(partnames[pindex]+"/dE over E"+i,cloudmax).fill(dE/p.getEnergy());
+ aida.cloud1D(partnames[pindex]+"/dE over rootE"+i,cloudmax).fill(dE/Math.sqrt(p.getEnergy()));
+ aida.cloud1D("Neutrals/dE"+i,cloudmax).fill(dE);
+ aida.cloud1D("Neutrals/dE over E"+i,cloudmax).fill(dE/p.getEnergy());
+ aida.cloud1D("Neutrals/dE over rootE"+i,cloudmax).fill(dE/Math.sqrt(p.getEnergy()));
+ }
+ }
+ else
+ {
+ emeas[0] = calc.getEnergy(cl);
+ for(int i=0;i<ncalculators;i++)
+ {
+ measured_neutral_energy[i] += emeas[i];
+ measured_nh_energy[i] += emeas[i];
+ double dE = emeas[i] - p.getEnergy();
+ aida.cloud1D(partnames[pindex]+"/dE"+i,cloudmax).fill(dE);
+ aida.cloud1D(partnames[pindex]+"/dE over E"+i,cloudmax).fill(dE/p.getEnergy());
+ aida.cloud1D(partnames[pindex]+"/dE over rootE"+i,cloudmax).fill(dE/Math.sqrt(p.getEnergy()));
+ aida.cloud1D("Neutrals/dE"+i,cloudmax).fill(dE);
+ aida.cloud1D("Neutrals/dE over E"+i,cloudmax).fill(dE/p.getEnergy());
+ aida.cloud1D("Neutrals/dE over rootE"+i,cloudmax).fill(dE/Math.sqrt(p.getEnergy()));
+ aida.cloud1D("NHadron/dE"+i,cloudmax).fill(dE);
+ aida.cloud1D("NHadron/dE over E"+i,cloudmax).fill(dE/p.getEnergy());
+ aida.cloud1D("NHadron/dE over rootE"+i,cloudmax).fill(dE/Math.sqrt(p.getEnergy()));
+ }
+ }
+ }
+ }
+ }
+ }
+ String hp = "Event/";
+ for(int i=0;i<ncalculators;i++)
+ {
+ double dh = measured_neutral_energy[i] - cut_neutral_energy;
+ double dhoe = dh/cut_neutral_energy;
+ double dhore = dh/Math.sqrt(cut_neutral_energy);
+ aida.cloud1D(hp+"Neutral dE"+i).fill(dh);
+ aida.cloud1D(hp+"Neutral dE"+i+" over E").fill(dhoe);
+ aida.cloud1D(hp+"Neutral dE"+i+" over rootE").fill(dhore);
+ dh = measured_nh_energy[i] - cut_nh_energy;
+ dhoe = dh/cut_nh_energy;
+ dhore = dh/Math.sqrt(cut_nh_energy);
+ aida.cloud1D(hp+"NHadron dE"+i).fill(dh);
+ aida.cloud1D(hp+"NHadron dE"+i+" over E").fill(dhoe);
+ aida.cloud1D(hp+"NHadron dE"+i+" over rootE").fill(dhore);
+ dh = measured_gamma_energy[i] - cut_gamma_energy;
+ dhoe = dh/cut_gamma_energy;
+ dhore = dh/Math.sqrt(cut_gamma_energy);
+ aida.cloud1D(hp+"Gamma dE"+i).fill(dh);
+ aida.cloud1D(hp+"Gamma dE"+i+" over E").fill(dhoe);
+ aida.cloud1D(hp+"Gamma dE"+i+" over rootE").fill(dhore);
+ }
+ aida.cloud1D(hp+"Neutral missed Energy").fill(missed_neutral_energy);
+ aida.cloud1D(hp+"NHadron missed Energy").fill(missed_nh_energy);
+ aida.cloud1D(hp+"Gamma missed Energy").fill(missed_gamma_energy);
+ }
+ }
+}
