lcsim-contrib/src/main/java/org/lcsim/contrib/SteveMagill
diff -N ChCalHitDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ChCalHitDriver.java 12 Aug 2009 18:31:10 -0000 1.1
@@ -0,0 +1,168 @@
+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;
+
+public class ChCalHitDriver extends Driver
+{
+ public ChCalHitDriver() {
+ }
+
+ private String[] _scnames;
+ private String[] _chnames;
+ private AIDA aida = AIDA.defaultInstance();
+
+ public void process(EventHeader event)
+ {
+
+ 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]);
+ // CalorimeterIDDecoder caldecoder = (CalorimeterIDDecoder) event.getMetaData(schits).getIDDecoder();
+ for (CalorimeterHit schit : schits)
+ {
+ scEsum += schit.getCorrectedEnergy();
+ cscEsum += 1.013*schit.getCorrectedEnergy(); // electron calibration
+ // don't have any muons - kludge to get mip E from pions
+ if (_scnames[i].equals("Edep_EcalBarrDigiHits") || _scnames[i].equals("Edep_EcalEndcapDigiHits"))
+ {
+ nemhits++;
+ int nlayer = schit.getMetaData().getIDDecoder().getLayer();
+ // int nlayer = caldecoder.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") || _scnames[i].equals("Edep_HcalEndcapDigiHits"))
+ {
+ nhadhits++;
+ int nlayer = schit.getMetaData().getIDDecoder().getLayer();
+ // int nlayer = caldecoder.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==8)
+ {
+ 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]);
+ // CalorimeterIDDecoder caldecoder = (CalorimeterIDDecoder) event.getMetaData(chhits).getIDDecoder();
+ for (CalorimeterHit chhit : chhits)
+ {
+// aida.cloud1D("Timing of Cherenkov hits").fill(chhit.getTime());
+ chEsum += chhit.getCorrectedEnergy();
+ cchEsum += 8326*chhit.getCorrectedEnergy(); // electron calibration
+ // don't have any muons - kludge to get mip E from pions
+ if (_chnames[i].equals("ChEBHits") || _chnames[i].equals("ChEECHits"))
+ {
+ int nlayer = chhit.getMetaData().getIDDecoder().getLayer();
+ // int nlayer = caldecoder.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("ChHBHits") || _chnames[i].equals("ChHECHits"))
+ {
+ int nlayer = chhit.getMetaData().getIDDecoder().getLayer();
+ // int nlayer = caldecoder.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.))
+// if (mbEsum<0.5 && scEsum>13. && scElhadl/scEsum<0.04*(scEsum-13.) && cscEsum>0.)
+// {
+ aida.cloud1D("Ratio of last HCAL layer to total ESum").fill(scElhadl/scEsum);
+ aida.cloud2D("Rat lhadlE vs E").fill(scEsum,scElhadl/scEsum);
+ double rat = cchEsum/cscEsum;
+ 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/20.);
+// aida.cloud2D("Scint over E vs Cher over Scint").fill(rat,cscEsum/20.);
+// if (rat>0.35 && rat<0.45) aida.cloud1D("Scint over E bin .4").fill(cscEsum/20.);
+// if (rat>0.45 && rat<0.55) aida.cloud1D("Scint over E bin .5").fill(cscEsum/20.);
+// if (rat>0.55 && rat<0.65) aida.cloud1D("Scint over E bin .6").fill(cscEsum/20.);
+// if (rat>0.65 && rat<0.75) aida.cloud1D("Scint over E bin .7").fill(cscEsum/20.);
+// if (rat>0.75 && rat<0.85) aida.cloud1D("Scint over E bin .8").fill(cscEsum/20.);
+// if (rat>0.85 && rat<0.95) aida.cloud1D("Scint over E bin .9").fill(cscEsum/20.);
+// if (rat>0.95 && rat<1.05) aida.cloud1D("Scint over E bin 1.0").fill(cscEsum/20.);
+// if (rat>1.05 && rat<1.15) aida.cloud1D("Scint over E bin 1.1").fill(cscEsum/20.);
+// if (rat>1.15 && rat<1.25) aida.cloud1D("Scint over E bin 1.2").fill(cscEsum/20.);
+ aida.cloud1D("Shifted Calib Scint ESum").fill(cscEsum/0.68);
+ aida.cloud1D("Ch Corr Scin ESum p2").fill(cscEsum/(.765-.719*rat+.794*rat*rat));
+ aida.cloud1D("Ch Corr Scin ESum p1").fill(cscEsum/(.395+.392*rat));
+ aida.cloud1D("Ch Corrected Scin ESum p3").fill(cscEsum/(.595+.529*rat-.524*rat*rat+.277*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>2. && cscEsum<14.) // 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;
+ }
+
+}