CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
EcalDigisValidation.cc
Go to the documentation of this file.
1 /*
2  * \file EcalDigisValidation.cc
3  *
4  * $Date: 2010/01/04 15:10:59 $
5  * $Revision: 1.30 $
6  * \author F. Cossutti
7  *
8 */
9 
16 
17 using namespace cms;
18 using namespace edm;
19 using namespace std;
20 
22  HepMCLabel(ps.getParameter<std::string>("moduleLabelMC")),
23  g4InfoLabel(ps.getParameter<std::string>("moduleLabelG4")),
24  EBdigiCollection_(ps.getParameter<edm::InputTag>("EBdigiCollection")),
25  EEdigiCollection_(ps.getParameter<edm::InputTag>("EEdigiCollection")),
26  ESdigiCollection_(ps.getParameter<edm::InputTag>("ESdigiCollection")){
27 
28 
29  // DQM ROOT output
30  outputFile_ = ps.getUntrackedParameter<string>("outputFile", "");
31 
32  if ( outputFile_.size() != 0 ) {
33  LogInfo("OutputInfo") << " Ecal Digi Task histograms will be saved to '" << outputFile_.c_str() << "'";
34  } else {
35  LogInfo("OutputInfo") << " Ecal Digi Task histograms will NOT be saved";
36  }
37 
38  // verbosity switch
39  verbose_ = ps.getUntrackedParameter<bool>("verbose", false);
40 
41  dbe_ = 0;
42 
43  // get hold of back-end interface
44  dbe_ = Service<DQMStore>().operator->();
45 
46  if ( dbe_ ) {
47  if ( verbose_ ) {
48  dbe_->setVerbose(1);
49  } else {
50  dbe_->setVerbose(0);
51  }
52  }
53 
54  if ( dbe_ ) {
55  if ( verbose_ ) dbe_->showDirStructure();
56  }
57 
58  gainConv_[1] = 1.;
59  gainConv_[2] = 2.;
60  gainConv_[3] = 12.;
61  gainConv_[0] = 12.; // saturated channels
62  barrelADCtoGeV_ = 0.035;
63  endcapADCtoGeV_ = 0.06;
64 
65  meGunEnergy_ = 0;
66  meGunEta_ = 0;
67  meGunPhi_ = 0;
68 
71 
74 
77 
78  Char_t histo[200];
79 
80 
81  if ( dbe_ ) {
82  dbe_->setCurrentFolder("EcalDigisV/EcalDigiTask");
83 
84  sprintf (histo, "EcalDigiTask Gun Momentum" ) ;
85  meGunEnergy_ = dbe_->book1D(histo, histo, 100, 0., 1000.);
86 
87  sprintf (histo, "EcalDigiTask Gun Eta" ) ;
88  meGunEta_ = dbe_->book1D(histo, histo, 700, -3.5, 3.5);
89 
90  sprintf (histo, "EcalDigiTask Gun Phi" ) ;
91  meGunPhi_ = dbe_->book1D(histo, histo, 360, 0., 360.);
92 
93  sprintf (histo, "EcalDigiTask Barrel maximum Digi over Sim ratio" ) ;
94  meEBDigiSimRatio_ = dbe_->book1D(histo, histo, 100, 0., 2.) ;
95 
96  sprintf (histo, "EcalDigiTask Endcap maximum Digi over Sim ratio" ) ;
97  meEEDigiSimRatio_ = dbe_->book1D(histo, histo, 100, 0., 2.) ;
98 
99  sprintf (histo, "EcalDigiTask Barrel maximum Digi over Sim ratio gt 10 ADC" ) ;
100  meEBDigiSimRatiogt10ADC_ = dbe_->book1D(histo, histo, 100, 0., 2.) ;
101 
102  sprintf (histo, "EcalDigiTask Endcap maximum Digi over Sim ratio gt 20 ADC" ) ;
103  meEEDigiSimRatiogt20ADC_ = dbe_->book1D(histo, histo, 100, 0., 2.) ;
104 
105  sprintf (histo, "EcalDigiTask Barrel maximum Digi over Sim ratio gt 100 ADC" ) ;
106  meEBDigiSimRatiogt100ADC_ = dbe_->book1D(histo, histo, 100, 0., 2.) ;
107 
108  sprintf (histo, "EcalDigiTask Endcap maximum Digi over Sim ratio gt 100 ADC" ) ;
109  meEEDigiSimRatiogt100ADC_ = dbe_->book1D(histo, histo, 100, 0., 2.) ;
110 
111  }
112 
113 }
114 
116 
117  if ( outputFile_.size() != 0 && dbe_ ) dbe_->save(outputFile_);
118 
119 }
120 
122 
124 
125 }
126 
128 
129 }
130 
132 
133  LogInfo("EventInfo") << " Run = " << e.id().run() << " Event = " << e.id().event();
134 
135  vector<SimTrack> theSimTracks;
136  vector<SimVertex> theSimVertexes;
137 
138  Handle<HepMCProduct> MCEvt;
141  Handle<CrossingFrame<PCaloHit> > crossingFrame;
142  Handle<EBDigiCollection> EcalDigiEB;
143  Handle<EEDigiCollection> EcalDigiEE;
144  Handle<ESDigiCollection> EcalDigiES;
145 
146  bool skipMC = false;
147  e.getByLabel(HepMCLabel, MCEvt);
148  if (!MCEvt.isValid()) { skipMC = true; }
149  e.getByLabel(g4InfoLabel,SimTk);
150  e.getByLabel(g4InfoLabel,SimVtx);
151 
152  const EBDigiCollection* EBdigis =0;
153  const EEDigiCollection* EEdigis =0;
154  const ESDigiCollection* ESdigis =0;
155 
156  bool isBarrel = true;
157  e.getByLabel( EBdigiCollection_, EcalDigiEB );
158  if (EcalDigiEB.isValid()) {
159  EBdigis = EcalDigiEB.product();
160  LogDebug("DigiInfo") << "total # EBdigis: " << EBdigis->size() ;
161  if ( EBdigis->size() == 0 ) isBarrel = false;
162  } else {
163  isBarrel = false;
164  }
165 
166  bool isEndcap = true;
167  e.getByLabel( EEdigiCollection_, EcalDigiEE );
168  if (EcalDigiEE.isValid()) {
169  EEdigis = EcalDigiEE.product();
170  LogDebug("DigiInfo") << "total # EEdigis: " << EEdigis->size() ;
171  if ( EEdigis->size() == 0 ) isEndcap = false;
172  } else {
173  isEndcap = false;
174  }
175 
176  bool isPreshower = true;
177  e.getByLabel( ESdigiCollection_, EcalDigiES );
178  if (EcalDigiES.isValid()) {
179  ESdigis = EcalDigiES.product();
180  LogDebug("DigiInfo") << "total # ESdigis: " << ESdigis->size() ;
181  if ( ESdigis->size() == 0 ) isPreshower = false;
182  } else {
183  isPreshower = false;
184  }
185 
186  theSimTracks.insert(theSimTracks.end(),SimTk->begin(),SimTk->end());
187  theSimVertexes.insert(theSimVertexes.end(),SimVtx->begin(),SimVtx->end());
188 
189  if ( ! skipMC ) {
190  double theGunEnergy = 0.;
191  for ( HepMC::GenEvent::particle_const_iterator p = MCEvt->GetEvent()->particles_begin();
192  p != MCEvt->GetEvent()->particles_end(); ++p ) {
193 
194  theGunEnergy = (*p)->momentum().e();
195  double htheta = (*p)->momentum().theta();
196  double heta = -log(tan(htheta * 0.5));
197  double hphi = (*p)->momentum().phi();
198  hphi = (hphi>=0) ? hphi : hphi+2*M_PI;
199  hphi = hphi / M_PI * 180.;
200  LogDebug("EventInfo") << "Particle gun type form MC = " << abs((*p)->pdg_id()) << "\n" << "Energy = "<< (*p)->momentum().e() << " Eta = " << heta << " Phi = " << hphi;
201 
202  if (meGunEnergy_) meGunEnergy_->Fill(theGunEnergy);
203  if (meGunEta_) meGunEta_->Fill(heta);
204  if (meGunPhi_) meGunPhi_->Fill(hphi);
205 
206  }
207  }
208 
209  int nvtx = 0;
210  for (vector<SimVertex>::iterator isimvtx = theSimVertexes.begin();
211  isimvtx != theSimVertexes.end(); ++isimvtx){
212  LogDebug("EventInfo") <<" Vertex index = " << nvtx << " event Id = " << isimvtx->eventId().rawId() << "\n" << " vertex dump: " << *isimvtx ;
213  ++nvtx;
214  }
215 
216  int ntrk = 0;
217  for (vector<SimTrack>::iterator isimtrk = theSimTracks.begin();
218  isimtrk != theSimTracks.end(); ++isimtrk){
219  LogDebug("EventInfo") <<" Track index = " << ntrk << " track Id = " << isimtrk->trackId() << " event Id = " << isimtrk->eventId().rawId() << "\n" << " track dump: " << *isimtrk ;
220  ++ntrk;
221  }
222 
223  // BARREL
224 
225  // loop over simHits
226 
227  if ( isBarrel ) {
228 
229  const std::string barrelHitsName(g4InfoLabel+"EcalHitsEB");
230  e.getByLabel("mix",barrelHitsName,crossingFrame);
231  std::auto_ptr<MixCollection<PCaloHit> >
232  barrelHits (new MixCollection<PCaloHit>(crossingFrame.product ()));
233 
234  MapType ebSimMap;
235 
236  for (MixCollection<PCaloHit>::MixItr hitItr = barrelHits->begin () ;
237  hitItr != barrelHits->end () ;
238  ++hitItr) {
239 
240  EBDetId ebid = EBDetId(hitItr->id()) ;
241 
242  LogDebug("HitInfo")
243  << " CaloHit " << hitItr->getName() << "\n"
244  << " DetID = " << hitItr->id()<< " EBDetId = " << ebid.ieta() << " " << ebid.iphi() << "\n"
245  << " Time = " << hitItr->time() << " Event id. = " << hitItr->eventId().rawId() << "\n"
246  << " Track Id = " << hitItr->geantTrackId() << "\n"
247  << " Energy = " << hitItr->energy();
248 
249  uint32_t crystid = ebid.rawId();
250  ebSimMap[crystid] += hitItr->energy();
251 
252  }
253 
254  // loop over Digis
255 
256  const EBDigiCollection * barrelDigi = EcalDigiEB.product () ;
257 
258  std::vector<double> ebAnalogSignal ;
259  std::vector<double> ebADCCounts ;
260  std::vector<double> ebADCGains ;
261  ebAnalogSignal.reserve(EBDataFrame::MAXSAMPLES);
262  ebADCCounts.reserve(EBDataFrame::MAXSAMPLES);
263  ebADCGains.reserve(EBDataFrame::MAXSAMPLES);
264 
265  for (unsigned int digis=0; digis<EcalDigiEB->size(); ++digis) {
266 
267  EBDataFrame ebdf=(*barrelDigi)[digis];
268  int nrSamples=ebdf.size();
269 
270  EBDetId ebid = ebdf.id () ;
271 
272  double Emax = 0. ;
273  int Pmax = 0 ;
274  double pedestalPreSample = 0.;
275  double pedestalPreSampleAnalog = 0.;
276 
277  for (int sample = 0 ; sample < nrSamples; ++sample) {
278  ebAnalogSignal[sample] = 0.;
279  ebADCCounts[sample] = 0.;
280  ebADCGains[sample] = -1.;
281  }
282 
283  for (int sample = 0 ; sample < nrSamples; ++sample) {
284 
286 
287  ebADCCounts[sample] = (mySample.adc()) ;
288  ebADCGains[sample] = (mySample.gainId()) ;
289  ebAnalogSignal[sample] = (ebADCCounts[sample]*gainConv_[(int)ebADCGains[sample]]*barrelADCtoGeV_);
290  if (Emax < ebAnalogSignal[sample] ) {
291  Emax = ebAnalogSignal[sample] ;
292  Pmax = sample ;
293  }
294  if ( sample < 3 ) {
295  pedestalPreSample += ebADCCounts[sample] ;
296  pedestalPreSampleAnalog += ebADCCounts[sample]*gainConv_[(int)ebADCGains[sample]]*barrelADCtoGeV_ ;
297  }
298  LogDebug("DigiInfo") << "EB sample " << sample << " ADC counts = " << ebADCCounts[sample] << " Gain Id = " << ebADCGains[sample] << " Analog eq = " << ebAnalogSignal[sample];
299  }
300 
301  pedestalPreSample /= 3. ;
302  pedestalPreSampleAnalog /= 3. ;
303  double Erec = Emax - pedestalPreSampleAnalog*gainConv_[(int)ebADCGains[Pmax]];
304 
305  if ( ebSimMap[ebid.rawId()] != 0. ) {
306  LogDebug("DigiInfo") << " Digi / Hit = " << Erec << " / " << ebSimMap[ebid.rawId()] << " gainConv " << gainConv_[(int)ebADCGains[Pmax]];
307  if ( meEBDigiSimRatio_ ) meEBDigiSimRatio_->Fill( Erec/ebSimMap[ebid.rawId()] ) ;
308  if ( Erec > 10.*barrelADCtoGeV_ && meEBDigiSimRatiogt10ADC_ ) meEBDigiSimRatiogt10ADC_->Fill( Erec/ebSimMap[ebid.rawId()] );
309  if ( Erec > 100.*barrelADCtoGeV_ && meEBDigiSimRatiogt100ADC_ ) meEBDigiSimRatiogt100ADC_->Fill( Erec/ebSimMap[ebid.rawId()] );
310 
311  }
312 
313  }
314 
315  }
316 
317  // ENDCAP
318 
319  // loop over simHits
320 
321  if ( isEndcap ) {
322 
323  const std::string endcapHitsName(g4InfoLabel+"EcalHitsEE");
324  e.getByLabel("mix",endcapHitsName,crossingFrame);
325  std::auto_ptr<MixCollection<PCaloHit> >
326  endcapHits (new MixCollection<PCaloHit>(crossingFrame.product ()));
327 
328  MapType eeSimMap;
329 
330  for (MixCollection<PCaloHit>::MixItr hitItr = endcapHits->begin () ;
331  hitItr != endcapHits->end () ;
332  ++hitItr) {
333 
334  EEDetId eeid = EEDetId(hitItr->id()) ;
335 
336  LogDebug("HitInfo")
337  << " CaloHit " << hitItr->getName() << "\n"
338  << " DetID = "<<hitItr->id()<< " EEDetId side = " << eeid.zside() << " = " << eeid.ix() << " " << eeid.iy() << "\n"
339  << " Time = " << hitItr->time() << " Event id. = " << hitItr->eventId().rawId() << "\n"
340  << " Track Id = " << hitItr->geantTrackId() << "\n"
341  << " Energy = " << hitItr->energy();
342 
343  uint32_t crystid = eeid.rawId();
344  eeSimMap[crystid] += hitItr->energy();
345 
346  }
347 
348  // loop over Digis
349 
350  const EEDigiCollection * endcapDigi = EcalDigiEE.product () ;
351 
352  std::vector<double> eeAnalogSignal ;
353  std::vector<double> eeADCCounts ;
354  std::vector<double> eeADCGains ;
355  eeAnalogSignal.reserve(EEDataFrame::MAXSAMPLES);
356  eeADCCounts.reserve(EEDataFrame::MAXSAMPLES);
357  eeADCGains.reserve(EEDataFrame::MAXSAMPLES);
358 
359  for (unsigned int digis=0; digis<EcalDigiEE->size(); ++digis) {
360 
361  EEDataFrame eedf=(*endcapDigi)[digis];
362  int nrSamples=eedf.size();
363 
364  EEDetId eeid = eedf.id () ;
365 
366  double Emax = 0. ;
367  int Pmax = 0 ;
368  double pedestalPreSample = 0.;
369  double pedestalPreSampleAnalog = 0.;
370 
371  for (int sample = 0 ; sample < nrSamples; ++sample) {
372  eeAnalogSignal[sample] = 0.;
373  eeADCCounts[sample] = 0.;
374  eeADCGains[sample] = -1.;
375  }
376 
377  for (int sample = 0 ; sample < nrSamples; ++sample) {
378 
380 
381  eeADCCounts[sample] = (mySample.adc()) ;
382  eeADCGains[sample] = (mySample.gainId()) ;
383  eeAnalogSignal[sample] = (eeADCCounts[sample]*gainConv_[(int)eeADCGains[sample]]*endcapADCtoGeV_);
384  if (Emax < eeAnalogSignal[sample] ) {
385  Emax = eeAnalogSignal[sample] ;
386  Pmax = sample ;
387  }
388  if ( sample < 3 ) {
389  pedestalPreSample += eeADCCounts[sample] ;
390  pedestalPreSampleAnalog += eeADCCounts[sample]*gainConv_[(int)eeADCGains[sample]]*endcapADCtoGeV_ ;
391  }
392  LogDebug("DigiInfo") << "EE sample " << sample << " ADC counts = " << eeADCCounts[sample] << " Gain Id = " << eeADCGains[sample] << " Analog eq = " << eeAnalogSignal[sample];
393  }
394  pedestalPreSample /= 3. ;
395  pedestalPreSampleAnalog /= 3. ;
396  double Erec = Emax - pedestalPreSampleAnalog*gainConv_[(int)eeADCGains[Pmax]];
397 
398  if (eeSimMap[eeid.rawId()] != 0. ) {
399  LogDebug("DigiInfo") << " Digi / Hit = " << Erec << " / " << eeSimMap[eeid.rawId()] << " gainConv " << gainConv_[(int)eeADCGains[Pmax]];
400  if ( meEEDigiSimRatio_) meEEDigiSimRatio_->Fill( Erec/eeSimMap[eeid.rawId()] ) ;
401  if ( Erec > 20.*endcapADCtoGeV_ && meEEDigiSimRatiogt20ADC_ ) meEEDigiSimRatiogt20ADC_->Fill( Erec/eeSimMap[eeid.rawId()] );
402  if ( Erec > 100.*endcapADCtoGeV_ && meEEDigiSimRatiogt100ADC_ ) meEEDigiSimRatiogt100ADC_->Fill( Erec/eeSimMap[eeid.rawId()] );
403  }
404 
405  }
406 
407  }
408 
409  if ( isPreshower) {
410 
411  const std::string preshowerHitsName(g4InfoLabel+"EcalHitsES");
412  e.getByLabel("mix",preshowerHitsName,crossingFrame);
413  std::auto_ptr<MixCollection<PCaloHit> >
414  preshowerHits (new MixCollection<PCaloHit>(crossingFrame.product ()));
415 
416  for (MixCollection<PCaloHit>::MixItr hitItr = preshowerHits->begin () ;
417  hitItr != preshowerHits->end () ;
418  ++hitItr) {
419 
420  ESDetId esid = ESDetId(hitItr->id()) ;
421 
422  LogDebug("HitInfo")
423  << " CaloHit " << hitItr->getName() << "\n"
424  << " DetID = " << hitItr->id()<< "ESDetId: z side " << esid.zside() << " plane " << esid.plane() << esid.six() << ',' << esid.siy() << ':' << esid.strip() << "\n"
425  << " Time = " << hitItr->time() << " Event id. = " << hitItr->eventId().rawId() << "\n"
426  << " Track Id = " << hitItr->geantTrackId() << "\n"
427  << " Energy = " << hitItr->energy();
428 
429  }
430 
431  }
432 
433 }
434 
436 {
437 
438  // ADC -> GeV Scale
440  eventSetup.get<EcalADCToGeVConstantRcd>().get(pAgc);
441  const EcalADCToGeVConstant* agc = pAgc.product();
442 
443  EcalMGPAGainRatio * defaultRatios = new EcalMGPAGainRatio();
444 
445  gainConv_[1] = 1.;
446  gainConv_[2] = defaultRatios->gain12Over6() ;
447  gainConv_[3] = gainConv_[2]*(defaultRatios->gain6Over1()) ;
448  gainConv_[0] = gainConv_[2]*(defaultRatios->gain6Over1()) ; // saturated channels
449 
450  LogDebug("EcalDigi") << " Gains conversions: " << "\n" << " g1 = " << gainConv_[1] << "\n" << " g2 = " << gainConv_[2] << "\n" << " g3 = " << gainConv_[3];
451  LogDebug("EcalDigi") << " Gains conversions: " << "\n" << " saturation = " << gainConv_[0];
452 
453  delete defaultRatios;
454 
455  const double barrelADCtoGeV_ = agc->getEBValue();
456  LogDebug("EcalDigi") << " Barrel GeV/ADC = " << barrelADCtoGeV_;
457  const double endcapADCtoGeV_ = agc->getEEValue();
458  LogDebug("EcalDigi") << " Endcap GeV/ADC = " << endcapADCtoGeV_;
459 
460 }
#define LogDebug(id)
RunNumber_t run() const
Definition: EventID.h:42
edm::InputTag EEdigiCollection_
EventNumber_t event() const
Definition: EventID.h:44
T getUntrackedParameter(std::string const &, T const &) const
int strip() const
Definition: ESDetId.h:41
int ix() const
Definition: EEDetId.h:71
key_type id() const
Definition: EBDataFrame.h:32
std::map< uint32_t, float, std::less< uint32_t > > MapType
void save(const std::string &filename, const std::string &path="", const std::string &pattern="", const std::string &rewrite="", SaveReferenceTag ref=SaveWithReference, int minStatus=dqm::qstatus::STATUS_OK, const std::string &fileupdate="RECREATE")
Definition: DQMStore.cc:2113
~EcalDigisValidation()
Destructor.
#define abs(x)
Definition: mlp_lapack.h:159
int six() const
Definition: ESDetId.h:37
MonitorElement * meEEDigiSimRatiogt100ADC_
std::map< int, double, std::less< int > > gainConv_
int gainId() const
get the gainId (2 bits)
int size() const
Definition: EcalDataFrame.h:27
MonitorElement * meGunEnergy_
void Fill(long long x)
int iphi() const
get the crystal iphi
Definition: EBDetId.h:46
uint32_t rawId() const
get the raw id
Definition: DetId.h:45
void beginRun(edm::Run const &, edm::EventSetup const &c)
MonitorElement * meEBDigiSimRatiogt10ADC_
MonitorElement * meEBDigiSimRatio_
int siy() const
Definition: ESDetId.h:39
void analyze(edm::Event const &e, edm::EventSetup const &c)
Analyze.
int zside() const
Definition: EEDetId.h:65
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
int iy() const
Definition: EEDetId.h:77
float gain6Over1() const
int ieta() const
get the crystal ieta
Definition: EBDetId.h:44
int zside() const
Definition: ESDetId.h:33
bool isValid() const
Definition: HandleBase.h:76
MonitorElement * meGunEta_
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:356
key_type id() const
Definition: EEDataFrame.h:29
#define M_PI
Definition: BFit3D.cc:3
const T & get() const
Definition: EventSetup.h:55
T const * product() const
Definition: ESHandle.h:62
float gain12Over6() const
EcalDigisValidation(const edm::ParameterSet &ps)
Constructor.
T const * product() const
Definition: Handle.h:74
MonitorElement * meEEDigiSimRatio_
edm::InputTag EBdigiCollection_
edm::InputTag ESdigiCollection_
edm::EventID id() const
Definition: EventBase.h:56
MonitorElement * meGunPhi_
int plane() const
Definition: ESDetId.h:35
MonitorElement * meEEDigiSimRatiogt20ADC_
void checkCalibrations(edm::EventSetup const &c)
MonitorElement * meEBDigiSimRatiogt100ADC_
static const int MAXSAMPLES
Definition: EcalDataFrame.h:49
Definition: Run.h:33
int adc() const
get the ADC sample (12 bits)