00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include <Validation/EcalDigis/interface/EcalDigisValidation.h>
00011 #include <DataFormats/EcalDetId/interface/EBDetId.h>
00012 #include <DataFormats/EcalDetId/interface/EEDetId.h>
00013 #include <DataFormats/EcalDetId/interface/ESDetId.h>
00014 #include "CalibCalorimetry/EcalTrivialCondModules/interface/EcalTrivialConditionRetriever.h"
00015 #include "DQMServices/Core/interface/DQMStore.h"
00016
00017 using namespace cms;
00018 using namespace edm;
00019 using namespace std;
00020
00021 EcalDigisValidation::EcalDigisValidation(const ParameterSet& ps):
00022 HepMCLabel(ps.getParameter<std::string>("moduleLabelMC")),
00023 g4InfoLabel(ps.getParameter<std::string>("moduleLabelG4")),
00024 EBdigiCollection_(ps.getParameter<edm::InputTag>("EBdigiCollection")),
00025 EEdigiCollection_(ps.getParameter<edm::InputTag>("EEdigiCollection")),
00026 ESdigiCollection_(ps.getParameter<edm::InputTag>("ESdigiCollection")){
00027
00028
00029
00030 outputFile_ = ps.getUntrackedParameter<string>("outputFile", "");
00031
00032 if ( outputFile_.size() != 0 ) {
00033 LogInfo("OutputInfo") << " Ecal Digi Task histograms will be saved to '" << outputFile_.c_str() << "'";
00034 } else {
00035 LogInfo("OutputInfo") << " Ecal Digi Task histograms will NOT be saved";
00036 }
00037
00038
00039 verbose_ = ps.getUntrackedParameter<bool>("verbose", false);
00040
00041 dbe_ = 0;
00042
00043
00044 dbe_ = Service<DQMStore>().operator->();
00045
00046 if ( dbe_ ) {
00047 if ( verbose_ ) {
00048 dbe_->setVerbose(1);
00049 } else {
00050 dbe_->setVerbose(0);
00051 }
00052 }
00053
00054 if ( dbe_ ) {
00055 if ( verbose_ ) dbe_->showDirStructure();
00056 }
00057
00058 gainConv_[1] = 1.;
00059 gainConv_[2] = 2.;
00060 gainConv_[3] = 12.;
00061 gainConv_[0] = 12.;
00062 barrelADCtoGeV_ = 0.035;
00063 endcapADCtoGeV_ = 0.06;
00064
00065 meGunEnergy_ = 0;
00066 meGunEta_ = 0;
00067 meGunPhi_ = 0;
00068
00069 meEBDigiSimRatio_ = 0;
00070 meEEDigiSimRatio_ = 0;
00071
00072 meEBDigiSimRatiogt10ADC_ = 0;
00073 meEEDigiSimRatiogt20ADC_ = 0;
00074
00075 meEBDigiSimRatiogt100ADC_ = 0;
00076 meEEDigiSimRatiogt100ADC_ = 0;
00077
00078 Char_t histo[200];
00079
00080
00081 if ( dbe_ ) {
00082 dbe_->setCurrentFolder("EcalDigisV/EcalDigiTask");
00083
00084 sprintf (histo, "EcalDigiTask Gun Momentum" ) ;
00085 meGunEnergy_ = dbe_->book1D(histo, histo, 100, 0., 1000.);
00086
00087 sprintf (histo, "EcalDigiTask Gun Eta" ) ;
00088 meGunEta_ = dbe_->book1D(histo, histo, 700, -3.5, 3.5);
00089
00090 sprintf (histo, "EcalDigiTask Gun Phi" ) ;
00091 meGunPhi_ = dbe_->book1D(histo, histo, 360, 0., 360.);
00092
00093 sprintf (histo, "EcalDigiTask Barrel maximum Digi over Sim ratio" ) ;
00094 meEBDigiSimRatio_ = dbe_->book1D(histo, histo, 100, 0., 2.) ;
00095
00096 sprintf (histo, "EcalDigiTask Endcap maximum Digi over Sim ratio" ) ;
00097 meEEDigiSimRatio_ = dbe_->book1D(histo, histo, 100, 0., 2.) ;
00098
00099 sprintf (histo, "EcalDigiTask Barrel maximum Digi over Sim ratio gt 10 ADC" ) ;
00100 meEBDigiSimRatiogt10ADC_ = dbe_->book1D(histo, histo, 100, 0., 2.) ;
00101
00102 sprintf (histo, "EcalDigiTask Endcap maximum Digi over Sim ratio gt 20 ADC" ) ;
00103 meEEDigiSimRatiogt20ADC_ = dbe_->book1D(histo, histo, 100, 0., 2.) ;
00104
00105 sprintf (histo, "EcalDigiTask Barrel maximum Digi over Sim ratio gt 100 ADC" ) ;
00106 meEBDigiSimRatiogt100ADC_ = dbe_->book1D(histo, histo, 100, 0., 2.) ;
00107
00108 sprintf (histo, "EcalDigiTask Endcap maximum Digi over Sim ratio gt 100 ADC" ) ;
00109 meEEDigiSimRatiogt100ADC_ = dbe_->book1D(histo, histo, 100, 0., 2.) ;
00110
00111 }
00112
00113 }
00114
00115 EcalDigisValidation::~EcalDigisValidation(){
00116
00117 if ( outputFile_.size() != 0 && dbe_ ) dbe_->save(outputFile_);
00118
00119 }
00120
00121 void EcalDigisValidation::beginRun(Run const &, EventSetup const & c){
00122
00123 checkCalibrations(c);
00124
00125 }
00126
00127 void EcalDigisValidation::endJob(){
00128
00129 }
00130
00131 void EcalDigisValidation::analyze(Event const & e, EventSetup const & c){
00132
00133 LogInfo("EventInfo") << " Run = " << e.id().run() << " Event = " << e.id().event();
00134
00135 vector<SimTrack> theSimTracks;
00136 vector<SimVertex> theSimVertexes;
00137
00138 Handle<HepMCProduct> MCEvt;
00139 Handle<SimTrackContainer> SimTk;
00140 Handle<SimVertexContainer> SimVtx;
00141 Handle<CrossingFrame<PCaloHit> > crossingFrame;
00142 Handle<EBDigiCollection> EcalDigiEB;
00143 Handle<EEDigiCollection> EcalDigiEE;
00144 Handle<ESDigiCollection> EcalDigiES;
00145
00146 bool skipMC = false;
00147 e.getByLabel(HepMCLabel, MCEvt);
00148 if (!MCEvt.isValid()) { skipMC = true; }
00149 e.getByLabel(g4InfoLabel,SimTk);
00150 e.getByLabel(g4InfoLabel,SimVtx);
00151
00152 const EBDigiCollection* EBdigis =0;
00153 const EEDigiCollection* EEdigis =0;
00154 const ESDigiCollection* ESdigis =0;
00155
00156 bool isBarrel = true;
00157 e.getByLabel( EBdigiCollection_, EcalDigiEB );
00158 if (EcalDigiEB.isValid()) {
00159 EBdigis = EcalDigiEB.product();
00160 LogDebug("DigiInfo") << "total # EBdigis: " << EBdigis->size() ;
00161 if ( EBdigis->size() == 0 ) isBarrel = false;
00162 } else {
00163 isBarrel = false;
00164 }
00165
00166 bool isEndcap = true;
00167 e.getByLabel( EEdigiCollection_, EcalDigiEE );
00168 if (EcalDigiEE.isValid()) {
00169 EEdigis = EcalDigiEE.product();
00170 LogDebug("DigiInfo") << "total # EEdigis: " << EEdigis->size() ;
00171 if ( EEdigis->size() == 0 ) isEndcap = false;
00172 } else {
00173 isEndcap = false;
00174 }
00175
00176 bool isPreshower = true;
00177 e.getByLabel( ESdigiCollection_, EcalDigiES );
00178 if (EcalDigiES.isValid()) {
00179 ESdigis = EcalDigiES.product();
00180 LogDebug("DigiInfo") << "total # ESdigis: " << ESdigis->size() ;
00181 if ( ESdigis->size() == 0 ) isPreshower = false;
00182 } else {
00183 isPreshower = false;
00184 }
00185
00186 theSimTracks.insert(theSimTracks.end(),SimTk->begin(),SimTk->end());
00187 theSimVertexes.insert(theSimVertexes.end(),SimVtx->begin(),SimVtx->end());
00188
00189 if ( ! skipMC ) {
00190 double theGunEnergy = 0.;
00191 for ( HepMC::GenEvent::particle_const_iterator p = MCEvt->GetEvent()->particles_begin();
00192 p != MCEvt->GetEvent()->particles_end(); ++p ) {
00193
00194 theGunEnergy = (*p)->momentum().e();
00195 double htheta = (*p)->momentum().theta();
00196 double heta = -log(tan(htheta * 0.5));
00197 double hphi = (*p)->momentum().phi();
00198 hphi = (hphi>=0) ? hphi : hphi+2*M_PI;
00199 hphi = hphi / M_PI * 180.;
00200 LogDebug("EventInfo") << "Particle gun type form MC = " << abs((*p)->pdg_id()) << "\n" << "Energy = "<< (*p)->momentum().e() << " Eta = " << heta << " Phi = " << hphi;
00201
00202 if (meGunEnergy_) meGunEnergy_->Fill(theGunEnergy);
00203 if (meGunEta_) meGunEta_->Fill(heta);
00204 if (meGunPhi_) meGunPhi_->Fill(hphi);
00205
00206 }
00207 }
00208
00209 int nvtx = 0;
00210 for (vector<SimVertex>::iterator isimvtx = theSimVertexes.begin();
00211 isimvtx != theSimVertexes.end(); ++isimvtx){
00212 LogDebug("EventInfo") <<" Vertex index = " << nvtx << " event Id = " << isimvtx->eventId().rawId() << "\n" << " vertex dump: " << *isimvtx ;
00213 ++nvtx;
00214 }
00215
00216 int ntrk = 0;
00217 for (vector<SimTrack>::iterator isimtrk = theSimTracks.begin();
00218 isimtrk != theSimTracks.end(); ++isimtrk){
00219 LogDebug("EventInfo") <<" Track index = " << ntrk << " track Id = " << isimtrk->trackId() << " event Id = " << isimtrk->eventId().rawId() << "\n" << " track dump: " << *isimtrk ;
00220 ++ntrk;
00221 }
00222
00223
00224
00225
00226
00227 if ( isBarrel ) {
00228
00229 const std::string barrelHitsName(g4InfoLabel+"EcalHitsEB");
00230 e.getByLabel("mix",barrelHitsName,crossingFrame);
00231 std::auto_ptr<MixCollection<PCaloHit> >
00232 barrelHits (new MixCollection<PCaloHit>(crossingFrame.product ()));
00233
00234 MapType ebSimMap;
00235
00236 for (MixCollection<PCaloHit>::MixItr hitItr = barrelHits->begin () ;
00237 hitItr != barrelHits->end () ;
00238 ++hitItr) {
00239
00240 EBDetId ebid = EBDetId(hitItr->id()) ;
00241
00242 LogDebug("HitInfo")
00243 << " CaloHit " << hitItr->getName() << "\n"
00244 << " DetID = " << hitItr->id()<< " EBDetId = " << ebid.ieta() << " " << ebid.iphi() << "\n"
00245 << " Time = " << hitItr->time() << " Event id. = " << hitItr->eventId().rawId() << "\n"
00246 << " Track Id = " << hitItr->geantTrackId() << "\n"
00247 << " Energy = " << hitItr->energy();
00248
00249 uint32_t crystid = ebid.rawId();
00250 ebSimMap[crystid] += hitItr->energy();
00251
00252 }
00253
00254
00255
00256 const EBDigiCollection * barrelDigi = EcalDigiEB.product () ;
00257
00258 std::vector<double> ebAnalogSignal ;
00259 std::vector<double> ebADCCounts ;
00260 std::vector<double> ebADCGains ;
00261 ebAnalogSignal.reserve(EBDataFrame::MAXSAMPLES);
00262 ebADCCounts.reserve(EBDataFrame::MAXSAMPLES);
00263 ebADCGains.reserve(EBDataFrame::MAXSAMPLES);
00264
00265 for (unsigned int digis=0; digis<EcalDigiEB->size(); ++digis) {
00266
00267 EBDataFrame ebdf=(*barrelDigi)[digis];
00268 int nrSamples=ebdf.size();
00269
00270 EBDetId ebid = ebdf.id () ;
00271
00272 double Emax = 0. ;
00273 int Pmax = 0 ;
00274 double pedestalPreSample = 0.;
00275 double pedestalPreSampleAnalog = 0.;
00276
00277 for (int sample = 0 ; sample < nrSamples; ++sample) {
00278 ebAnalogSignal[sample] = 0.;
00279 ebADCCounts[sample] = 0.;
00280 ebADCGains[sample] = -1.;
00281 }
00282
00283 for (int sample = 0 ; sample < nrSamples; ++sample) {
00284
00285 EcalMGPASample mySample = ebdf[sample];
00286
00287 ebADCCounts[sample] = (mySample.adc()) ;
00288 ebADCGains[sample] = (mySample.gainId()) ;
00289 ebAnalogSignal[sample] = (ebADCCounts[sample]*gainConv_[(int)ebADCGains[sample]]*barrelADCtoGeV_);
00290 if (Emax < ebAnalogSignal[sample] ) {
00291 Emax = ebAnalogSignal[sample] ;
00292 Pmax = sample ;
00293 }
00294 if ( sample < 3 ) {
00295 pedestalPreSample += ebADCCounts[sample] ;
00296 pedestalPreSampleAnalog += ebADCCounts[sample]*gainConv_[(int)ebADCGains[sample]]*barrelADCtoGeV_ ;
00297 }
00298 LogDebug("DigiInfo") << "EB sample " << sample << " ADC counts = " << ebADCCounts[sample] << " Gain Id = " << ebADCGains[sample] << " Analog eq = " << ebAnalogSignal[sample];
00299 }
00300
00301 pedestalPreSample /= 3. ;
00302 pedestalPreSampleAnalog /= 3. ;
00303 double Erec = Emax - pedestalPreSampleAnalog*gainConv_[(int)ebADCGains[Pmax]];
00304
00305 if ( ebSimMap[ebid.rawId()] != 0. ) {
00306 LogDebug("DigiInfo") << " Digi / Hit = " << Erec << " / " << ebSimMap[ebid.rawId()] << " gainConv " << gainConv_[(int)ebADCGains[Pmax]];
00307 if ( meEBDigiSimRatio_ ) meEBDigiSimRatio_->Fill( Erec/ebSimMap[ebid.rawId()] ) ;
00308 if ( Erec > 10.*barrelADCtoGeV_ && meEBDigiSimRatiogt10ADC_ ) meEBDigiSimRatiogt10ADC_->Fill( Erec/ebSimMap[ebid.rawId()] );
00309 if ( Erec > 100.*barrelADCtoGeV_ && meEBDigiSimRatiogt100ADC_ ) meEBDigiSimRatiogt100ADC_->Fill( Erec/ebSimMap[ebid.rawId()] );
00310
00311 }
00312
00313 }
00314
00315 }
00316
00317
00318
00319
00320
00321 if ( isEndcap ) {
00322
00323 const std::string endcapHitsName(g4InfoLabel+"EcalHitsEE");
00324 e.getByLabel("mix",endcapHitsName,crossingFrame);
00325 std::auto_ptr<MixCollection<PCaloHit> >
00326 endcapHits (new MixCollection<PCaloHit>(crossingFrame.product ()));
00327
00328 MapType eeSimMap;
00329
00330 for (MixCollection<PCaloHit>::MixItr hitItr = endcapHits->begin () ;
00331 hitItr != endcapHits->end () ;
00332 ++hitItr) {
00333
00334 EEDetId eeid = EEDetId(hitItr->id()) ;
00335
00336 LogDebug("HitInfo")
00337 << " CaloHit " << hitItr->getName() << "\n"
00338 << " DetID = "<<hitItr->id()<< " EEDetId side = " << eeid.zside() << " = " << eeid.ix() << " " << eeid.iy() << "\n"
00339 << " Time = " << hitItr->time() << " Event id. = " << hitItr->eventId().rawId() << "\n"
00340 << " Track Id = " << hitItr->geantTrackId() << "\n"
00341 << " Energy = " << hitItr->energy();
00342
00343 uint32_t crystid = eeid.rawId();
00344 eeSimMap[crystid] += hitItr->energy();
00345
00346 }
00347
00348
00349
00350 const EEDigiCollection * endcapDigi = EcalDigiEE.product () ;
00351
00352 std::vector<double> eeAnalogSignal ;
00353 std::vector<double> eeADCCounts ;
00354 std::vector<double> eeADCGains ;
00355 eeAnalogSignal.reserve(EEDataFrame::MAXSAMPLES);
00356 eeADCCounts.reserve(EEDataFrame::MAXSAMPLES);
00357 eeADCGains.reserve(EEDataFrame::MAXSAMPLES);
00358
00359 for (unsigned int digis=0; digis<EcalDigiEE->size(); ++digis) {
00360
00361 EEDataFrame eedf=(*endcapDigi)[digis];
00362 int nrSamples=eedf.size();
00363
00364 EEDetId eeid = eedf.id () ;
00365
00366 double Emax = 0. ;
00367 int Pmax = 0 ;
00368 double pedestalPreSample = 0.;
00369 double pedestalPreSampleAnalog = 0.;
00370
00371 for (int sample = 0 ; sample < nrSamples; ++sample) {
00372 eeAnalogSignal[sample] = 0.;
00373 eeADCCounts[sample] = 0.;
00374 eeADCGains[sample] = -1.;
00375 }
00376
00377 for (int sample = 0 ; sample < nrSamples; ++sample) {
00378
00379 EcalMGPASample mySample = eedf[sample];
00380
00381 eeADCCounts[sample] = (mySample.adc()) ;
00382 eeADCGains[sample] = (mySample.gainId()) ;
00383 eeAnalogSignal[sample] = (eeADCCounts[sample]*gainConv_[(int)eeADCGains[sample]]*endcapADCtoGeV_);
00384 if (Emax < eeAnalogSignal[sample] ) {
00385 Emax = eeAnalogSignal[sample] ;
00386 Pmax = sample ;
00387 }
00388 if ( sample < 3 ) {
00389 pedestalPreSample += eeADCCounts[sample] ;
00390 pedestalPreSampleAnalog += eeADCCounts[sample]*gainConv_[(int)eeADCGains[sample]]*endcapADCtoGeV_ ;
00391 }
00392 LogDebug("DigiInfo") << "EE sample " << sample << " ADC counts = " << eeADCCounts[sample] << " Gain Id = " << eeADCGains[sample] << " Analog eq = " << eeAnalogSignal[sample];
00393 }
00394 pedestalPreSample /= 3. ;
00395 pedestalPreSampleAnalog /= 3. ;
00396 double Erec = Emax - pedestalPreSampleAnalog*gainConv_[(int)eeADCGains[Pmax]];
00397
00398 if (eeSimMap[eeid.rawId()] != 0. ) {
00399 LogDebug("DigiInfo") << " Digi / Hit = " << Erec << " / " << eeSimMap[eeid.rawId()] << " gainConv " << gainConv_[(int)eeADCGains[Pmax]];
00400 if ( meEEDigiSimRatio_) meEEDigiSimRatio_->Fill( Erec/eeSimMap[eeid.rawId()] ) ;
00401 if ( Erec > 20.*endcapADCtoGeV_ && meEEDigiSimRatiogt20ADC_ ) meEEDigiSimRatiogt20ADC_->Fill( Erec/eeSimMap[eeid.rawId()] );
00402 if ( Erec > 100.*endcapADCtoGeV_ && meEEDigiSimRatiogt100ADC_ ) meEEDigiSimRatiogt100ADC_->Fill( Erec/eeSimMap[eeid.rawId()] );
00403 }
00404
00405 }
00406
00407 }
00408
00409 if ( isPreshower) {
00410
00411 const std::string preshowerHitsName(g4InfoLabel+"EcalHitsES");
00412 e.getByLabel("mix",preshowerHitsName,crossingFrame);
00413 std::auto_ptr<MixCollection<PCaloHit> >
00414 preshowerHits (new MixCollection<PCaloHit>(crossingFrame.product ()));
00415
00416 for (MixCollection<PCaloHit>::MixItr hitItr = preshowerHits->begin () ;
00417 hitItr != preshowerHits->end () ;
00418 ++hitItr) {
00419
00420 ESDetId esid = ESDetId(hitItr->id()) ;
00421
00422 LogDebug("HitInfo")
00423 << " CaloHit " << hitItr->getName() << "\n"
00424 << " DetID = " << hitItr->id()<< "ESDetId: z side " << esid.zside() << " plane " << esid.plane() << esid.six() << ',' << esid.siy() << ':' << esid.strip() << "\n"
00425 << " Time = " << hitItr->time() << " Event id. = " << hitItr->eventId().rawId() << "\n"
00426 << " Track Id = " << hitItr->geantTrackId() << "\n"
00427 << " Energy = " << hitItr->energy();
00428
00429 }
00430
00431 }
00432
00433 }
00434
00435 void EcalDigisValidation::checkCalibrations(edm::EventSetup const & eventSetup)
00436 {
00437
00438
00439 edm::ESHandle<EcalADCToGeVConstant> pAgc;
00440 eventSetup.get<EcalADCToGeVConstantRcd>().get(pAgc);
00441 const EcalADCToGeVConstant* agc = pAgc.product();
00442
00443 EcalMGPAGainRatio * defaultRatios = new EcalMGPAGainRatio();
00444
00445 gainConv_[1] = 1.;
00446 gainConv_[2] = defaultRatios->gain12Over6() ;
00447 gainConv_[3] = gainConv_[2]*(defaultRatios->gain6Over1()) ;
00448 gainConv_[0] = gainConv_[2]*(defaultRatios->gain6Over1()) ;
00449
00450 LogDebug("EcalDigi") << " Gains conversions: " << "\n" << " g1 = " << gainConv_[1] << "\n" << " g2 = " << gainConv_[2] << "\n" << " g3 = " << gainConv_[3];
00451 LogDebug("EcalDigi") << " Gains conversions: " << "\n" << " saturation = " << gainConv_[0];
00452
00453 delete defaultRatios;
00454
00455 const double barrelADCtoGeV_ = agc->getEBValue();
00456 LogDebug("EcalDigi") << " Barrel GeV/ADC = " << barrelADCtoGeV_;
00457 const double endcapADCtoGeV_ = agc->getEEValue();
00458 LogDebug("EcalDigi") << " Endcap GeV/ADC = " << endcapADCtoGeV_;
00459
00460 }
1.7.3