CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HcalSimHitStudy.cc
Go to the documentation of this file.
3 
5 
7 
8  g4Label = ps.getUntrackedParameter<std::string>("moduleLabel","g4SimHits");
9  hcalHits = ps.getUntrackedParameter<std::string>("HitCollection","HcalHits");
10  outFile_ = ps.getUntrackedParameter<std::string>("outputFile", "hcHit.root");
11  verbose_ = ps.getUntrackedParameter<bool>("Verbose", false);
12  checkHit_= true;
13 
14  tok_hits_ = consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label,hcalHits));
15 
16  edm::LogInfo("HcalSim") << "Module Label: " << g4Label << " Hits: "
17  << hcalHits << " / "<< checkHit_
18  << " Output: " << outFile_;
19 
20 }
21 
23 
25 {
26  ib.setCurrentFolder("HcalHitsV/HcalSimHitsTask");
27 
28  //Histograms for Hits
29  if (checkHit_) {
30  meAllNHit_ = ib.book1D("Hit01","Number of Hits in HCal",20000,0.,20000.);
31  meBadDetHit_= ib.book1D("Hit02","Hits with wrong Det", 100,0.,100.);
32  meBadSubHit_= ib.book1D("Hit03","Hits with wrong Subdet",100,0.,100.);
33  meBadIdHit_ = ib.book1D("Hit04","Hits with wrong ID", 100,0.,100.);
34  meHBNHit_ = ib.book1D("Hit05","Number of Hits in HB",20000,0.,20000.);
35  meHENHit_ = ib.book1D("Hit06","Number of Hits in HE",10000,0.,10000.);
36  meHONHit_ = ib.book1D("Hit07","Number of Hits in HO",10000,0.,10000.);
37  meHFNHit_ = ib.book1D("Hit08","Number of Hits in HF",10000,0.,10000.);
38  meDetectHit_= ib.book1D("Hit09","Detector ID", 50,0.,50.);
39  meSubdetHit_= ib.book1D("Hit10","Subdetectors in HCal", 50,0.,50.);
40  meDepthHit_ = ib.book1D("Hit11","Depths in HCal", 20,0.,20.);
41  meEtaHit_ = ib.book1D("Hit12","Eta in HCal", 101,-50.5,50.5);
42  //KC: There are different phi segmentation schemes, this plot uses wider bins to represent the most sparse segmentation
43  mePhiHit_ = ib.book1D("Hit13","Phi in HCal (HB,HO)", 72,0.5,72.5);
44  mePhiHitb_ = ib.book1D("Hit13b","Phi in HCal (HE,HF)", 72,0.5,72.5);
45  meEnergyHit_= ib.book1D("Hit14","Energy in HCal", 2000,0.,20.);
46  meTimeHit_ = ib.book1D("Hit15","Time in HCal", 528,0.,528.);
47  meTimeWHit_ = ib.book1D("Hit16","Time in HCal (E wtd)", 528,0.,528.);
48  meHBDepHit_ = ib.book1D("Hit17","Depths in HB", 20,0.,20.);
49  meHEDepHit_ = ib.book1D("Hit18","Depths in HE", 20,0.,20.);
50  meHODepHit_ = ib.book1D("Hit19","Depths in HO", 20,0.,20.);
51  meHFDepHit_ = ib.book1D("Hit20","Depths in HF", 20,0.,20.);
52  meHBEtaHit_ = ib.book1D("Hit21","Eta in HB", 101,-50.5,50.5);
53  meHEEtaHit_ = ib.book1D("Hit22","Eta in HE", 101,-50.5,50.5);
54  meHOEtaHit_ = ib.book1D("Hit23","Eta in HO", 101,-50.5,50.5);
55  meHFEtaHit_ = ib.book1D("Hit24","Eta in HF", 101,-50.5,50.5);
56  meHBPhiHit_ = ib.book1D("Hit25","Phi in HB", 72,0.5,72.5);
57  meHEPhiHit_ = ib.book1D("Hit26","Phi in HE", 72,0.5,72.5);
58  meHOPhiHit_ = ib.book1D("Hit27","Phi in HO", 72,0.5,72.5);
59  meHFPhiHit_ = ib.book1D("Hit28","Phi in HF", 72,0.5,72.5);
60  meHBEneHit_ = ib.book1D("Hit29","Energy in HB", 2000,0.,20.);
61  meHEEneHit_ = ib.book1D("Hit30","Energy in HE", 500,0.,5.);
62  meHOEneHit_ = ib.book1D("Hit31","Energy in HO", 500,0.,5.);
63  meHFEneHit_ = ib.book1D("Hit32","Energy in HF", 1000,0.5,1000.5);
64  meHBTimHit_ = ib.book1D("Hit33","Time in HB", 528,0.,528.);
65  meHETimHit_ = ib.book1D("Hit34","Time in HE", 528,0.,528.);
66  meHOTimHit_ = ib.book1D("Hit35","Time in HO", 528,0.,528.);
67  meHFTimHit_ = ib.book1D("Hit36","Time in HF", 528,0.,528.);
68  //These are the zoomed in energy ranges
69  meHBEneHit2_ = ib.book1D("Hit37","Energy in HB 2", 100,0.,0.0001);
70  meHEEneHit2_ = ib.book1D("Hit38","Energy in HE 2", 100,0.,0.0001);
71  meHOEneHit2_ = ib.book1D("Hit39","Energy in HO 2", 100,0.,0.0001);
72  meHFEneHit2_ = ib.book1D("Hit40","Energy in HF 2", 100,0.5,100.5);
73  meHBL10Ene_ = ib.book1D("Hit41","Log10Energy in HB", 140, -10., 4. );
74  meHEL10Ene_ = ib.book1D("Hit42","Log10Energy in HE", 140, -10., 4. );
75  meHFL10Ene_ = ib.book1D("Hit43","Log10Energy in HF", 50, -1., 4. );
76  meHOL10Ene_ = ib.book1D("Hit44","Log10Energy in HO", 140, -10., 4. );
77  meHBL10EneP_ = ib.bookProfile("Hit45","Log10Energy in HB vs Hit contribution", 140, -10., 4., 100, 0., 1. );
78  meHEL10EneP_ = ib.bookProfile("Hit46","Log10Energy in HE vs Hit contribution", 140, -10., 4., 100, 0., 1. );
79  meHFL10EneP_ = ib.bookProfile("Hit47","Log10Energy in HF vs Hit contribution", 140, -10., 4., 100, 0., 1. );
80  meHOL10EneP_ = ib.bookProfile("Hit48","Log10Energy in HO vs Hit contribution", 140, -10., 4., 100, 0., 1. );
81  }
82 
83 }
84 
85 
86 /*void HcalSimHitStudy::endJob() {
87  if (dbe_ && outFile_.size() > 0) dbe_->save(outFile_);
88 }*/
89 
91 
92  LogDebug("HcalSim") << "Run = " << e.id().run() << " Event = "
93  << e.id().event();
94 
95  std::vector<PCaloHit> caloHits;
97 
98  bool getHits = false;
99  if (checkHit_) {
100  e.getByToken(tok_hits_,hitsHcal);
101  if (hitsHcal.isValid()) getHits = true;
102  }
103 
104  LogDebug("HcalSim") << "HcalValidation: Input flags Hits " << getHits;
105 
106  if (getHits) {
107  caloHits.insert(caloHits.end(),hitsHcal->begin(),hitsHcal->end());
108  LogDebug("HcalSim") << "HcalValidation: Hit buffer "
109  << caloHits.size();
110  analyzeHits (caloHits);
111  }
112 }
113 
114 void HcalSimHitStudy::analyzeHits (std::vector<PCaloHit>& hits) {
115 
116  int nHit = hits.size();
117  int nHB=0, nHE=0, nHO=0, nHF=0, nBad1=0, nBad2=0, nBad=0;
118  std::vector<double> encontHB(140, 0.);
119  std::vector<double> encontHE(140, 0.);
120  std::vector<double> encontHF(140, 0.);
121  std::vector<double> encontHO(140, 0.);
122  double entotHB = 0, entotHE = 0, entotHF = 0, entotHO = 0;
123 
124  for (int i=0; i<nHit; i++) {
125  double energy = hits[i].energy();
126  double log10en = log10(energy);
127  int log10i = int( (log10en+10.)*10. );
128  double time = hits[i].time();
129  unsigned int id_ = hits[i].id();
130  HcalDetId id = HcalDetId(id_);
131  int det = id.det();
132  int subdet = id.subdet();
133  int depth = id.depth();
134  int eta = id.ieta();
135  int phi = id.iphi();
136  LogDebug("HcalSim") << "Hit[" << i << "] ID " << std::hex << id_
137  << std::dec << " Det " << det << " Sub "
138  << subdet << " depth " << depth << " Eta " << eta
139  << " Phi " << phi << " E " << energy << " time "
140  << time;
141  if (det == 4) { // Check DetId.h
142  if (subdet == static_cast<int>(HcalBarrel)) nHB++;
143  else if (subdet == static_cast<int>(HcalEndcap)) nHE++;
144  else if (subdet == static_cast<int>(HcalOuter)) nHO++;
145  else if (subdet == static_cast<int>(HcalForward)) nHF++;
146  else { nBad++; nBad2++;}
147  } else { nBad++; nBad1++;}
148 
149  meDetectHit_->Fill(double(det));
150  if (det == 4) {
151  meSubdetHit_->Fill(double(subdet));
152  meDepthHit_->Fill(double(depth));
153  meEtaHit_->Fill(double(eta));
154 
155  //We will group the phi plots by HB,HO and HE,HF since these groups share similar segmentation schemes
156  if (subdet == static_cast<int>(HcalBarrel)) mePhiHit_->Fill(double(phi));
157  else if (subdet == static_cast<int>(HcalEndcap)) mePhiHitb_->Fill(double(phi));
158  else if (subdet == static_cast<int>(HcalOuter)) mePhiHit_->Fill(double(phi));
159  else if (subdet == static_cast<int>(HcalForward)) mePhiHitb_->Fill(double(phi));
160 
161 
162  //KC: HF energy is in photoelectrons rather than eV, so it will not be included in total HCal energy
163  if (subdet != static_cast<int>(HcalForward)){
164  meEnergyHit_->Fill(energy);
165 
166  //Since the HF energy is a different scale it does not make sense to include it in the Energy Weighted Plot
167  meTimeWHit_->Fill(double(time),energy);
168  }
169  meTimeHit_->Fill(time);
170 
171  if (subdet == static_cast<int>(HcalBarrel)) {
172  meHBDepHit_->Fill(double(depth));
173  meHBEtaHit_->Fill(double(eta));
174  meHBPhiHit_->Fill(double(phi));
175  meHBEneHit_->Fill(energy);
176  meHBEneHit2_->Fill(energy);
177  meHBTimHit_->Fill(time);
178  meHBL10Ene_->Fill(log10en);
179  if( log10i >=0 && log10i < 140 ) encontHB[log10i] += energy;
180  entotHB += energy;
181  } else if (subdet == static_cast<int>(HcalEndcap)) {
182  meHEDepHit_->Fill(double(depth));
183  meHEEtaHit_->Fill(double(eta));
184  meHEPhiHit_->Fill(double(phi));
185  meHEEneHit_->Fill(energy);
186  meHEEneHit2_->Fill(energy);
187  meHETimHit_->Fill(time);
188  meHEL10Ene_->Fill(log10en);
189  if( log10i >=0 && log10i < 140 ) encontHE[log10i] += energy;
190  entotHE += energy;
191  } else if (subdet == static_cast<int>(HcalOuter)) {
192  meHODepHit_->Fill(double(depth));
193  meHOEtaHit_->Fill(double(eta));
194  meHOPhiHit_->Fill(double(phi));
195  meHOEneHit_->Fill(energy);
196  meHOEneHit2_->Fill(energy);
197  meHOTimHit_->Fill(time);
198  meHOL10Ene_->Fill(log10en);
199  if( log10i >=0 && log10i < 140 ) encontHO[log10i] += energy;
200  entotHO += energy;
201  } else if (subdet == static_cast<int>(HcalForward)) {
202  meHFDepHit_->Fill(double(depth));
203  meHFEtaHit_->Fill(double(eta));
204  meHFPhiHit_->Fill(double(phi));
205  meHFEneHit_->Fill(energy);
206  meHFEneHit2_->Fill(energy);
207  meHFTimHit_->Fill(time);
208  meHFL10Ene_->Fill(log10en);
209  if( log10i >=0 && log10i < 140 ) encontHF[log10i] += energy;
210  entotHF += energy;
211  }
212  }
213  }
214  if( entotHB != 0 ) for( int i=0; i<140; i++ ) meHBL10EneP_->Fill( -10.+(float(i)+0.5)/10., encontHB[i]/entotHB );
215  if( entotHE != 0 ) for( int i=0; i<140; i++ ) meHEL10EneP_->Fill( -10.+(float(i)+0.5)/10., encontHE[i]/entotHE );
216  if( entotHF != 0 ) for( int i=0; i<140; i++ ) meHFL10EneP_->Fill( -10.+(float(i)+0.5)/10., encontHF[i]/entotHF );
217  if( entotHO != 0 ) for( int i=0; i<140; i++ ) meHOL10EneP_->Fill( -10.+(float(i)+0.5)/10., encontHO[i]/entotHO );
218 
219  meAllNHit_->Fill(double(nHit));
220  meBadDetHit_->Fill(double(nBad1));
221  meBadSubHit_->Fill(double(nBad2));
222  meBadIdHit_->Fill(double(nBad));
223  meHBNHit_->Fill(double(nHB));
224  meHENHit_->Fill(double(nHE));
225  meHONHit_->Fill(double(nHO));
226  meHFNHit_->Fill(double(nHF));
227 
228  LogDebug("HcalSim") << "HcalSimHitStudy::analyzeHits: HB " << nHB
229  << " HE " << nHE << " HO " << nHO << " HF " << nHF
230  << " Bad " << nBad << " All " << nHit;
231 
232 }
233 
#define LogDebug(id)
RunNumber_t run() const
Definition: EventID.h:39
MonitorElement * mePhiHit_
edm::EDGetTokenT< edm::PCaloHitContainer > tok_hits_
EventNumber_t event() const
Definition: EventID.h:41
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
MonitorElement * meHEDepHit_
MonitorElement * meHOL10Ene_
int ib
Definition: cuy.py:660
MonitorElement * meHOPhiHit_
MonitorElement * bookProfile(Args &&...args)
Definition: DQMStore.h:157
MonitorElement * meDetectHit_
MonitorElement * meHFNHit_
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:462
MonitorElement * meHFEtaHit_
MonitorElement * meHFL10EneP_
MonitorElement * meHFEneHit2_
MonitorElement * meHOL10EneP_
MonitorElement * meHBNHit_
MonitorElement * meHFTimHit_
MonitorElement * meHBEneHit2_
MonitorElement * meHONHit_
void Fill(long long x)
MonitorElement * meAllNHit_
MonitorElement * meHBL10Ene_
MonitorElement * meBadDetHit_
MonitorElement * meHFL10Ene_
MonitorElement * mePhiHitb_
MonitorElement * meHETimHit_
MonitorElement * meTimeHit_
void analyzeHits(std::vector< PCaloHit > &)
MonitorElement * meHFDepHit_
MonitorElement * meHBTimHit_
MonitorElement * book1D(Args &&...args)
Definition: DQMStore.h:115
MonitorElement * meTimeWHit_
bool isValid() const
Definition: HandleBase.h:75
MonitorElement * meHBDepHit_
void analyze(const edm::Event &e, const edm::EventSetup &c)
MonitorElement * meBadSubHit_
virtual void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &)
MonitorElement * meHENHit_
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:276
MonitorElement * meBadIdHit_
MonitorElement * meHEL10EneP_
MonitorElement * meEnergyHit_
std::string hcalHits
MonitorElement * meHOEtaHit_
MonitorElement * meEtaHit_
MonitorElement * meHBEneHit_
MonitorElement * meHOTimHit_
MonitorElement * meHBL10EneP_
edm::EventID id() const
Definition: EventBase.h:59
MonitorElement * meHEEneHit2_
MonitorElement * meHFPhiHit_
std::string outFile_
MonitorElement * meDepthHit_
MonitorElement * meHBEtaHit_
MonitorElement * meHODepHit_
MonitorElement * meHEL10Ene_
MonitorElement * meHEEneHit_
std::string g4Label
MonitorElement * meHEPhiHit_
HcalSimHitStudy(const edm::ParameterSet &ps)
MonitorElement * meSubdetHit_
MonitorElement * meHBPhiHit_
Definition: Run.h:43
MonitorElement * meHFEneHit_
MonitorElement * meHEEtaHit_
MonitorElement * meHOEneHit_
MonitorElement * meHOEneHit2_