3232#include " Geometry/Records/interface/CaloGeometryRecord.h"
3333
3434#include < TH1F.h>
35+ #include < TH2F.h>
3536#include < cmath>
3637#include < iostream>
3738#include < fstream>
4041#include < string>
4142#include < vector>
4243
43- // #define EDM_ML_DEBUG
44+ #define EDM_ML_DEBUG
4445
4546class EcalSimHitStudy : public edm ::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
4647public:
@@ -74,6 +75,7 @@ class EcalSimHitStudy : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::on
7475 TH1F *etot_[ndets_], *etotg_[ndets_], *edepAll_[ndets_];
7576 TH1F *r1by9_[ndets_], *r1by25_[ndets_], *r9by25_[ndets_];
7677 TH1F *sEtaEta_ [ndets_], *sEtaPhi_ [ndets_], *sPhiPhi_ [ndets_];
78+ TH2F *poszp_[ndets_], *poszn_[ndets_];
7779};
7880
7981EcalSimHitStudy::EcalSimHitStudy (const edm::ParameterSet& ps) {
@@ -134,96 +136,119 @@ void EcalSimHitStudy::beginJob() {
134136 for (int i = 0 ; i < ndets_; i++) {
135137 sprintf (name, " Hit%d" , i);
136138 sprintf (title, " Number of hits in %s" , dets[i].c_str ());
137- hit_[i] = tfile->make <TH1F>(name, title , 1000 , 0 ., 20000 .);
139+ hit_[i] = tfile->make <TH1F>(name, dets[i]. c_str () , 1000 , 0 ., 20000 .);
138140 hit_[i]->GetXaxis ()->SetTitle (title);
139141 hit_[i]->GetYaxis ()->SetTitle (" Events" );
140142 hit_[i]->Sumw2 ();
141143 sprintf (name, " Time%d" , i);
142144 sprintf (title, " Time of the hit (ns) in %s" , dets[i].c_str ());
143- time_[i] = tfile->make <TH1F>(name, title , 1000 , 0 ., 1000 .);
145+ time_[i] = tfile->make <TH1F>(name, dets[i]. c_str () , 1000 , 0 ., 1000 .);
144146 time_[i]->GetXaxis ()->SetTitle (title);
145147 time_[i]->GetYaxis ()->SetTitle (" Hits" );
146148 time_[i]->Sumw2 ();
147149 sprintf (name, " TimeAll%d" , i);
148150 sprintf (title, " Hit time (ns) in %s (for first hit in the cell)" , dets[i].c_str ());
149- timeAll_[i] = tfile->make <TH1F>(name, title , 1000 , 0 ., 1000 .);
151+ timeAll_[i] = tfile->make <TH1F>(name, dets[i]. c_str () , 1000 , 0 ., 1000 .);
150152 timeAll_[i]->GetXaxis ()->SetTitle (title);
151153 timeAll_[i]->GetYaxis ()->SetTitle (" Hits" );
152154 timeAll_[i]->Sumw2 ();
153155 sprintf (name, " Edep%d" , i);
154156 sprintf (title, " Energy deposit (GeV) in %s" , dets[i].c_str ());
155- edep_[i] = tfile->make <TH1F>(name, title , 5000 , 0 ., maxEnergy_);
157+ edep_[i] = tfile->make <TH1F>(name, dets[i]. c_str () , 5000 , 0 ., maxEnergy_);
156158 edep_[i]->GetXaxis ()->SetTitle (title);
157159 edep_[i]->GetYaxis ()->SetTitle (" Hits" );
158160 edep_[i]->Sumw2 ();
159161 sprintf (name, " EdepAll%d" , i);
160162 sprintf (title, " Total Energy deposit in the cell (GeV) in %s" , dets[i].c_str ());
161- edepAll_[i] = tfile->make <TH1F>(name, title , 5000 , 0 ., maxEnergy_);
163+ edepAll_[i] = tfile->make <TH1F>(name, dets[i]. c_str () , 5000 , 0 ., maxEnergy_);
162164 edepAll_[i]->GetXaxis ()->SetTitle (title);
163165 edepAll_[i]->GetYaxis ()->SetTitle (" Hits" );
164166 edepAll_[i]->Sumw2 ();
165167 sprintf (name, " EdepEM%d" , i);
166168 sprintf (title, " Energy deposit (GeV) by EM particles in %s" , dets[i].c_str ());
167- edepEM_[i] = tfile->make <TH1F>(name, title , 5000 , 0 ., maxEnergy_);
169+ edepEM_[i] = tfile->make <TH1F>(name, dets[i]. c_str () , 5000 , 0 ., maxEnergy_);
168170 edepEM_[i]->GetXaxis ()->SetTitle (title);
169171 edepEM_[i]->GetYaxis ()->SetTitle (" Hits" );
170172 edepEM_[i]->Sumw2 ();
171173 sprintf (name, " EdepHad%d" , i);
172174 sprintf (title, " Energy deposit (GeV) by hadrons in %s" , dets[i].c_str ());
173- edepHad_[i] = tfile->make <TH1F>(name, title , 5000 , 0 ., maxEnergy_);
175+ edepHad_[i] = tfile->make <TH1F>(name, dets[i]. c_str () , 5000 , 0 ., maxEnergy_);
174176 edepHad_[i]->GetXaxis ()->SetTitle (title);
175177 edepHad_[i]->GetYaxis ()->SetTitle (" Hits" );
176178 edepHad_[i]->Sumw2 ();
177179 sprintf (name, " Etot%d" , i);
178180 sprintf (title, " Total energy deposit (GeV) in %s" , dets[i].c_str ());
179- etot_[i] = tfile->make <TH1F>(name, title , 5000 , 0 ., maxEnergy_);
181+ etot_[i] = tfile->make <TH1F>(name, dets[i]. c_str () , 5000 , 0 ., maxEnergy_);
180182 etot_[i]->GetXaxis ()->SetTitle (title);
181183 etot_[i]->GetYaxis ()->SetTitle (" Events" );
182184 etot_[i]->Sumw2 ();
183185 sprintf (name, " EtotG%d" , i);
184186 sprintf (title, " Total energy deposit (GeV) in %s (t < 100 ns)" , dets[i].c_str ());
185- etotg_[i] = tfile->make <TH1F>(name, title , 5000 , 0 ., maxEnergy_);
187+ etotg_[i] = tfile->make <TH1F>(name, dets[i]. c_str () , 5000 , 0 ., maxEnergy_);
186188 etotg_[i]->GetXaxis ()->SetTitle (title);
187189 etotg_[i]->GetYaxis ()->SetTitle (" Events" );
188190 etotg_[i]->Sumw2 ();
189191 sprintf (name, " r1by9%d" , i);
190192 sprintf (title, " E1/E9 in %s" , dets[i].c_str ());
191- r1by9_[i] = tfile->make <TH1F>(name, title , 100 , 0.0 , 1.0 );
193+ r1by9_[i] = tfile->make <TH1F>(name, dets[i]. c_str () , 100 , 0.0 , 1.0 );
192194 r1by9_[i]->GetXaxis ()->SetTitle (title);
193195 r1by9_[i]->GetYaxis ()->SetTitle (" Events" );
194196 r1by9_[i]->Sumw2 ();
195197 sprintf (name, " r1by25%d" , i);
196198 sprintf (title, " E1/E25 in %s" , dets[i].c_str ());
197- r1by25_[i] = tfile->make <TH1F>(name, title , 100 , 0.0 , 1.0 );
199+ r1by25_[i] = tfile->make <TH1F>(name, dets[i]. c_str () , 100 , 0.0 , 1.0 );
198200 r1by25_[i]->GetXaxis ()->SetTitle (title);
199201 r1by25_[i]->GetYaxis ()->SetTitle (" Events" );
200202 r1by25_[i]->Sumw2 ();
201203 sprintf (name, " r9by25%d" , i);
202204 sprintf (title, " E9/E25 in %s" , dets[i].c_str ());
203- r9by25_[i] = tfile->make <TH1F>(name, title , 100 , 0.0 , 1.0 );
205+ r9by25_[i] = tfile->make <TH1F>(name, dets[i]. c_str () , 100 , 0.0 , 1.0 );
204206 r9by25_[i]->GetXaxis ()->SetTitle (title);
205207 r9by25_[i]->GetYaxis ()->SetTitle (" Events" );
206208 r9by25_[i]->Sumw2 ();
207209 double ymax = (i == 0 ) ? 0.0005 : 0.005 ;
208210 sprintf (name, " sEtaEta%d" , i);
209211 sprintf (title, " Cov(#eta,#eta) in %s" , dets[i].c_str ());
210- sEtaEta_ [i] = tfile->make <TH1F>(name, title , 1000 , 0.0 , ymax);
212+ sEtaEta_ [i] = tfile->make <TH1F>(name, dets[i]. c_str () , 1000 , 0.0 , ymax);
211213 sEtaEta_ [i]->GetXaxis ()->SetTitle (title);
212214 sEtaEta_ [i]->GetYaxis ()->SetTitle (" Events" );
213215 sEtaEta_ [i]->Sumw2 ();
214216 sprintf (name, " sEtaPhi%d" , i);
215217 sprintf (title, " Cov(#eta,#phi) in %s" , dets[i].c_str ());
216- sEtaPhi_ [i] = tfile->make <TH1F>(name, title , 1000 , 0.0 , ymax);
218+ sEtaPhi_ [i] = tfile->make <TH1F>(name, dets[i]. c_str () , 1000 , 0.0 , ymax);
217219 sEtaPhi_ [i]->GetXaxis ()->SetTitle (title);
218220 sEtaPhi_ [i]->GetYaxis ()->SetTitle (" Events" );
219221 sEtaPhi_ [i]->Sumw2 ();
220222 ymax = (i == 0 ) ? 0.001 : 0.01 ;
221223 sprintf (name, " sPhiPhi%d" , i);
222224 sprintf (title, " Cov(#phi,#phi) in %s" , dets[i].c_str ());
223- sPhiPhi_ [i] = tfile->make <TH1F>(name, title , 1000 , 0.0 , ymax);
225+ sPhiPhi_ [i] = tfile->make <TH1F>(name, dets[i]. c_str () , 1000 , 0.0 , ymax);
224226 sPhiPhi_ [i]->GetXaxis ()->SetTitle (title);
225227 sPhiPhi_ [i]->GetYaxis ()->SetTitle (" Events" );
226228 sPhiPhi_ [i]->Sumw2 ();
229+ if (i == 0 ) {
230+ sprintf (title, " %s+" , dets[i].c_str ());
231+ poszp_[i] = tfile->make <TH2F>(" poszp0" , title, 100 , 0 , 100 , 360 , 0 , 360 );
232+ poszp_[i]->GetXaxis ()->SetTitle (" i#eta" );
233+ poszp_[i]->GetYaxis ()->SetTitle (" i#phi" );
234+ sprintf (title, " %s-" , dets[i].c_str ());
235+ poszn_[i] = tfile->make <TH2F>(" poszn0" , title, 100 , 0 , 100 , 360 , 0 , 360 );
236+ poszn_[i]->GetXaxis ()->SetTitle (" i#eta" );
237+ poszn_[i]->GetYaxis ()->SetTitle (" i#phi" );
238+ } else {
239+ sprintf (title, " %s+" , dets[i].c_str ());
240+ poszp_[i] = tfile->make <TH2F>(" poszp1" , title, 100 , -200 , 200 , 100 , -200 , 200 );
241+ poszp_[i]->GetXaxis ()->SetTitle (" x (cm)" );
242+ poszp_[i]->GetYaxis ()->SetTitle (" y (cm)" );
243+ sprintf (title, " %s-" , dets[i].c_str ());
244+ poszn_[i] = tfile->make <TH2F>(" poszn1" , title, 100 , -200 , 200 , 100 , -200 , 200 );
245+ poszn_[i]->GetXaxis ()->SetTitle (" x (cm)" );
246+ poszn_[i]->GetYaxis ()->SetTitle (" y (cm)" );
247+ }
248+ poszp_[i]->GetYaxis ()->SetTitleOffset (1.2 );
249+ poszp_[i]->Sumw2 ();
250+ poszn_[i]->GetYaxis ()->SetTitleOffset (1.2 );
251+ poszn_[i]->Sumw2 ();
227252 }
228253}
229254
@@ -236,34 +261,37 @@ void EcalSimHitStudy::analyze(const edm::Event& e, const edm::EventSetup& iS) {
236261 iS.get <CaloGeometryRecord>().get (pG);
237262 geometry_ = pG.product ();
238263
264+ double eInc = 0 , etaInc = 0 , phiInc = 0 ;
265+ int type (-1 );
239266 edm::Handle<edm::HepMCProduct> EvtHandle;
240267 e.getByToken (tok_evt_, EvtHandle);
241- const HepMC::GenEvent* myGenEvent = EvtHandle->GetEvent ();
268+ if (EvtHandle.isValid ()) {
269+ const HepMC::GenEvent* myGenEvent = EvtHandle->GetEvent ();
242270
243- double eInc = 0 , etaInc = 0 , phiInc = 0 ;
244- HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin ();
245- if (p != myGenEvent->particles_end ()) {
246- eInc = (*p)->momentum ().e ();
247- etaInc = (*p)->momentum ().eta ();
248- phiInc = (*p)->momentum ().phi ();
271+ HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin ();
272+ if (p != myGenEvent->particles_end ()) {
273+ eInc = (*p)->momentum ().e ();
274+ etaInc = (*p)->momentum ().eta ();
275+ phiInc = (*p)->momentum ().phi ();
276+ }
277+ double ptInc = eInc / std::cosh (etaInc);
278+ ptInc_->Fill (ptInc);
279+ eneInc_->Fill (eInc);
280+ etaInc_->Fill (etaInc);
281+ phiInc_->Fill (phiInc);
282+
283+ if (std::abs (etaInc) < 1.46 )
284+ type = 0 ;
285+ else if (std::abs (etaInc) > 1.49 && std::abs (etaInc) < 3.0 )
286+ type = 1 ;
249287 }
250- double ptInc = eInc / std::cosh (etaInc);
251- ptInc_->Fill (ptInc);
252- eneInc_->Fill (eInc);
253- etaInc_->Fill (etaInc);
254- phiInc_->Fill (phiInc);
255288
256- int type (-1 );
257- if (std::abs (etaInc) < 1.46 )
258- type = 0 ;
259- else if (std::abs (etaInc) > 1.49 && std::abs (etaInc) < 3.0 )
260- type = 1 ;
261- if (type >= 0 ) {
262- bool getHits (false );
289+ int typeMin = (type < 0 ) ? 0 : type;
290+ int typeMax = (type < 0 ) ? 1 : type;
291+ for (int type = typeMin; type <= typeMax; ++type) {
263292 edm::Handle<edm::PCaloHitContainer> hitsCalo;
264293 e.getByToken (toks_calo_[type], hitsCalo);
265- if (hitsCalo.isValid ())
266- getHits = true ;
294+ bool getHits = (hitsCalo.isValid ());
267295#ifdef EDM_ML_DEBUG
268296 edm::LogVerbatim (" HitStudy" ) << " EcalSimHitStudy: Input flags Hits " << getHits << " with " << hitsCalo->size ()
269297 << " hits" ;
@@ -340,6 +368,19 @@ void EcalSimHitStudy::analyzeHits(std::vector<PCaloHit>& hits, int indx) {
340368 for (auto it : hitMap) {
341369 timeAll_[indx]->Fill ((it.second ).time );
342370 edepAll_[indx]->Fill ((it.second ).energy );
371+ DetId id (it.first );
372+ if (indx == 0 ) {
373+ if (EBDetId (id).zside () >= 0 )
374+ poszp_[indx]->Fill (EBDetId (id).ietaAbs (), EBDetId (id).iphi ());
375+ else
376+ poszn_[indx]->Fill (EBDetId (id).ietaAbs (), EBDetId (id).iphi ());
377+ } else {
378+ GlobalPoint gpos = geometry_->getGeometry (id)->getPosition ();
379+ if (EEDetId (id).zside () >= 0 )
380+ poszp_[indx]->Fill (gpos.x (), gpos.y ());
381+ else
382+ poszn_[indx]->Fill (gpos.x (), gpos.y ());
383+ }
343384 }
344385
345386 math::XYZVector meanPosition (0.0 , 0.0 , 0.0 );
0 commit comments