lcsim-contrib/src/main/java/org/lcsim/contrib/SteveMagill
diff -N ChCorrHistDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ChCorrHistDriver.java 12 Aug 2009 18:34:19 -0000 1.1
@@ -0,0 +1,242 @@
+package org.lcsim.contrib.SteveMagill;
+
+import java.util.*;
+import org.lcsim.util.aida.AIDA;
+import org.lcsim.event.EventHeader;
+import org.lcsim.util.Driver;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.ReconstructedParticle;
+
+public class ChCorrHistDriver extends Driver
+{
+ private double _escalib = 1.0;
+ private double _eccalib = 8000;
+ private String[] _scnames;
+ private String[] _chnames;
+ private String _rpname;
+ private AIDA aida = AIDA.defaultInstance();
+
+ public ChCorrHistDriver(double escalib, double eccalib)
+ {
+ _escalib = escalib;
+ _eccalib = eccalib;
+ }
+
+ public void process(EventHeader event)
+ {
+ // get energy of reco particle for this calibration
+ double parE = 0.;
+ List<ReconstructedParticle> recops = event.get(ReconstructedParticle.class,_rpname);
+ // don't do anything unless only 1 particle
+ if (recops.size() == 1)
+ {
+ for (ReconstructedParticle rp : recops)
+ {
+ parE = rp.getEnergy();
+ }
+
+ double scEsum = 0.;
+ double cscEsum = 0.;
+ double chEsum = 0.;
+ double cchEsum = 0.;
+ double scElhadl = 0.;
+
+ // first make scintillator plots
+ double nlhadl = 0;
+ double nhadhits = 0;
+ double nemhits = 0;
+ double nhadrat = 0.;
+ double nhademh = 0.;
+ for (int i=0; i<_scnames.length; i++)
+ {
+ List<CalorimeterHit> schits = event.get(CalorimeterHit.class,_scnames[i]);
+ for (CalorimeterHit schit : schits)
+ {
+ long shid = schit.getCellID();
+ double sce = _escalib*schit.getCorrectedEnergy();
+ scEsum += schit.getCorrectedEnergy();
+ cscEsum += _escalib*schit.getCorrectedEnergy(); // electron calibration
+ if (_scnames[i].equals("Edep_EcalBarrDigiHits"))
+ {
+ List<CalorimeterHit> emcbhits = event.get(CalorimeterHit.class,"Ceren_EcalBarrDigiHits");
+ for (CalorimeterHit echit : emcbhits)
+ {
+ long chid = echit.getCellID();
+ double ece = _eccalib*echit.getCorrectedEnergy();
+ if (shid == chid) aida.cloud1D("Cerenkov over Scin for cell").fill(ece/sce);
+ }
+ nemhits++;
+ int nlayer = schit.getMetaData().getIDDecoder().getLayer();
+ aida.cloud1D("Timing of EM scint hits").fill(schit.getTime());
+ if (nlayer<3) aida.cloud1D("First EM layers E dist Scin").fill(schit.getCorrectedEnergy());
+ }
+ if (_scnames[i].equals("Edep_EcalEndcapDigiHits"))
+ {
+ List<CalorimeterHit> emcechits = event.get(CalorimeterHit.class,"Ceren_EcalEndcapDigiHits");
+ for (CalorimeterHit echit : emcechits)
+ {
+ long chid = echit.getCellID();
+ double ece = _eccalib*echit.getCorrectedEnergy();
+ if (shid == chid) aida.cloud1D("Cerenkov over Scin for cell").fill(ece/sce);
+ }
+ nemhits++;
+ int nlayer = schit.getMetaData().getIDDecoder().getLayer();
+ aida.cloud1D("Timing of EM scint hits").fill(schit.getTime());
+ if (nlayer<3) aida.cloud1D("First EM layers E dist Scin").fill(schit.getCorrectedEnergy());
+ }
+ // add up energy and number of hits in last HCAL layer to check for containment
+ if (_scnames[i].equals("Edep_HcalBarrDigiHits"))
+ {
+ List<CalorimeterHit> hbhits = event.get(CalorimeterHit.class,"Ceren_HcalBarrDigiHits");
+ for (CalorimeterHit hbhit : hbhits)
+ {
+ long chid = hbhit.getCellID();
+ double hce = _eccalib*hbhit.getCorrectedEnergy();
+ if (shid == chid) aida.cloud1D("Cerenkov over Scin for cell").fill(hce/sce);
+ }
+ nhadhits++;
+ int nlayer = schit.getMetaData().getIDDecoder().getLayer();
+ aida.cloud1D("Timing of HAD scint hits").fill(schit.getTime());
+ if (nlayer<3) aida.cloud1D("First HAD layers E dist Scin").fill(schit.getCorrectedEnergy());
+ if (nlayer==16)
+ {
+ nlhadl++;
+ scElhadl += schit.getCorrectedEnergy();
+ }
+ }
+ if (_scnames[i].equals("Edep_HcalEndcapDigiHits"))
+ {
+ List<CalorimeterHit> hechits = event.get(CalorimeterHit.class,"Ceren_HcalEndcapDigiHits");
+ for (CalorimeterHit hechit : hechits)
+ {
+ long chid = hechit.getCellID();
+ double hce = _eccalib*hechit.getCorrectedEnergy();
+ if (shid == chid) aida.cloud1D("Cerenkov over Scin for cell").fill(hce/sce);
+ }
+ nhadhits++;
+ int nlayer = schit.getMetaData().getIDDecoder().getLayer();
+ aida.cloud1D("Timing of HAD scint hits").fill(schit.getTime());
+ if (nlayer<3) aida.cloud1D("First HAD layers E dist Scin").fill(schit.getCorrectedEnergy());
+ if (nlayer==16)
+ {
+ nlhadl++;
+ scElhadl += schit.getCorrectedEnergy();
+ }
+ }
+ }
+ }
+ aida.cloud1D("Scintillator Hit ESum").fill(scEsum);
+ double avghl = nhadhits/9.;
+ if (avghl>0)
+ {
+ nhadrat = nlhadl/avghl;
+ aida.cloud1D("Ratio of Nlhadl to Avg Nhadl").fill(nhadrat);
+ nhademh = nhadhits/(nhadhits+nemhits);
+ aida.cloud1D("Ratio of Nhad to Total hits").fill(nhademh);
+ }
+
+ // make cherenkov plots
+ for (int i=0; i<_chnames.length; i++)
+ {
+ List<CalorimeterHit> chhits = event.get(CalorimeterHit.class,_chnames[i]);
+ for (CalorimeterHit chhit : chhits)
+ {
+// aida.cloud1D("Timing of Cherenkov hits").fill(chhit.getTime());
+ chEsum += chhit.getCorrectedEnergy();
+ cchEsum += _eccalib*chhit.getCorrectedEnergy(); // electron calibration
+ // don't have any muons - kludge to get mip E from pions
+ if (_chnames[i].equals("Ceren_EcalBarrDigiHits") || _chnames[i].equals("Ceren_EcalEndcapDigiHits"))
+ {
+ int nlayer = chhit.getMetaData().getIDDecoder().getLayer();
+ aida.cloud1D("Timing of EM ceren hits").fill(chhit.getTime());
+ if (nlayer<3) aida.cloud1D("First EM layers E dist Ceren").fill(chhit.getCorrectedEnergy());
+ }
+ if (_chnames[i].equals("Ceren_HcalBarrDigiHits") || _chnames[i].equals("Ceren_HcalEndcapDigiHits"))
+ {
+ int nlayer = chhit.getMetaData().getIDDecoder().getLayer();
+ aida.cloud1D("Timing of HAD ceren hits").fill(chhit.getTime());
+ if (nlayer<3) aida.cloud1D("First HAD layers E dist Ceren").fill(chhit.getCorrectedEnergy());
+ }
+ }
+ }
+ aida.cloud1D("Cerenkov Hit ESum").fill(chEsum);
+
+ double mbEsum = 0.;
+ // get muon hit ESum for leakage check
+ List<CalorimeterHit> mbhits = event.get(CalorimeterHit.class,"MuonBarrHits");
+ aida.cloud1D("Number of muon hits").fill(mbhits.size());
+ for (CalorimeterHit mbhit : mbhits)
+ {
+ mbEsum += mbhit.getCorrectedEnergy();
+ }
+
+ aida.cloud1D("Muon hit ESum").fill(mbEsum);
+ aida.cloud2D("ScinEsum vs Muon Esum").fill(mbEsum,scEsum);
+ aida.cloud2D("Scin ESum vs Num Muon Hits").fill(mbhits.size(),scEsum);
+ // if (scElhadl/scEsum<0.05*(scEsum-12.) && scEsum>12. && scEsum<20. && scElhadl/scEsum<-0.05*(scEsum-20.))
+ // only plot for minimal leakage - low muon E, deposited ESum>.65, etc
+ if (cscEsum>0.)
+ {
+ aida.cloud1D("Total Esum C over S ratio").fill(cchEsum/cscEsum);
+ double rat = cchEsum/cscEsum;
+ if (mbhits.size()<5 && scEsum/parE>0.4)
+ {
+// aida.cloud1D("Ratio of last HCAL layer to total ESum").fill(scElhadl/scEsum);
+// aida.cloud2D("Rat lhadlE vs E").fill(scEsum,scElhadl/scEsum);
+ aida.cloud1D("eCalib Cher Hit ESum").fill(cchEsum);
+ aida.cloud1D("eCalib Scint Hit ESum").fill(cscEsum);
+ aida.cloud2D("Scint vs Cher Hit E").fill(chEsum,scEsum);
+ aida.cloud2D("Scint over E vs Cherenkov").fill(cchEsum,cscEsum/parE);
+ aida.cloud2D("Scint over E vs Cher over Scint").fill(rat,cscEsum/parE);
+ aida.cloud1D("C over S").fill(rat);
+ if (rat>0.0 && rat<0.1) aida.cloud1D("Scint over E bin .05").fill(cscEsum/parE);
+ if (rat>0.1 && rat<0.2) aida.cloud1D("Scint over E bin .15").fill(cscEsum/parE);
+ if (rat>0.2 && rat<0.3) aida.cloud1D("Scint over E bin .25").fill(cscEsum/parE);
+ if (rat>0.3 && rat<0.4) aida.cloud1D("Scint over E bin .35").fill(cscEsum/parE);
+ if (rat>0.4 && rat<0.5) aida.cloud1D("Scint over E bin .45").fill(cscEsum/parE);
+ if (rat>0.5 && rat<0.6) aida.cloud1D("Scint over E bin .55").fill(cscEsum/parE);
+ if (rat>0.6 && rat<0.7) aida.cloud1D("Scint over E bin .65").fill(cscEsum/parE);
+ if (rat>0.7 && rat<0.8) aida.cloud1D("Scint over E bin .75").fill(cscEsum/parE);
+ if (rat>0.8 && rat<0.9) aida.cloud1D("Scint over E bin .85").fill(cscEsum/parE);
+ if (rat>0.9 && rat<1.0) aida.cloud1D("Scint over E bin .95").fill(cscEsum/parE);
+ if (rat>1.0) aida.cloud1D("Scint over E bin gt 1.0").fill(cscEsum/parE);
+ }
+ aida.cloud1D("Average Calib Scint ESum").fill(cscEsum/0.711);
+ aida.cloud1D("Ch Corrected Scin ESum p2").fill(cscEsum/(.614-.303*rat+.684*rat*rat));
+ aida.cloud1D("Ch Corr Scin ESum p1").fill(cscEsum/(.345+.627*rat));
+ aida.cloud1D("Ch Corrected Scin ESum p3").fill(cscEsum/(.480+.628*rat-1.085*rat*rat+.975*rat*rat*rat));
+ }
+ if (mbhits.size()>4 && nhadrat>1.4 && nhademh>0.88) aida.cloud1D("E of Leak").fill(cscEsum);
+ if (mbhits.size()<5 && nhadrat<1.4 && nhademh<0.88) aida.cloud1D("E of No Leak").fill(cscEsum);
+ if (cscEsum/parE>0.1 && cscEsum/parE<0.7) // probably leakage, make plot
+ {
+ if (mbhits.size()>4) aida.cloud1D("Leakage Cands").fill(1.);
+ if (nhadrat>1.4) aida.cloud1D("Leakage Cands").fill(2.);
+ if (nhademh>0.88) aida.cloud1D("Leakage Cands").fill(3.);
+ }
+ if (cscEsum>14.)
+ {
+ if (mbhits.size()>4) aida.cloud1D("No Leakage Cands").fill(1.5);
+ if (nhadrat>1.4) aida.cloud1D("No Leakage Cands").fill(2.5);
+ if (nhademh>0.88) aida.cloud1D("No Leakage Cands").fill(3.5);
+ }
+ }
+ }
+
+ public void setScHitNames(String[] scnames)
+ {
+ _scnames = scnames;
+ }
+
+ public void setChHitNames(String[] chnames)
+ {
+ _chnames = chnames;
+ }
+
+ public void setMCRecoPName(String rpname)
+ {
+ _rpname = rpname;
+ }
+
+
+}