00001
00002
00003
00004
00005
00006
00007
00008
00009 #include "Validation/EcalDigis/interface/EcalSelectiveReadoutValidation.h"
00010 #include "Validation/EcalDigis/src/EcalSRPCompat.h"
00011
00012 #include "Validation/EcalDigis/src/ecalDccMap.h"
00013
00014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00015 #include "FWCore/Framework/interface/ESHandle.h"
00016 #include "FWCore/Framework/interface/Event.h"
00017 #include "FWCore/Framework/interface/EventSetup.h"
00018 #include "FWCore/ServiceRegistry/interface/Service.h"
00019
00020 #include "Geometry/EcalMapping/interface/EcalElectronicsMapping.h"
00021 #include "Geometry/EcalMapping/interface/EcalMappingRcd.h"
00022 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00023 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
00024 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
00025
00026 #include "DQMServices/Core/interface/DQMStore.h"
00027
00028 #include "CalibCalorimetry/EcalTPGTools/interface/EcalTPGScale.h"
00029
00030 #include <string.h>
00031 #include "DQMServices/Core/interface/MonitorElement.h"
00032
00033 #if (CMSSW_COMPAT_VERSION>=210)
00034 # include "Geometry/Records/interface/CaloGeometryRecord.h"
00035 typedef CaloGeometryRecord MyCaloGeometryRecord;
00036 #else
00037 # include "Geometry/Records/interface/IdealGeometryRecord.h"
00038 typedef IdealGeometryRecord MyCaloGeometryRecord;
00039 #endif
00040
00041 #define ML_DEBUG
00042
00043 using namespace cms;
00044 using namespace edm;
00045 using namespace std;
00046
00047 const double EcalSelectiveReadoutValidation::rad2deg = 45./atan(1.);
00048
00049 const int EcalSelectiveReadoutValidation::nDccRus_[nDccs_] ={
00050
00051
00052 34, 32, 33, 33, 32, 34, 33, 34, 33,
00053
00054
00055 68, 68, 68, 68, 68, 68, 68, 68, 68, 68, 68, 68, 68, 68, 68, 68, 68, 68,
00056
00057
00058 68, 68, 68, 68, 68, 68, 68, 68, 68, 68, 68, 68, 68, 68, 68, 68, 68, 68,
00059
00060
00061 32, 33, 33, 32, 34, 33, 34, 33, 34
00062 };
00063
00064
00065
00066 #ifdef DO_TIMING
00067 # include <sys/time.h>
00068 struct PgTiming{
00069 PgTiming(const std::string& mess): mess_(mess), stopped_(false){
00070 gettimeofday(&start_, 0);
00071 }
00072
00073 ~PgTiming(){
00074 if(!stopped_) stop();
00075 }
00076 void stop(){
00077 timeval t;
00078 gettimeofday(&t, 0);
00079 std::cout << "[PgTiming] " << mess_ << " "
00080 << ((t.tv_sec-start_.tv_sec)*1.e3
00081 + (t.tv_usec-start_.tv_usec)*1.e-3)
00082 << " ms.\n";
00083 stopped_ = true;
00084 }
00085 timeval start_;
00086 std::string mess_;
00087 bool stopped_;
00088 };
00089 #else //DO_TIMING is not defined
00090 struct PgTiming{
00091 PgTiming(const std::string& mess){}
00092 void stop(){}
00093 };
00094 #endif //DO_TIMING defined condition
00095
00096 EcalSelectiveReadoutValidation::EcalSelectiveReadoutValidation(const ParameterSet& ps):
00097 collNotFoundWarn_(ps.getUntrackedParameter<bool>("warnIfCollectionNotFound", true)),
00098 ebDigis_(ps.getParameter<edm::InputTag>("EbDigiCollection"), false,
00099 collNotFoundWarn_),
00100 eeDigis_(ps.getParameter<edm::InputTag>("EeDigiCollection"), false,
00101 collNotFoundWarn_),
00102 ebNoZsDigis_(ps.getParameter<edm::InputTag>("EbUnsuppressedDigiCollection"),
00103 false, false),
00104 eeNoZsDigis_(ps.getParameter<edm::InputTag>("EeUnsuppressedDigiCollection"),
00105 false, false),
00106 ebSrFlags_(ps.getParameter<edm::InputTag>("EbSrFlagCollection"), false,
00107 collNotFoundWarn_),
00108 eeSrFlags_(ps.getParameter<edm::InputTag>("EeSrFlagCollection"), false,
00109 collNotFoundWarn_),
00110 ebComputedSrFlags_(ps.getParameter<edm::InputTag>("EbSrFlagFromTTCollection"), false,
00111 false),
00112 eeComputedSrFlags_(ps.getParameter<edm::InputTag>("EeSrFlagFromTTCollection"), false,
00113 false),
00114 ebSimHits_(ps.getParameter<edm::InputTag>("EbSimHitCollection"), false,
00115 false),
00116 eeSimHits_(ps.getParameter<edm::InputTag>("EeSimHitCollection"), false,
00117 false),
00118 tps_(ps.getParameter<edm::InputTag>("TrigPrimCollection"), false,
00119 collNotFoundWarn_),
00120 ebRecHits_(ps.getParameter<edm::InputTag>("EbRecHitCollection"), false,
00121 false),
00122 eeRecHits_(ps.getParameter<edm::InputTag>("EeRecHitCollection"), false,
00123 false),
00124 fedRaw_(ps.getParameter<edm::InputTag>("FEDRawCollection"), false,
00125 false),
00126 tmax(0),
00127 tmin(numeric_limits<int64_t>::max()),
00128 l1aOfTmin(0),
00129 l1aOfTmax(0),
00130 triggerTowerMap_(0),
00131 localReco_(ps.getParameter<bool>("LocalReco")),
00132 weights_(ps.getParameter<vector<double> >("weights")),
00133 tpInGeV_(ps.getParameter<bool>("tpInGeV")),
00134 firstFIRSample_(ps.getParameter<int>("ecalDccZs1stSample")),
00135 useEventRate_(ps.getParameter<bool>("useEventRate")),
00136 logErrForDccs_(nDccs_, false),
00137 ievt_(0),
00138 allHists_(false),
00139 histDir_(ps.getParameter<string>("histDir")),
00140 withEeSimHit_(false),
00141 withEbSimHit_(false){
00142
00143 PgTiming t("EcalSelectiveReadoutValidation ctor");
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160 double ebZsThr = ps.getParameter<double>("ebZsThrADCCount");
00161 double eeZsThr = ps.getParameter<double>("eeZsThrADCCount");
00162
00163 ebZsThr_ = lround(ebZsThr*4);
00164 eeZsThr_ = lround(eeZsThr*4);
00165
00166
00167 srpAlgoErrorLogFileName_
00168 = ps.getUntrackedParameter<string>("srpAlgoErrorLogFile","");
00169 logSrpAlgoErrors_ = (srpAlgoErrorLogFileName_.size()!=0);
00170
00171
00172 srApplicationErrorLogFileName_
00173 = ps.getUntrackedParameter<string>("srApplicationErrorLogFile","");
00174 logSrApplicationErrors_ = (srApplicationErrorLogFileName_.size()!=0);
00175
00176
00177 configFirWeights(ps.getParameter<vector<double> >("dccWeights"));
00178
00179
00180 outputFile_ = ps.getUntrackedParameter<string>("outputFile", "");
00181
00182 if(outputFile_.size() != 0){
00183 LogInfo("OutputInfo") << " Ecal Digi Task histograms will be saved to '"
00184 << outputFile_.c_str() << "'";
00185 } else{
00186 LogInfo("OutputInfo") << " Ecal Digi Task histograms will NOT be saved";
00187 }
00188
00189
00190 verbose_ = ps.getUntrackedParameter<bool>("verbose", false);
00191
00192
00193 dbe_ = Service<DQMStore>().operator->();
00194
00195 if(verbose_){
00196 dbe_->setVerbose(1);
00197 } else{
00198 dbe_->setVerbose(0);
00199 }
00200
00201 if(verbose_) dbe_->showDirStructure();
00202
00203 dbe_->setCurrentFolder(histDir_);
00204
00205 vector<string>
00206 hists(ps.getUntrackedParameter<vector<string> >("histograms",
00207 vector<string>(1, "all")));
00208
00209 for(vector<string>::iterator it = hists.begin();
00210 it!=hists.end(); ++it) histList_.insert(*it);
00211 if(histList_.find("all") != histList_.end()) allHists_ = true;
00212
00213
00214 meEbFixedPayload_ = bookFloat("ebFixedVol");
00215 double ebFixed = getDccOverhead(EB)*nEbDccs;
00216 double eeFixed = getDccOverhead(EE)*nEeDccs;
00217 meEbFixedPayload_->Fill(ebFixed);
00218 meEeFixedPayload_ = bookFloat("eeFixedVol");
00219 meEbFixedPayload_->Fill(eeFixed);
00220 meFixedPayload_ = bookFloat("fixedVol");
00221 meFixedPayload_->Fill(ebFixed+eeFixed);
00222
00223 meL1aRate_ = bookFloat("l1aRate_");
00224
00225 meDccVol_ = bookProfile("hDccVol",
00226 "ECAL DCC event fragment size;Dcc id; "
00227 "<Event size> (kB)", nDccs_, .5, .5+nDccs_);
00228
00229 meDccLiVol_ = bookProfile("hDccLiVol",
00230 "LI channel payload per DCC;Dcc id; "
00231 "<Event size> (kB)", nDccs_, .5, .5+nDccs_);
00232
00233 meDccHiVol_ = bookProfile("hDccHiVol",
00234 "HI channel payload per DCC;Dcc id; "
00235 "<Event size> (kB)", nDccs_, .5, .5+nDccs_);
00236
00237 meDccVolFromData_ = bookProfile("hDccVolFromData",
00238 "ECAL DCC event fragment size;Dcc id; "
00239 "<Event size> (kB)", nDccs_, .5, .5+nDccs_);
00240
00241 meVolBLI_ = book1D("hVolBLI",
00242 "ECAL Barrel low interest crystal data payload;"
00243 "Event size (kB);Nevts",
00244 100, 0., 200.);
00245
00246 meVolELI_ = book1D("hVolELI",
00247 "Endcap low interest crystal data payload;"
00248 "Event size (kB);Nevts",
00249 100, 0., 200.);
00250
00251 meVolLI_ = book1D("hVolLI",
00252 "ECAL low interest crystal data payload;"
00253 "Event size (kB);Nevts",
00254 100, 0., 200.);
00255
00256 meVolBHI_ = book1D("hVolBHI",
00257 "Barrel high interest crystal data payload;"
00258 "Event size (kB);Nevts",
00259 100, 0., 200.);
00260
00261 meVolEHI_ = book1D("hVolEHI",
00262 "Endcap high interest crystal data payload;"
00263 "Event size (kB);Nevts",
00264 100, 0., 200.);
00265
00266 meVolHI_ = book1D("hVolHI",
00267 "ECAL high interest crystal data payload;"
00268 "Event size (kB);Nevts",
00269 100, 0., 200.);
00270
00271 meVolB_ = book1D("hVolB",
00272 "Barrel data volume;Event size (kB);Nevts",
00273 100, 0., 200.);
00274
00275 meVolE_ = book1D("hVolE",
00276 "Endcap data volume;Event size (kB);Nevts",
00277 100, 0., 200.);
00278
00279 meVol_ = book1D("hVol",
00280 "ECAL data volume;Event size (kB);Nevts",
00281 100, 0., 200.);
00282
00283 meChOcc_ = book2D("h2ChOcc",
00284 "ECAL crystal channel occupancy after zero suppression;"
00285 "iX -200 / iEta / iX + 100;"
00286 "iY / iPhi (starting from -10^{o}!);"
00287 "Event count",
00288 401, -200.5, 200.5,
00289 360, .5, 360.5);
00290
00291
00292 string tpUnit;
00293 if(tpInGeV_) tpUnit = string("GeV"); else tpUnit = string("TP hw unit");
00294 string title;
00295 title = string("Trigger primitive TT E_{T};E_{T} ")
00296 + tpUnit + string(";Event Count");
00297 meTp_ = book1D("hTp",
00298 title.c_str(),
00299 (tpInGeV_?100:40), 0., (tpInGeV_?10.:40.));
00300
00301 meTtf_ = book1D("hTtf",
00302 "Trigger primitive TT flag;Flag number;Event count",
00303 8, -.5, 7.5);
00304
00305 title = string("Trigger tower flag vs TP;E_{T}(TT) (")
00306 + tpUnit + string(");Flag number");
00307 meTtfVsTp_ = book2D("h2TtfVsTp",
00308 title.c_str(),
00309 100, 0., (tpInGeV_?10.:40.),
00310 8, -.5, 7.5);
00311
00312 meTtfVsEtSum_ = book2D("h2TtfVsEtSum",
00313 "Trigger tower flag vs #sumE_{T};"
00314 "E_{T}(TT) (GeV);"
00315 "TTF",
00316 100, 0., 10.,
00317 8, -.5, 7.5);
00318 title = string("Trigger primitive Et (TP) vs #sumE_{T};"
00319 "E_{T} (sum) (GeV);"
00320 "E_{T} (TP) (") + tpUnit + string (")");
00321
00322 meTpVsEtSum_ = book2D("h2TpVsEtSum",
00323 title.c_str(),
00324 100, 0., 10.,
00325 100, 0., (tpInGeV_?10.:40.));
00326
00327 title = string("Trigger primitive E_{T};"
00328 "iEta;"
00329 "iPhi;"
00330 "E_{T} (TP) (") + tpUnit + string (")");
00331 meTpMap_ = bookProfile2D("h2Tp",
00332 title.c_str(),
00333 57, -28.5, 28.5,
00334 72, .5, 72.5);
00335
00336
00337 meFullRoRu_ = book2D("h2FRORu",
00338 "Full Read-out readout unit;"
00339 "iX - 40 / iEta / iX + 20;"
00340 "iY / iPhi (iPhi = 1 at phi = 0 rad);"
00341 "Event count",
00342 80, -39.5, 40.5,
00343 72, .5, 72.5);
00344
00345 meFullRoCnt_ = book1D("hFROCnt",
00346 "Number of Full-readout-flagged readout units;"
00347 "FRO RU count;Event count",
00348 300, -.5, 299.5);
00349
00350 meEbFullRoCnt_ = book1D("hEbFROCnt",
00351 "Number of EB Full-readout-flagged readout units;"
00352 "FRO RU count;Event count",
00353 200, -.5, 199.5);
00354
00355 meEeFullRoCnt_ = book1D("hEeFROCnt",
00356 "Number of EE Full-readout-flagged readout units;"
00357 "FRO RU count;Event count",
00358 200, -.5, 199.5);
00359
00360 meZs1Ru_ = book2D("h2Zs1Ru",
00361 "Readout unit with ZS-thr-1 flag;"
00362 "iX - 40 / iEta / iX + 20;"
00363 "iY0 / iPhi0 (iPhi = 1 at phi = 0 rad);"
00364 "Event count",
00365 80, -39.5, 40.5,
00366 72, .5, 72.5);
00367
00368 meForcedRu_ = book2D("h2ForcedRu",
00369 "ECAL readout unit with forced bit of SR flag on;"
00370 "iX - 40 / iEta / iX + 20;"
00371 "iY / iPhi (iPhi = 1 at phi = 0 rad);"
00372 "Event count",
00373 80, -39.5, 40.5,
00374 72, .5, 72.5);
00375
00376 meLiTtf_ = book2D("h2LiTtf",
00377 "Low interest trigger tower flags;"
00378 "iEta;"
00379 "iPhi;"
00380 "Event count",
00381 57, -28.5, 28.5,
00382 72, .5, 72.5);
00383
00384 meMiTtf_ = book2D("h2MiTtf",
00385 "Mid interest trigger tower flags;"
00386 "iEta;"
00387 "iPhi;"
00388 "Event count",
00389 57, -28.5, 28.5,
00390 72, .5, 72.5);
00391
00392 meHiTtf_ = book2D("h2HiTtf",
00393 "High interest trigger tower flags;"
00394 "iEta;"
00395 "iPhi;"
00396 "Event count",
00397 57, -28.5, 28.5,
00398 72, .5, 72.5);
00399
00400 meForcedTtf_ = book2D("h2ForcedTtf",
00401 "Trigger tower flags with forced bit set;"
00402 "iEta;"
00403 "iPhi;"
00404 "Event count",
00405 57, -28.5, 28.5,
00406 72, .5, 72.5);
00407
00408
00409 const float ebMinNoise = -1.;
00410 const float ebMaxNoise = 1.;
00411
00412 const float eeMinNoise = -1.;
00413 const float eeMaxNoise = 1.;
00414
00415 #if 0
00416 const float ebMinE = 0.;
00417 const float ebMaxE = 120.;
00418
00419 const float eeMinE = 0.;
00420 const float eeMaxE = 120.;
00421 #else
00422 const float ebMinE = ebMinNoise;
00423 const float ebMaxE = ebMaxNoise;
00424
00425 const float eeMinE = eeMinNoise;
00426 const float eeMaxE = eeMaxNoise;
00427 #endif
00428
00429
00430 const int evtMax = 500;
00431
00432 meEbRecE_ = book1D("hEbRecE",
00433 "Crystal reconstructed energy;E (GeV);Event count",
00434 100, ebMinE, ebMaxE);
00435
00436 meEbEMean_ = bookProfile("hEbEMean",
00437 "EE <E_hit>;event #;<E_hit> (GeV)",
00438 evtMax, .5, evtMax + .5);
00439
00440 meEbNoise_ = book1D("hEbNoise",
00441 "Crystal noise "
00442 "(rec E of crystal without deposited energy)"
00443 ";Rec E (GeV);Event count",
00444 100, ebMinNoise, ebMaxNoise);
00445
00446 meEbLiZsFir_ = book1D("zsEbLiFIRemu",
00447 "Emulated ouput of ZS FIR filter for EB "
00448 "low interest crystals;"
00449 "ADC count*4;"
00450 "Event count",
00451 60, -30, 30);
00452
00453 meEbHiZsFir_ = book1D("zsEbHiFIRemu",
00454 "Emulated ouput of ZS FIR filter for EB "
00455 "high interest crystals;"
00456 "ADC count*4;"
00457 "Event count",
00458 60, -30, 30);
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468 meEbSimE_ = book1D("hEbSimE", "EB hit crystal simulated energy",
00469 100, ebMinE, ebMaxE);
00470
00471 meEbRecEHitXtal_ = book1D("hEbRecEHitXtal",
00472 "EB rec energy of hit crystals",
00473 100, ebMinE, ebMaxE);
00474
00475 meEbRecVsSimE_ = book2D("hEbRecVsSimE",
00476 "Crystal simulated vs reconstructed energy;"
00477 "Esim (GeV);Erec GeV);Event count",
00478 100, ebMinE, ebMaxE,
00479 100, ebMinE, ebMaxE);
00480
00481 meEbNoZsRecVsSimE_ = book2D("hEbNoZsRecVsSimE",
00482 "Crystal no-zs simulated vs reconstructed "
00483 "energy;"
00484 "Esim (GeV);Erec GeV);Event count",
00485 100, ebMinE, ebMaxE,
00486 100, ebMinE, ebMaxE);
00487
00488 meEeRecE_ = book1D("hEeRecE",
00489 "EE crystal reconstructed energy;E (GeV);"
00490 "Event count",
00491 100, eeMinE, eeMaxE);
00492
00493 meEeEMean_ = bookProfile("hEeEMean",
00494 "<E_{EE hit}>;event;<E_{hit}> (GeV)",
00495 evtMax, .5, evtMax + .5);
00496
00497
00498 meEeNoise_ = book1D("hEeNoise",
00499 "EE crystal noise "
00500 "(rec E of crystal without deposited energy);"
00501 "E (GeV);Event count",
00502 200, eeMinNoise, eeMaxNoise);
00503
00504 meEeLiZsFir_ = book1D("zsEeLiFIRemu",
00505 "Emulated ouput of ZS FIR filter for EE "
00506 "low interest crystals;"
00507 "ADC count*4;"
00508 "Event count",
00509 60, -30, 30);
00510
00511 meEeHiZsFir_ = book1D("zsEeHiFIRemu",
00512 "Emulated ouput of ZS FIR filter for EE "
00513 "high interest crystals;"
00514 "ADC count*4;"
00515 "Event count",
00516 60, -30, 30);
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527 meEeSimE_ = book1D("hEeSimE", "EE hit crystal simulated energy",
00528 100, eeMinE, eeMaxE);
00529
00530 meEeRecEHitXtal_ = book1D("hEeRecEHitXtal",
00531 "EE rec energy of hit crystals",
00532 100, eeMinE, eeMaxE);
00533
00534 meEeRecVsSimE_ = book2D("hEeRecVsSimE",
00535 "EE crystal simulated vs reconstructed energy;"
00536 "Esim (GeV);Erec GeV);Event count",
00537 100, eeMinE, eeMaxE,
00538 100, eeMinE, eeMaxE);
00539
00540 meEeNoZsRecVsSimE_ = book2D("hEeNoZsRecVsSimE",
00541 "EE crystal no-zs simulated vs "
00542 "reconstructed "
00543 "energy;Esim (GeV);Erec GeV);Event count",
00544 100, eeMinE, eeMaxE,
00545 100, eeMinE, eeMaxE);
00546
00547 meSRFlagsConsistency_ = book2D("hSRAlgoErrorMap",
00548 "TTFlags and SR Flags mismatch;"
00549 "iX - 40 / iEta / iX + 20;"
00550 "iY / iPhi (iPhi = 1 at phi = 0 rad);"
00551 "Event count",
00552 80, -39.5, 40.5,
00553 72, .5, 72.5);
00554
00555
00556 meIncompleteFROMap_ = book2D("hIncompleteFROMap",
00557 "Incomplete full-readout-flagged readout units;"
00558 "iX - 40 / iEta / iX + 20;"
00559 "iY / iPhi (iPhi = 1 at phi = 0 rad);"
00560 "Event count",
00561 80, -39.5, 40.5,
00562 72, .5, 72.5);
00563
00564 meIncompleteFROCnt_ = book1D("hIncompleteFROCnt",
00565 "Number of incomplete full-readout-flagged "
00566 "readout units;"
00567 "Number of RUs;Event count;",
00568 200, -.5, 199.5);
00569
00570 meIncompleteFRORateMap_
00571 = bookProfile2D("hIncompleteFRORateMap",
00572 "Incomplete full-readout-flagged readout units;"
00573 "iX - 40 / iEta / iX + 20;"
00574 "iY / iPhi (iPhi = 1 at phi = 0 rad);"
00575 "Incomplete error rate",
00576 80, -39.5, 40.5,
00577 72, .5, 72.5);
00578
00579
00580 meDroppedFROMap_ = book2D("hDroppedFROMap",
00581 "Dropped full-readout-flagged readout units;"
00582 "iX - 40 / iEta / iX + 20;"
00583 "iY / iPhi (iPhi = 1 at phi = 0 rad);"
00584 "Event count",
00585 80, -39.5, 40.5,
00586 72, .5, 72.5);
00587
00588 meDroppedFROCnt_ = book1D("hDroppedFROCnt",
00589 "Number of dropped full-readout-flagged "
00590 "RU count;RU count;Event count",
00591 200, -.5, 199.5);
00592
00593 meCompleteZSCnt_ = book1D("hCompleteZsCnt",
00594 "Number of zero-suppressed-flagged RU fully "
00595 "readout;"
00596 "RU count;Event count",
00597 200, -.5, 199.5);
00598
00599 stringstream buf;
00600 buf << "Number of LI EB channels below the " << ebZsThr_/4. << " ADC count ZS threshold;"
00601 "Channel count;Event count",
00602 meEbZsErrCnt_ = book1D("hEbZsErrCnt",
00603 buf.str().c_str(),
00604 200, -.5, 199.5);
00605
00606 buf.str("");
00607 buf << "Number of LI EE channels below the " << eeZsThr_/4. << " ADC count ZS theshold;"
00608 "Channel count;Event count",
00609 meEeZsErrCnt_ = book1D("hEeZsErrCnt",
00610 buf.str().c_str(),
00611 200, -.5, 199.5);
00612
00613 meZsErrCnt_ = book1D("hZsErrCnt",
00614 "Number of LI channels below the ZS threshold;"
00615 "Channel count;Event count",
00616 200, -.5, 199.5);
00617
00618 meEbZsErrType1Cnt_ = book1D("hEbZsErrType1Cnt",
00619 "Number of EB channels below the ZS "
00620 "threshold in a LI but fully readout RU;"
00621 "Channel count;Event count;",
00622 200, -.5, 199.5);
00623
00624 meEeZsErrType1Cnt_ = book1D("hEeZsErrType1Cnt",
00625 "Number EE channels below the ZS threshold"
00626 " in a LI but fully readout RU;"
00627 "Channel count;Event count",
00628 200, -.5, 199.5);
00629
00630 meZsErrType1Cnt_ = book1D("hZsErrType1Cnt",
00631 "Number of LI channels below the ZS threshold "
00632 "in a LI but fully readout RU;"
00633 "Channel count;Event count",
00634 200, -.5, 199.5);
00635
00636
00637 meDroppedFRORateMap_
00638 = bookProfile2D("hDroppedFRORateMap",
00639 "Dropped full-readout-flagged readout units"
00640 "iX - 40 / iEta / iX + 20;"
00641 "iY / iPhi (iPhi = 1 at phi = 0 rad);"
00642 "Dropping rate",
00643 80, -39.5, 40.5,
00644 72, .5, 72.5);
00645
00646 meCompleteZSMap_ = book2D("hCompleteZSMap",
00647 "Complete zero-suppressed-flagged readout units;"
00648 "iX - 40 / iEta / iX + 20;"
00649 "iY / iPhi (iPhi = 1 at phi = 0 rad);"
00650 "Event count",
00651 80, -39.5, 40.5,
00652 72, .5, 72.5);
00653
00654 meCompleteZSRateMap_
00655 = bookProfile2D("hCompleteZSRate",
00656 "Complete zero-suppressed-flagged readout units;"
00657 "iX - 40 / iEta / iX + 20;"
00658 "iY / iPhi (iPhi = 1 at phi = 0 rad);"
00659 "Completeness rate",
00660 80, -39.5, 40.5,
00661 72, .5, 72.5);
00662
00663
00664
00665 printAvailableHists();
00666
00667
00668 stringstream s;
00669 for(set<string>::iterator it = histList_.begin();
00670 it!=histList_.end();
00671 ++it){
00672 if(*it!=string("all")
00673 && availableHistList_.find(*it)==availableHistList_.end()){
00674 s << (s.str().size()==0?"":", ") << *it;
00675 }
00676 }
00677 if(s.str().size()!=0){
00678 LogWarning("Configuration")
00679 << "Parameter 'histList' contains some unknown histogram(s). "
00680 "Check spelling. Following name were not found: "
00681 << s.str();
00682 }
00683 }
00684
00685
00686 void EcalSelectiveReadoutValidation::updateL1aRate(const edm::Event& event){
00687 const int32_t bx = event.bunchCrossing();
00688 if(bx<1 || bx > 3564) return;
00689
00690
00691 int64_t t = event.bunchCrossing() + (event.orbitNumber()-1)*3564;
00692
00693 if(t<tmin){
00694 tmin = t;
00695 l1aOfTmin = event.id().event();
00696 }
00697
00698 if(t>tmax){
00699 tmax = t;
00700 l1aOfTmax = event.id().event();
00701 }
00702 }
00703
00704 double EcalSelectiveReadoutValidation::getL1aRate() const{
00705 LogDebug("EcalSrValid") << __FILE__ << ":" << __LINE__ << ": "
00706 << "Tmax = " << tmax << " x 25ns; Tmin = " << tmin
00707 << " x 25ns; L1A(Tmax) = " << l1aOfTmax << "; L1A(Tmin) = "
00708 << l1aOfTmin << "\n";
00709 return (double)(l1aOfTmax - l1aOfTmin) / ((tmax-tmin) * 25e-9);
00710 }
00711
00712 void EcalSelectiveReadoutValidation::analyze(Event const & event,
00713 EventSetup const & es){
00714
00715 updateL1aRate(event);
00716
00717 {
00718 PgTiming t("collection readout");
00719
00720
00721 readAllCollections(event);
00722
00723 }
00724
00725 withEeSimHit_ = (eeSimHits_->size()!=0);
00726 withEbSimHit_ = (ebSimHits_->size()!=0);
00727
00728 if(ievt_<10){
00729 edm::LogInfo("EcalSrValid") << "Size of TP collection: " << tps_->size() << "\n"
00730 << "Size of EB SRF collection read from data: "
00731 << ebSrFlags_->size() << "\n"
00732 << "Size of EB SRF collection computed from data TTFs: "
00733 << ebComputedSrFlags_->size() << "\n"
00734 << "Size of EE SRF collection read from data: "
00735 << eeSrFlags_->size() << "\n"
00736 << "Size of EE SRF collection computed from data TTFs: "
00737 << eeComputedSrFlags_->size() << "\n";
00738 }
00739
00740 if(ievt_==0){
00741 selectFedsForLog();
00742 }
00743
00744
00745 setTtEtSums(es, *ebNoZsDigis_, *eeNoZsDigis_);
00746
00747 {
00748 PgTiming t("data volume analysis");
00749
00750
00751 analyzeDataVolume(event, es);
00752 }
00753
00754 {
00755 PgTiming t("EB analysis");
00756
00757
00758
00759 analyzeEB(event, es);
00760 }
00761
00762 {
00763 PgTiming t("EE analysis");
00764
00765
00766
00767 analyzeEE(event, es);
00768 }
00769
00770 fill(meFullRoCnt_, nEeFROCnt_+nEbFROCnt_);
00771 fill(meEbFullRoCnt_, nEbFROCnt_);
00772 fill(meEeFullRoCnt_, nEeFROCnt_);
00773
00774 fill(meEbZsErrCnt_, nEbZsErrors_);
00775 fill(meEeZsErrCnt_, nEeZsErrors_);
00776 fill(meZsErrCnt_, nEbZsErrors_ + nEeZsErrors_);
00777
00778 fill(meEbZsErrType1Cnt_, nEbZsErrorsType1_);
00779 fill(meEeZsErrType1Cnt_, nEeZsErrorsType1_);
00780 fill(meZsErrType1Cnt_, nEbZsErrorsType1_ + nEeZsErrorsType1_);
00781
00782 {
00783 PgTiming t("TP analysis");
00784
00785 analyzeTP(event, es);
00786 }
00787
00788
00789
00790 if(ebComputedSrFlags_->size()){
00791 compareSrfColl(event, *ebSrFlags_, *ebComputedSrFlags_);
00792 }
00793 if(eeComputedSrFlags_->size()){
00794 compareSrfColl(event, *eeSrFlags_, *eeComputedSrFlags_);
00795 }
00796 nDroppedFRO_ = 0;
00797 nIncompleteFRO_ = 0;
00798 nCompleteZS_ = 0;
00799 checkSrApplication(event, *ebSrFlags_);
00800 checkSrApplication(event, *eeSrFlags_);
00801 fill(meDroppedFROCnt_, nDroppedFRO_);
00802 fill(meIncompleteFROCnt_, nIncompleteFRO_);
00803 fill(meCompleteZSCnt_, nCompleteZS_);
00804 ++ievt_;
00805 }
00806
00807
00808 void EcalSelectiveReadoutValidation::analyzeEE(const edm::Event& event,
00809 const edm::EventSetup& es){
00810 bool eventError = false;
00811 nEeZsErrors_ = 0;
00812
00813 {
00814 PgTiming t("analyzeEE: init");
00815 for(int iZ0=0; iZ0<nEndcaps; ++iZ0){
00816 for(int iX0=0; iX0<nEeX; ++iX0){
00817 for(int iY0=0; iY0<nEeY; ++iY0){
00818 eeEnergies[iZ0][iX0][iY0].noZsRecE = -numeric_limits<double>::max();
00819 eeEnergies[iZ0][iX0][iY0].recE = -numeric_limits<double>::max();
00820 eeEnergies[iZ0][iX0][iY0].simE = 0;
00821 eeEnergies[iZ0][iX0][iY0].simHit = 0;
00822 eeEnergies[iZ0][iX0][iY0].gain12 = false;
00823 }
00824 }
00825 }
00826 }
00827
00828
00829 edm::ESHandle<CaloGeometry> geoHandle;
00830 es.get<MyCaloGeometryRecord>().get(geoHandle);
00831 const CaloSubdetectorGeometry *geometry_p
00832 = (*geoHandle).getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
00833 CaloSubdetectorGeometry const& geometry = *geometry_p;
00834
00835 {
00836 PgTiming t("analyzeEE: unsupressed digis");
00837
00838 for (unsigned int digis=0; digis<eeNoZsDigis_->size(); ++digis){
00839
00840 EEDataFrame frame = (*eeNoZsDigis_)[digis];
00841 int iX0 = iXY2cIndex(frame.id().ix());
00842 int iY0 = iXY2cIndex(frame.id().iy());
00843 int iZ0 = frame.id().zside()>0?1:0;
00844
00845 if(iX0<0 || iX0>=nEeX){
00846 edm::LogError("EcalSrValid") << "iX0 (= " << iX0 << ") is out of range ("
00847 << "[0," << nEeX -1 << "]\n";
00848 }
00849 if(iY0<0 || iY0>=nEeY){
00850 edm::LogError("EcalSrValid") << "iY0 (= " << iY0 << ") is out of range ("
00851 << "[0," << nEeY -1 << "]\n";
00852 }
00853
00854 eeEnergies[iZ0][iX0][iY0].noZsRecE = frame2Energy(frame);
00855
00856 eeEnergies[iZ0][iX0][iY0].gain12 = true;
00857 for(int i = 0; i< frame.size(); ++i){
00858 const int gain12Code = 0x1;
00859 if(frame[i].gainId()!=gain12Code) eeEnergies[iZ0][iX0][iY0].gain12 = false;
00860 }
00861
00862 const GlobalPoint xtalPos
00863 = geometry.getGeometry(frame.id())->getPosition();
00864
00865 eeEnergies[iZ0][iX0][iY0].phi = rad2deg*((double)xtalPos.phi());
00866 eeEnergies[iZ0][iX0][iY0].eta = xtalPos.eta();
00867 }
00868 }
00869
00870 {
00871 PgTiming t("analyzeEE:rec hits");
00872
00873 if(!localReco_){
00874 for(RecHitCollection::const_iterator it
00875 = eeRecHits_->begin();
00876 it != eeRecHits_->end(); ++it){
00877 const RecHit& hit = *it;
00878 int iX0 = iXY2cIndex(static_cast<const EEDetId&>(hit.id()).ix());
00879 int iY0 = iXY2cIndex(static_cast<const EEDetId&>(hit.id()).iy());
00880 int iZ0 = static_cast<const EEDetId&>(hit.id()).zside()>0?1:0;
00881
00882 if(iX0<0 || iX0>=nEeX){
00883 LogError("EcalSrValid") << "iX0 (= " << iX0 << ") is out of range ("
00884 << "[0," << nEeX -1 << "]\n";
00885 }
00886 if(iY0<0 || iY0>=nEeY){
00887 LogError("EcalSrValid") << "iY0 (= " << iY0 << ") is out of range ("
00888 << "[0," << nEeY -1 << "]\n";
00889 }
00890
00891 eeEnergies[iZ0][iX0][iY0].recE = hit.energy();
00892 }
00893 }
00894 }
00895
00896 {
00897 PgTiming t("analyzeEE:sim hits");
00898
00899 for(vector<PCaloHit>::const_iterator it = eeSimHits_->begin();
00900 it != eeSimHits_->end(); ++it){
00901 const PCaloHit& simHit = *it;
00902 EEDetId detId(simHit.id());
00903 int iX = detId.ix();
00904 int iX0 =iXY2cIndex(iX);
00905 int iY = detId.iy();
00906 int iY0 = iXY2cIndex(iY);
00907 int iZ0 = detId.zside()>0?1:0;
00908 eeEnergies[iZ0][iX0][iY0].simE += simHit.energy();
00909 ++eeEnergies[iZ0][iX0][iY0].simHit;
00910 }
00911 }
00912
00913 {
00914 PgTiming t("analyzeEE: suppressed digis");
00915
00916
00917 for(EEDigiCollection::const_iterator it = eeDigis_->begin();
00918 it != eeDigis_->end(); ++it){
00919 const EEDataFrame& frame = *it;
00920 int iX0 = iXY2cIndex(static_cast<const EEDetId&>(frame.id()).ix());
00921 int iY0 = iXY2cIndex(static_cast<const EEDetId&>(frame.id()).iy());
00922 int iZ0 = static_cast<const EEDetId&>(frame.id()).zside()>0?1:0;
00923 if(iX0<0 || iX0>=nEeX){
00924 LogError("EcalSrValid") << "iX0 (= " << iX0 << ") is out of range ("
00925 << "[0," << nEeX -1 << "]\n";
00926 }
00927 if(iY0<0 || iY0>=nEeY){
00928 LogError("EcalSrValid") << "iY0 (= " << iY0 << ") is out of range ("
00929 << "[0," << nEeY -1 << "]\n";
00930 }
00931
00932 if(localReco_){
00933 eeEnergies[iZ0][iX0][iY0].recE = frame2Energy(frame);
00934 }
00935
00936 eeEnergies[iZ0][iX0][iY0].gain12 = true;
00937 for(int i = 0; i< frame.size(); ++i){
00938 const int gain12Code = 0x1;
00939 if(frame[i].gainId()!=gain12Code){
00940 eeEnergies[iZ0][iX0][iY0].gain12 = false;
00941 }
00942 }
00943
00944 fill(meChOcc_, xtalGraphX(frame.id()), xtalGraphY(frame.id()));
00945
00946 EESrFlagCollection::const_iterator srf
00947 = eeSrFlags_->find(readOutUnitOf(frame.id()));
00948
00949 bool highInterest = false;
00950
00951
00952 if(srf==eeSrFlags_->end()) continue;
00953
00954 if(srf!=eeSrFlags_->end()){
00955 highInterest = ((srf->value() & ~EcalSrFlag::SRF_FORCED_MASK)
00956 == EcalSrFlag::SRF_FULL);
00957 }
00958
00959 if(highInterest){
00960 fill(meEeHiZsFir_, dccZsFIR(frame, firWeights_, firstFIRSample_, 0));
00961 } else{
00962 int v = dccZsFIR(frame, firWeights_, firstFIRSample_, 0);
00963 fill(meEeLiZsFir_, v);
00964 if(v < eeZsThr_){
00965 eventError = true;
00966 ++nEeZsErrors_;
00967 pair<int,int> ru = dccCh(frame.id());
00968 if(isRuComplete_[ru.first][ru.second-1]) ++nEeZsErrorsType1_;
00969 if(nEeZsErrors_ < 3){
00970 srApplicationErrorLog_ << event.id() << ", "
00971 << "RU " << frame.id() << ", "
00972 << "DCC " << ru.first
00973 << " Ch : " << ru.second << ": "
00974 << "LI channel under ZS threshold.\n";
00975 }
00976 if(nEeZsErrors_==3){
00977 srApplicationErrorLog_ << event.id() << ": "
00978 << "more ZS errors for this event...\n";
00979 }
00980 }
00981 }
00982 }
00983 }
00984
00985 {
00986 PgTiming t("analyzeEE: energies");
00987
00988 for(int iZ0=0; iZ0<nEndcaps; ++iZ0){
00989 for(int iX0=0; iX0<nEeX; ++iX0){
00990 for(int iY0=0; iY0<nEeY; ++iY0){
00991 double recE = eeEnergies[iZ0][iX0][iY0].recE;
00992 if(recE==-numeric_limits<double>::max()) continue;
00993 fill(meEeRecE_, eeEnergies[iZ0][iX0][iY0].recE);
00994
00995 fill(meEeEMean_, ievt_+1,
00996 eeEnergies[iZ0][iX0][iY0].recE);
00997
00998 if(withEeSimHit_){
00999 if(!eeEnergies[iZ0][iX0][iY0].simHit){
01000 fill(meEeNoise_, eeEnergies[iZ0][iX0][iY0].noZsRecE);
01001 } else{
01002 fill(meEeSimE_, eeEnergies[iZ0][iX0][iY0].simE);
01003 fill(meEeRecEHitXtal_, eeEnergies[iZ0][iX0][iY0].recE);
01004 }
01005 fill(meEeRecVsSimE_, eeEnergies[iZ0][iX0][iY0].simE,
01006 eeEnergies[iZ0][iX0][iY0].recE);
01007 fill(meEeNoZsRecVsSimE_, eeEnergies[iZ0][iX0][iY0].simE,
01008 eeEnergies[iZ0][iX0][iY0].noZsRecE);
01009 }
01010 }
01011 }
01012 }
01013 }
01014
01015 {
01016 PgTiming t("analyzeEE: RU");
01017
01018 nEeFROCnt_ = 0;
01019 char eeSrfMark[2][100][100];
01020 bzero(eeSrfMark, sizeof(eeSrfMark));
01021
01022 for(EESrFlagCollection::const_iterator it = eeSrFlags_->begin();
01023 it != eeSrFlags_->end(); ++it){
01024 const EESrFlag& srf = *it;
01025 int iX = srf.id().ix();
01026 int iY = srf.id().iy();
01027 int iZ = srf.id().zside();
01028 if(iX<1 || iY > 100) throw cms::Exception("EcalSelectiveReadoutValidation")
01029 << "Found an endcap SRF with an invalid det ID: " << srf.id() << ".\n";
01030 ++eeSrfMark[iZ>0?1:0][iX-1][iY-1];
01031 if(eeSrfMark[iZ>0?1:0][iX-1][iY-1] > 1) throw cms::Exception("EcalSelectiveReadoutValidation")
01032 << "Duplicate SRF for supercrystal " << srf.id() << ".\n";
01033 int flag = srf.value() & ~EcalSrFlag::SRF_FORCED_MASK;
01034 if(flag == EcalSrFlag::SRF_ZS1){
01035 fill(meZs1Ru_, ruGraphX(srf.id()), ruGraphY(srf.id()));
01036 }
01037
01038 if(flag == EcalSrFlag::SRF_FULL){
01039 fill(meFullRoRu_, ruGraphX(srf.id()), ruGraphY(srf.id()));
01040 ++nEeFROCnt_;
01041 }
01042
01043 if(srf.value() & EcalSrFlag::SRF_FORCED_MASK){
01044 fill(meForcedRu_, ruGraphX(srf.id()), ruGraphY(srf.id()));
01045 }
01046 }
01047 }
01048
01049 {
01050 PgTiming t("analyzeEE: SR appli error log");
01051
01052 if(eventError) srApplicationErrorLog_ << event.id()
01053 << ": " << nEeZsErrors_
01054 << " ZS-flagged EE channels under "
01055 "the ZS threshold, whose " << nEeZsErrorsType1_
01056 << " in a complete RU.\n";
01057 }
01058 }
01059
01060 void
01061 EcalSelectiveReadoutValidation::analyzeEB(const edm::Event& event,
01062 const edm::EventSetup& es){
01063
01064 bool eventError = false;
01065 nEbZsErrors_ = 0;
01066 vector<pair<int,int> > xtalEtaPhi;
01067
01068 {
01069 PgTiming t("analyzeEB: init");
01070
01071 xtalEtaPhi.reserve(nEbPhi*nEbEta);
01072 for(int iEta0=0; iEta0<nEbEta; ++iEta0){
01073 for(int iPhi0=0; iPhi0<nEbPhi; ++iPhi0){
01074 ebEnergies[iEta0][iPhi0].noZsRecE = -numeric_limits<double>::max();
01075 ebEnergies[iEta0][iPhi0].recE = -numeric_limits<double>::max();
01076 ebEnergies[iEta0][iPhi0].simE = 0;
01077 ebEnergies[iEta0][iPhi0].simHit = 0;
01078 ebEnergies[iEta0][iPhi0].gain12 = false;
01079 xtalEtaPhi.push_back(pair<int,int>(iEta0, iPhi0));
01080 }
01081 }
01082 }
01083
01084
01085 edm::ESHandle<CaloGeometry> geoHandle;
01086
01087 PgTiming t1("analyzeEB: geomRetrieval");
01088 es.get<MyCaloGeometryRecord>().get(geoHandle);
01089 const CaloSubdetectorGeometry *geometry_p
01090 = (*geoHandle).getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
01091 CaloSubdetectorGeometry const& geometry = *geometry_p;
01092 t1.stop();
01093
01094
01095 {
01096 PgTiming t("analyzeEB: unsuppressed digi loop");
01097
01098 for(EBDigiCollection::const_iterator it = ebNoZsDigis_->begin();
01099 it != ebNoZsDigis_->end(); ++it){
01100 const EBDataFrame& frame = *it;
01101 int iEta0 = iEta2cIndex(static_cast<const EBDetId&>(frame.id()).ieta());
01102 int iPhi0 = iPhi2cIndex(static_cast<const EBDetId&>(frame.id()).iphi());
01103 if(iEta0<0 || iEta0>=nEbEta){
01104 stringstream s;
01105 s << "EcalSelectiveReadoutValidation: "
01106 << "iEta0 (= " << iEta0 << ") is out of range ("
01107 << "[0," << nEbEta -1 << "]\n";
01108 throw cms::Exception(s.str());
01109 }
01110 if(iPhi0<0 || iPhi0>=nEbPhi){
01111 stringstream s;
01112 s << "EcalSelectiveReadoutValidation: "
01113 << "iPhi0 (= " << iPhi0 << ") is out of range ("
01114 << "[0," << nEbPhi -1 << "]\n";
01115 throw cms::Exception(s.str());
01116 }
01117
01118 ebEnergies[iEta0][iPhi0].noZsRecE = frame2Energy(frame);
01119 ebEnergies[iEta0][iPhi0].gain12 = true;
01120 for(int i = 0; i< frame.size(); ++i){
01121 const int gain12Code = 0x1;
01122 if(frame[i].gainId()!=gain12Code) ebEnergies[iEta0][iPhi0].gain12 = false;
01123 }
01124
01125 const GlobalPoint xtalPos
01126 = geometry.getGeometry(frame.id())->getPosition();
01127
01128 ebEnergies[iEta0][iPhi0].phi = rad2deg*((double)xtalPos.phi());
01129 ebEnergies[iEta0][iPhi0].eta = xtalPos.eta();
01130 }
01131 }
01132
01133
01134 {
01135 PgTiming t("analyzeEB: simHit loop");
01136
01137 for(vector<PCaloHit>::const_iterator it = ebSimHits_->begin();
01138 it != ebSimHits_->end(); ++it){
01139 const PCaloHit& simHit = *it;
01140 EBDetId detId(simHit.id());
01141 int iEta = detId.ieta();
01142 int iEta0 =iEta2cIndex(iEta);
01143 int iPhi = detId.iphi();
01144 int iPhi0 = iPhi2cIndex(iPhi);
01145 ebEnergies[iEta0][iPhi0].simE += simHit.energy();
01146 ++ebEnergies[iEta0][iPhi0].simHit;
01147 }
01148 }
01149
01150 bool crystalShot[nEbEta][nEbPhi];
01151 {
01152 PgTiming t("analyzeEB: suppressed digi loop init");
01153
01154 for(int iEta0=0; iEta0<nEbEta; ++iEta0){
01155 for(int iPhi0=0; iPhi0<nEbPhi; ++iPhi0){
01156 crystalShot[iEta0][iPhi0] = false;
01157 }
01158 }
01159 }
01160
01161 int nEbDigi = 0;
01162
01163 {
01164 PgTiming t("analyzeEB: suppressed digi loop");
01165
01166 for(EBDigiCollection::const_iterator it = ebDigis_->begin();
01167 it != ebDigis_->end(); ++it){
01168 ++nEbDigi;
01169 const EBDataFrame& frame = *it;
01170 int iEta = static_cast<const EBDetId&>(frame.id()).ieta();
01171 int iPhi = static_cast<const EBDetId&>(frame.id()).iphi();
01172 int iEta0 = iEta2cIndex(iEta);
01173 int iPhi0 = iPhi2cIndex(iPhi);
01174 if(iEta0<0 || iEta0>=nEbEta){
01175 throw (cms::Exception("EcalSelectiveReadoutValidation")
01176 << "iEta0 (= " << iEta0 << ") is out of range ("
01177 << "[0," << nEbEta -1 << "]");
01178 }
01179 if(iPhi0<0 || iPhi0>=nEbPhi){
01180 throw (cms::Exception("EcalSelectiveReadoutValidation")
01181 << "iPhi0 (= " << iPhi0 << ") is out of range ("
01182 << "[0," << nEbPhi -1 << "]");
01183 }
01184 assert(iEta0>=0 && iEta0<nEbEta);
01185 assert(iPhi0>=0 && iPhi0<nEbPhi);
01186 if(!crystalShot[iEta0][iPhi0]){
01187 crystalShot[iEta0][iPhi0] = true;
01188 } else{
01189 cout << "Error: several digi for same crystal!";
01190 abort();
01191 }
01192 if(localReco_){
01193 ebEnergies[iEta0][iPhi0].recE = frame2Energy(frame);
01194 }
01195
01196 ebEnergies[iEta0][iPhi0].gain12 = true;
01197 for(int i = 0; i< frame.size(); ++i){
01198 const int gain12Code = 0x1;
01199 if(frame[i].gainId()!=gain12Code){
01200 ebEnergies[iEta0][iPhi0].gain12 = false;
01201 }
01202 }
01203
01204 fill(meChOcc_, xtalGraphX(frame.id()), xtalGraphY(frame.id()));
01205 EBSrFlagCollection::const_iterator srf
01206 = ebSrFlags_->find(readOutUnitOf(frame.id()));
01207
01208 bool highInterest = false;
01209
01210
01211
01212
01213
01214
01215 if(srf != ebSrFlags_->end()){
01216 highInterest = ((srf->value() & ~EcalSrFlag::SRF_FORCED_MASK)
01217 == EcalSrFlag::SRF_FULL);
01218 }
01219
01220
01221 if(highInterest){
01222 fill(meEbHiZsFir_, dccZsFIR(frame, firWeights_, firstFIRSample_, 0));
01223 } else{
01224 int v = dccZsFIR(frame, firWeights_, firstFIRSample_, 0);
01225 fill(meEbLiZsFir_, v);
01226 if(v < ebZsThr_){
01227 eventError = true;
01228 ++nEbZsErrors_;
01229 pair<int,int> ru = dccCh(frame.id());
01230 if(isRuComplete_[ru.first][ru.second-1]) ++nEbZsErrorsType1_;
01231 if(nEbZsErrors_ < 3){
01232 srApplicationErrorLog_ << event.id() << ", "
01233 << "RU " << frame.id() << ", "
01234 << "DCC " << ru.first
01235 << " Ch : " << ru.second << ": "
01236 << "LI channel under ZS threshold.\n";
01237 }
01238 if(nEbZsErrors_==3){
01239 srApplicationErrorLog_ << event.id() << ": "
01240 << "more ZS errors for this event...\n";
01241 }
01242 }
01243 }
01244 }
01245 }
01246
01247
01248 {
01249 PgTiming t("analyzeEB: rec hit loop");
01250
01251 if(!localReco_){
01252 for(RecHitCollection::const_iterator it
01253 = ebRecHits_->begin();
01254 it != ebRecHits_->end(); ++it){
01255 ++nEbDigi;
01256 const RecHit& hit = *it;
01257 int iEta = static_cast<const EBDetId&>(hit.id()).ieta();
01258 int iPhi = static_cast<const EBDetId&>(hit.id()).iphi();
01259 int iEta0 = iEta2cIndex(iEta);
01260 int iPhi0 = iPhi2cIndex(iPhi);
01261 if(iEta0<0 || iEta0>=nEbEta){
01262 LogError("EcalSrValid") << "iEta0 (= " << iEta0 << ") is out of range ("
01263 << "[0," << nEbEta -1 << "]\n";
01264 }
01265 if(iPhi0<0 || iPhi0>=nEbPhi){
01266 LogError("EcalSrValid") << "iPhi0 (= " << iPhi0 << ") is out of range ("
01267 << "[0," << nEbPhi -1 << "]\n";
01268 }
01269 ebEnergies[iEta0][iPhi0].recE = hit.energy();
01270 }
01271 }
01272 }
01273
01274
01275
01276
01277
01278
01279
01280
01281
01282
01283
01284
01285 {
01286 PgTiming t("analyzeEB: loop on energies");
01287
01288 for(unsigned int i=0; i<xtalEtaPhi.size(); ++i){
01289 int iEta0 = xtalEtaPhi[i].first;
01290 int iPhi0= xtalEtaPhi[i].second;
01291 energiesEb_t& energies = ebEnergies[iEta0][iPhi0];
01292
01293 double recE = energies.recE;
01294 if(recE!=-numeric_limits<double>::max()){
01295 fill(meEbRecE_, ebEnergies[iEta0][iPhi0].recE);
01296 fill(meEbEMean_, ievt_+1, recE);
01297 }
01298
01299 if(withEbSimHit_){
01300 if(!energies.simHit){
01301 fill(meEbNoise_, energies.noZsRecE);
01302 } else{
01303 fill(meEbSimE_, energies.simE);
01304 fill(meEbRecEHitXtal_, energies.recE);
01305 }
01306 fill(meEbRecVsSimE_, energies.simE, energies.recE);
01307 fill(meEbNoZsRecVsSimE_, energies.simE, energies.noZsRecE);
01308 }
01309 }
01310 }
01311
01312 {
01313 PgTiming t("analyzeEB: SRF");
01314
01315 nEbFROCnt_ = 0;
01316 char ebSrfMark[2][17][72];
01317 bzero(ebSrfMark, sizeof(ebSrfMark));
01318
01319 for(EBSrFlagCollection::const_iterator it = ebSrFlags_->begin();
01320 it != ebSrFlags_->end(); ++it){
01321 const EBSrFlag& srf = *it;
01322 int iEtaAbs = srf.id().ietaAbs();
01323 int iPhi = srf.id().iphi();
01324 int iZ = srf.id().zside();
01325
01326
01327
01328
01329 if(iEtaAbs < 1 || iEtaAbs > 17
01330 || iPhi < 1 || iPhi > 72) throw cms::Exception("EcalSelectiveReadoutValidation")
01331 << "Found a barrel SRF with an invalid det ID: " << srf.id() << ".\n";
01332 ++ebSrfMark[iZ>0?1:0][iEtaAbs-1][iPhi-1];
01333 if(ebSrfMark[iZ>0?1:0][iEtaAbs-1][iPhi-1] > 1) throw cms::Exception("EcalSelectiveReadoutValidation")
01334 << "Duplicate SRF for RU " << srf.id() << ".\n";
01335 int flag = srf.value() & ~EcalSrFlag::SRF_FORCED_MASK;
01336 if(flag == EcalSrFlag::SRF_ZS1){
01337 fill(meZs1Ru_, ruGraphX(srf.id()), ruGraphY(srf.id()));
01338 }
01339 if(flag == EcalSrFlag::SRF_FULL){
01340 fill(meFullRoRu_, ruGraphX(srf.id()), ruGraphY(srf.id()));
01341 ++nEbFROCnt_;
01342 }
01343 if(srf.value() & EcalSrFlag::SRF_FORCED_MASK){
01344 fill(meForcedRu_, ruGraphX(srf.id()), ruGraphY(srf.id()));
01345 }
01346 }
01347 }
01348
01349 {
01350 PgTiming t("analyzeEB: logSRerror");
01351
01352 if(eventError) srApplicationErrorLog_ << event.id()
01353 << ": " << nEbZsErrors_
01354 << " ZS-flagged EB channels under "
01355 "the ZS threshold, whose " << nEbZsErrorsType1_
01356 << " in a complete RU.\n";
01357 }
01358 }
01359
01360 EcalSelectiveReadoutValidation::~EcalSelectiveReadoutValidation(){
01361 }
01362
01363 void EcalSelectiveReadoutValidation::beginRun(const edm::Run& r, const edm::EventSetup& es){
01364
01365 edm::ESHandle<EcalTrigTowerConstituentsMap> hTriggerTowerMap;
01366 es.get<IdealGeometryRecord>().get(hTriggerTowerMap);
01367 triggerTowerMap_ = hTriggerTowerMap.product();
01368
01369
01370 ESHandle< EcalElectronicsMapping > ecalmapping;
01371 es.get< EcalMappingRcd >().get(ecalmapping);
01372 elecMap_ = ecalmapping.product();
01373
01374 initAsciiFile();
01375 }
01376
01377 void EcalSelectiveReadoutValidation::endRun(const edm::Run& r, const edm::EventSetup& es){
01378 meL1aRate_->Fill(getL1aRate());
01379 if(useEventRate_) normalizeHists(ievt_);
01380 if(outputFile_.size()!=0) dbe_->save(outputFile_);
01381 }
01382
01383 void
01384 EcalSelectiveReadoutValidation::analyzeTP(edm::Event const & event,
01385 edm::EventSetup const & es){
01386 EcalTPGScale ecalScale;
01387 #if (CMSSW_COMPAT_VERSION>=210)
01388 ecalScale.setEventSetup(es) ;
01389 #endif
01390
01391
01392
01393
01394 for(EcalTrigPrimDigiCollection::const_iterator it = tps_->begin();
01395 it != tps_->end(); ++it){
01396
01397
01398
01399
01400
01401
01402
01403
01404
01405
01406
01407
01408 double tpEt;
01409 if(tpInGeV_){
01410 #if (CMSSW_COMPAT_VERSION<210)
01411 tpEt = ecalScale.getTPGInGeV(es, *it);
01412 #else
01413 tpEt = ecalScale.getTPGInGeV(it->compressedEt(), it->id()) ;
01414 #endif
01415 } else{
01416 tpEt = it->compressedEt();
01417 }
01418 int iEta = it->id().ieta();
01419 int iEta0 = iTtEta2cIndex(iEta);
01420 int iPhi = it->id().iphi();
01421 int iPhi0 = iTtEta2cIndex(iPhi);
01422 double etSum = ttEtSums[iEta0][iPhi0];
01423 fill(meTp_, tpEt);
01424 fill(meTpVsEtSum_, etSum, tpEt);
01425 fill(meTtf_, it->ttFlag());
01426 if((it->ttFlag() & 0x3) == 0){
01427 fill(meLiTtf_, iEta, iPhi);
01428 }
01429 if((it->ttFlag() & 0x3) == 1){
01430 fill(meMiTtf_, iEta, iPhi);
01431 }
01432 if((it->ttFlag() & 0x3) == 3){
01433 fill(meHiTtf_, iEta, iPhi);
01434 }
01435 if((it->ttFlag() & 0x4)){
01436 fill(meForcedTtf_, iEta, iPhi);
01437 }
01438
01439 fill(meTtfVsTp_, tpEt, it->ttFlag());
01440 fill(meTtfVsEtSum_, etSum, it->ttFlag());
01441 fill(meTpMap_, iEta, iPhi, tpEt, 1.);
01442 }
01443 }
01444
01445 void EcalSelectiveReadoutValidation::analyzeDataVolume(const Event& e,
01446 const EventSetup& es){
01447
01448 anaDigiInit();
01449
01450
01451
01452 for(int iDcc = minDccId_; iDcc <= maxDccId_; ++iDcc){
01453 for(int iCh = 1; iCh < nDccRus_[iDcc-minDccId_]; ++iCh){
01454 isRuComplete_[iDcc-minDccId_][iCh-1]
01455 = (nPerRu_[iDcc-minDccId_][iCh-1]==getCrystalCount(iDcc, iCh));
01456 }
01457 }
01458
01459
01460
01461 for (unsigned int digis=0; digis<ebDigis_->size(); ++digis){
01462 EBDataFrame ebdf = (*ebDigis_)[digis];
01463 anaDigi(ebdf, *ebSrFlags_);
01464 }
01465
01466
01467 for (unsigned int digis=0; digis<eeDigis_->size(); ++digis){
01468 EEDataFrame eedf = (*eeDigis_)[digis];
01469 anaDigi(eedf, *eeSrFlags_);
01470 }
01471
01472
01473 for(unsigned iDcc0 = 0; iDcc0 < nDccs_; ++iDcc0){
01474 fill(meDccVol_, iDcc0+1, getDccEventSize(iDcc0, nPerDcc_[iDcc0])/kByte_);
01475 fill(meDccLiVol_, iDcc0+1,
01476 getDccSrDependentPayload(iDcc0, nLiRuPerDcc_[iDcc0],
01477 nLiPerDcc_[iDcc0])/kByte_);
01478 fill(meDccHiVol_, iDcc0+1,
01479 getDccSrDependentPayload(iDcc0, nHiRuPerDcc_[iDcc0],
01480 nHiPerDcc_[iDcc0])/kByte_);
01481 const FEDRawDataCollection& raw = *fedRaw_;
01482 fill(meDccVolFromData_, iDcc0+1,
01483 ((double)raw.FEDData(601+iDcc0).size())/kByte_);
01484 }
01485
01486
01487
01488 double a = nEbLI_*getBytesPerCrystal()/kByte_;
01489 fill(meVolBLI_, a);
01490 double b = nEeLI_*getBytesPerCrystal()/kByte_;
01491 fill(meVolELI_, b);
01492 fill(meVolLI_, a+b);
01493
01494
01495 a = nEbHI_*getBytesPerCrystal()/kByte_;
01496 fill(meVolBHI_, a);
01497 b = nEeHI_*getBytesPerCrystal()/kByte_;
01498 fill(meVolEHI_, b);
01499 fill(meVolHI_, a+b);
01500
01501
01502 a = getEbEventSize(nEb_)/kByte_;
01503 fill(meVolB_, a);
01504 b = getEeEventSize(nEe_)/kByte_;
01505 fill(meVolE_, b);
01506 fill(meVol_, a+b);
01507 }
01508
01509
01510 template<class T, class U>
01511 void EcalSelectiveReadoutValidation::anaDigi(const T& frame,
01512 const U& srFlagColl){
01513 const DetId& xtalId = frame.id();
01514 typedef typename U::key_type RuDetId;
01515 const RuDetId& ruId = readOutUnitOf(frame.id());
01516 typename U::const_iterator srf = srFlagColl.find(ruId);
01517
01518 bool highInterest = false;
01519 int flag = 0;
01520
01521 if(srf != srFlagColl.end()){
01522
01523
01524
01525
01526 flag = srf->value() & ~EcalSrFlag::SRF_FORCED_MASK;
01527
01528 highInterest = (flag == EcalSrFlag::SRF_FULL);
01529
01530 }
01531
01532 bool barrel = (xtalId.subdetId()==EcalBarrel);
01533
01534 pair<int,int> ch = dccCh(xtalId);
01535
01536 if(barrel){
01537 ++nEb_;
01538 if(highInterest){
01539 ++nEbHI_;
01540 } else{
01541 ++nEbLI_;
01542 }
01543 int iEta0 = iEta2cIndex(static_cast<const EBDetId&>(xtalId).ieta());
01544 int iPhi0 = iPhi2cIndex(static_cast<const EBDetId&>(xtalId).iphi());
01545 if(!ebRuActive_[iEta0/ebTtEdge][iPhi0/ebTtEdge]){
01546 ++nRuPerDcc_[ch.first-minDccId_];
01547 if(highInterest){
01548
01549 ++nHiRuPerDcc_[ch.first-minDccId_];
01550 } else{
01551 ++nLiRuPerDcc_[ch.first-minDccId_];
01552 }
01553
01554
01555
01556
01557
01558
01559 ebRuActive_[iEta0/ebTtEdge][iPhi0/ebTtEdge] = true;
01560 }
01561 } else{
01562 ++nEe_;
01563 if(highInterest){
01564 ++nEeHI_;
01565 } else{
01566 ++nEeLI_;
01567 }
01568 int iX0 = iXY2cIndex(static_cast<const EEDetId&>(frame.id()).ix());
01569 int iY0 = iXY2cIndex(static_cast<const EEDetId&>(frame.id()).iy());
01570 int iZ0 = static_cast<const EEDetId&>(frame.id()).zside()>0?1:0;
01571
01572 if(!eeRuActive_[iZ0][iX0/scEdge][iY0/scEdge]){
01573 ++nRuPerDcc_[ch.first-minDccId_];
01574 if(highInterest){
01575
01576 ++nHiRuPerDcc_[ch.first-minDccId_];
01577 } else{
01578 ++nLiRuPerDcc_[ch.first-minDccId_];
01579 }
01580
01581
01582
01583
01584
01585
01586 eeRuActive_[iZ0][iX0/scEdge][iY0/scEdge] = true;
01587 }
01588 }
01589
01590 if(ch.second < 1 || ch.second > 68){
01591 throw cms::Exception("EcalSelectiveReadoutValidation")
01592 << "Error in DCC channel retrieval for crystal with detId "
01593 << xtalId.rawId() << "DCC channel out of allowed range [1..68]\n";
01594 }
01595 ++nPerDcc_[ch.first-minDccId_];
01596 ++nPerRu_[ch.first-minDccId_][ch.second-1];
01597 if(highInterest){
01598 ++nHiPerDcc_[ch.first-minDccId_];
01599 } else{
01600 ++nLiPerDcc_[ch.first-minDccId_];
01601 }
01602 }
01603
01604 void EcalSelectiveReadoutValidation::anaDigiInit(){
01605 nEb_ = 0;
01606 nEe_ = 0;
01607 nEeLI_ = 0;
01608 nEeHI_ = 0;
01609 nEbLI_ = 0;
01610 nEbHI_ = 0;
01611 bzero(nPerDcc_, sizeof(nPerDcc_));
01612 bzero(nLiPerDcc_, sizeof(nLiPerDcc_));
01613 bzero(nHiPerDcc_, sizeof(nHiPerDcc_));
01614 bzero(nRuPerDcc_, sizeof(nRuPerDcc_));
01615 bzero(ebRuActive_, sizeof(ebRuActive_));
01616 bzero(eeRuActive_, sizeof(eeRuActive_));
01617 bzero(nPerRu_, sizeof(nPerRu_));
01618 bzero(nLiRuPerDcc_, sizeof(nLiRuPerDcc_));
01619 bzero(nHiRuPerDcc_, sizeof(nHiRuPerDcc_));
01620 }
01621
01622 double EcalSelectiveReadoutValidation::frame2Energy(const EcalDataFrame& frame) const{
01623 static bool firstCall = true;
01624 if(firstCall){
01625 stringstream buf;
01626 buf << "Weights:";
01627 for(unsigned i=0; i<weights_.size();++i){
01628 buf << "\t" << weights_[i];
01629 }
01630 edm::LogInfo("EcalSrValid") << buf.str() << "\n";
01631 firstCall = false;
01632 }
01633 double adc2GeV = 0.;
01634
01635 if(typeid(EBDataFrame)==typeid(frame)){
01636 adc2GeV = .035;
01637 } else if(typeid(EEDataFrame)==typeid(frame)){
01638 adc2GeV = 0.06;
01639 } else{
01640 assert(false);
01641 }
01642
01643 double acc = 0;
01644
01645 const int n = min(frame.size(), (int)weights_.size());
01646
01647 double gainInv[] = {12., 1., 6., 12.};
01648
01649 for(int i=0; i < n; ++i){
01650 acc += weights_[i]*frame[i].adc()*gainInv[frame[i].gainId()]*adc2GeV;
01651 }
01652 return acc;
01653 }
01654
01655 int EcalSelectiveReadoutValidation::getRuCount(int iDcc0) const{
01656
01657
01658
01659
01660
01661
01662
01663
01664
01665 return nRuPerDcc_[iDcc0];
01666 }
01667
01668 pair<int,int> EcalSelectiveReadoutValidation::dccCh(const DetId& detId) const{
01669 if(detId.det()!=DetId::Ecal){
01670 throw cms::Exception("InvalidParameter")
01671 << "Wrong type of DetId passed to the "
01672 "EcalSelectiveReadoutValidation::dccCh(const DetId&). "
01673 "An ECAL DetId was expected.\n";
01674 }
01675
01676 DetId xtalId;
01677 switch(detId.subdetId()){
01678 case EcalTriggerTower:
01679 {
01680 const EcalTrigTowerDetId tt = detId;
01681
01682
01683
01684 const int iTtPhi0 = iTtPhi2cIndex(tt.iphi());
01685 const int iTtEta0 = iTtEta2cIndex(tt.ieta());
01686 const int oneXtalPhi0 = iTtPhi0 * 5;
01687 const int oneXtalEta0 = (iTtEta0 - nOneEeTtEta) * 5;
01688
01689 xtalId = EBDetId(cIndex2iEta(oneXtalEta0),
01690 cIndex2iPhi(oneXtalPhi0));
01691 }
01692 break;
01693 case EcalEndcap:
01694 if(detId.rawId() & 0x8000){
01695 return elecMap_->getDCCandSC(EcalScDetId(detId));
01696
01697
01698
01699
01700
01701 } else {
01702 xtalId = detId;
01703 }
01704 break;
01705 case EcalBarrel:
01706 xtalId = detId;
01707 break;
01708 default:
01709 throw cms::Exception("InvalidParameter")
01710 << "Wrong type of DetId passed to the method "
01711 "EcalSelectiveReadoutValidation::dccCh(const DetId&). "
01712 "A valid EcalTriggerTower, EcalBarrel or EcalEndcap DetId was expected. "
01713 "detid = " << xtalId.rawId() << ".\n";
01714 }
01715
01716 const EcalElectronicsId& EcalElecId = elecMap_->getElectronicsId(xtalId);
01717
01718 pair<int,int> result;
01719 result.first = EcalElecId.dccId();
01720
01721 if(result.first < minDccId_ || result.second > maxDccId_){
01722 throw cms::Exception("OutOfRange")
01723 << "Got an invalid DCC ID, DCCID = " << result.first
01724 << " for DetId 0x" << hex << detId.rawId()
01725 << " and 0x" << xtalId.rawId() << dec << "\n";
01726 }
01727
01728 result.second = EcalElecId.towerId();
01729
01730 if(result.second < 1 || result.second > 68){
01731 throw cms::Exception("OutOfRange")
01732 << "Got an invalid DCC channel ID, DCC_CH = " << result.second
01733 << " for DetId 0x" << hex << detId.rawId()
01734 << " and 0x" << xtalId.rawId() << dec << "\n";
01735 }
01736
01737 return result;
01738 }
01739
01740 EcalScDetId
01741 EcalSelectiveReadoutValidation::superCrystalOf(const EEDetId& xtalId) const
01742 {
01743
01744 const int scEdge = 5;
01745 EcalScDetId id = EcalScDetId((xtalId.ix()-1)/scEdge+1,
01746 (xtalId.iy()-1)/scEdge+1,
01747 xtalId.zside());
01748 return id;
01749
01750
01751
01752
01753
01754
01755
01756
01757
01758
01759
01760
01761
01762
01763
01764
01765
01766
01767
01768
01769
01770
01771
01772
01773
01774
01775
01776 }
01777
01778
01779 EcalTrigTowerDetId
01780 EcalSelectiveReadoutValidation::readOutUnitOf(const EBDetId& xtalId) const{
01781 return triggerTowerMap_->towerOf(xtalId);
01782 }
01783
01784 EcalScDetId
01785 EcalSelectiveReadoutValidation::readOutUnitOf(const EEDetId& xtalId) const{
01786
01787 const EcalElectronicsId& EcalElecId = elecMap_->getElectronicsId(xtalId);
01788 int iDCC= EcalElecId.dccId();
01789 int iDccChan = EcalElecId.towerId();
01790 const bool ignoreSingle = true;
01791 const vector<EcalScDetId> id = elecMap_->getEcalScDetId(iDCC, iDccChan, ignoreSingle);
01792 return id.size()>0?id[0]:EcalScDetId();
01793 }
01794
01795 void
01796 EcalSelectiveReadoutValidation::setTtEtSums(const edm::EventSetup& es,
01797 const EBDigiCollection& ebDigis,
01798 const EEDigiCollection& eeDigis){
01799
01800 static const CaloSubdetectorGeometry* eeGeometry = 0;
01801 static const CaloSubdetectorGeometry* ebGeometry = 0;
01802 if(eeGeometry==0 || ebGeometry==0){
01803 edm::ESHandle<CaloGeometry> geoHandle;
01804 es.get<MyCaloGeometryRecord>().get(geoHandle);
01805 eeGeometry
01806 = (*geoHandle).getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
01807 ebGeometry
01808 = (*geoHandle).getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
01809 }
01810
01811
01812 for(int iEta0 = 0; iEta0 < nTtEta; ++iEta0){
01813 for(int iPhi0 = 0; iPhi0 < nTtPhi; ++iPhi0){
01814 ttEtSums[iEta0][iPhi0] = 0.;
01815 }
01816 }
01817
01818 for(EBDigiCollection::const_iterator it = ebDigis_->begin();
01819 it != ebDigis_->end(); ++it){
01820 const EBDataFrame& frame = *it;
01821 const EcalTrigTowerDetId& ttId = triggerTowerMap_->towerOf(frame.id());
01822
01823
01824
01825
01826 const int iTtEta0 = iTtEta2cIndex(ttId.ieta());
01827 const int iTtPhi0 = iTtPhi2cIndex(ttId.iphi());
01828 double theta = ebGeometry->getGeometry(frame.id())->getPosition().theta();
01829 double e = frame2EnergyForTp(frame);
01830 if((frame2EnergyForTp(frame,-1) < e) && (frame2EnergyForTp(frame, 1) < e)){
01831 ttEtSums[iTtEta0][iTtPhi0] += e*sin(theta);
01832 }
01833 }
01834
01835 for(EEDigiCollection::const_iterator it = eeDigis.begin();
01836 it != eeDigis.end(); ++it){
01837 const EEDataFrame& frame = *it;
01838 const EcalTrigTowerDetId& ttId = triggerTowerMap_->towerOf(frame.id());
01839 const int iTtEta0 = iTtEta2cIndex(ttId.ieta());
01840 const int iTtPhi0 = iTtPhi2cIndex(ttId.iphi());
01841
01842
01843
01844
01845 double theta = eeGeometry->getGeometry(frame.id())->getPosition().theta();
01846 double e = frame2EnergyForTp(frame);
01847 if((frame2EnergyForTp(frame,-1) < e) && (frame2EnergyForTp(frame, 1) < e)){
01848 ttEtSums[iTtEta0][iTtPhi0] += e*sin(theta);
01849 }
01850 }
01851
01852
01853 int innerTTEtas[] = {0, 1, 54, 55};
01854 for(unsigned iRing = 0; iRing < sizeof(innerTTEtas)/sizeof(innerTTEtas[0]);
01855 ++iRing){
01856 int iTtEta0 = innerTTEtas[iRing];
01857
01858
01859
01860
01861
01862
01863
01864
01865
01866
01867
01868
01869
01870
01871
01872 for(unsigned iTtPhi0 = 0; iTtPhi0 < nTtPhi-1; iTtPhi0 += 2){
01873 double et = .5*(ttEtSums[iTtEta0][iTtPhi0]
01874 +ttEtSums[iTtEta0][iTtPhi0+1]);
01875
01876
01877
01878 ttEtSums[iTtEta0][iTtPhi0] = et;
01879 ttEtSums[iTtEta0][iTtPhi0+1] = et;
01880 }
01881 }
01882 }
01883
01884 template<class T>
01885 double EcalSelectiveReadoutValidation::frame2EnergyForTp(const T& frame,
01886 int offset) const{
01887
01888
01889 double weights[] = {0., -1/3., -1/3., -1/3., 0., 1.};
01890
01891 double adc2GeV = 0.;
01892 if(typeid(frame) == typeid(EBDataFrame)){
01893 adc2GeV = 0.035;
01894 } else if(typeid(frame) == typeid(EEDataFrame)){
01895 adc2GeV = 0.060;
01896 } else{
01897
01898 throw cms::Exception("Severe Error")
01899 << __FILE__ << ":" << __LINE__ << ": "
01900 << "this is a bug. Please report it.\n";
01901 }
01902
01903 double acc = 0;
01904
01905 const int n = min<int>(frame.size(), sizeof(weights)/sizeof(weights[0]));
01906
01907 double gainInv[] = {12., 1., 6., 12};
01908
01909 for(int i=offset; i < n; ++i){
01910 int iframe = i + offset;
01911 if(iframe>=0 && iframe<frame.size()){
01912 acc += weights[i]*frame[iframe].adc()
01913 *gainInv[frame[iframe].gainId()]*adc2GeV;
01914
01915
01916
01917 }
01918 }
01919
01920 return acc;
01921 }
01922
01923 MonitorElement* EcalSelectiveReadoutValidation::bookFloat(const std::string& name){
01924 if(!registerHist(name, "")) return 0;
01925 MonitorElement* result = dbe_->bookFloat(name);
01926 if(result==0){
01927 throw cms::Exception("DQM")
01928 << "Failed to book integer DQM monitor element" << name;
01929 }
01930 return result;
01931 }
01932
01933
01934 MonitorElement* EcalSelectiveReadoutValidation::book1D(const std::string& name, const std::string& title, int nbins, double xmin, double xmax){
01935 if(!registerHist(name, title)) return 0;
01936 MonitorElement* result = dbe_->book1D(name, title, nbins, xmin, xmax);
01937 if(result==0){
01938 throw cms::Exception("Histo")
01939 << "Failed to book histogram " << name;
01940 }
01941 return result;
01942 }
01943
01944 MonitorElement* EcalSelectiveReadoutValidation::book2D(const std::string& name, const std::string& title, int nxbins, double xmin, double xmax, int nybins, double ymin, double ymax){
01945 if(!registerHist(name, title)) return 0;
01946 MonitorElement* result = dbe_->book2D(name, title, nxbins, xmin, xmax,
01947 nybins, ymin, ymax);
01948 if(result==0){
01949 throw cms::Exception("Histo")
01950 << "Failed to book histogram " << name;
01951 }
01952 return result;
01953 }
01954
01955 MonitorElement* EcalSelectiveReadoutValidation::bookProfile(const std::string& name, const std::string& title, int nbins, double xmin, double xmax){
01956 if(!registerHist(name, title)) return 0;
01957 MonitorElement* result = dbe_->bookProfile(name, title, nbins, xmin, xmax,
01958 0, 0, 0);
01959 if(result==0){
01960 throw cms::Exception("Histo")
01961 << "Failed to book histogram " << name;
01962 }
01963 return result;
01964 }
01965
01966 MonitorElement* EcalSelectiveReadoutValidation::bookProfile2D(const std::string& name, const std::string& title, int nbinx, double xmin, double xmax, int nbiny, double ymin, double ymax, const char* option){
01967 if(!registerHist(name, title)) return 0;
01968 MonitorElement* result
01969 = dbe_->bookProfile2D(name,
01970 title,
01971 nbinx, xmin, xmax,
01972 nbiny, ymin, ymax,
01973 0, 0, 0,
01974 option);
01975 if(result==0){
01976 throw cms::Exception("Histo")
01977 << "Failed to book histogram " << name;
01978 }
01979 return result;
01980 }
01981
01982 bool EcalSelectiveReadoutValidation::registerHist(const std::string& name,
01983 const std::string& title){
01984 availableHistList_.insert(pair<string, string>(name, title));
01985 return allHists_ || histList_.find(name)!=histList_.end();
01986 }
01987
01988 void EcalSelectiveReadoutValidation::readAllCollections(const edm::Event& event){
01989 ebRecHits_.read(event);
01990 eeRecHits_.read(event);
01991 ebDigis_.read(event);
01992 eeDigis_.read(event);
01993 ebNoZsDigis_.read(event);
01994 eeNoZsDigis_.read(event);
01995 ebSrFlags_.read(event);
01996 eeSrFlags_.read(event);
01997 ebComputedSrFlags_.read(event);
01998 eeComputedSrFlags_.read(event);
01999 ebSimHits_.read(event);
02000 eeSimHits_.read(event);
02001 tps_.read(event);
02002 fedRaw_.read(event);
02003 }
02004
02005 void EcalSelectiveReadoutValidation::printAvailableHists(){
02006 LogInfo log("HistoList");
02007 log << "Avalailable histograms (DQM monitor elements): \n";
02008 for(map<string, string>::iterator it = availableHistList_.begin();
02009 it != availableHistList_.end();
02010 ++it){
02011 log << it->first << ": " << it->second << "\n";
02012 }
02013 log << "\nTo include an histogram add its name in the vstring parameter "
02014 "'histograms' of the EcalSelectiveReadoutValidation module\n";
02015 }
02016
02017 double EcalSelectiveReadoutValidation::getEbEventSize(double nReadXtals) const{
02018 double ruHeaderPayload = 0.;
02019 const int firstEbDcc0 = nEeDccs/2;
02020 for(int iDcc0 = firstEbDcc0; iDcc0 < firstEbDcc0 + nEbDccs; ++iDcc0){
02021 ruHeaderPayload += getRuCount(iDcc0)*8.;
02022 }
02023
02024 return getDccOverhead(EB)*nEbDccs + nReadXtals*getBytesPerCrystal()
02025 + ruHeaderPayload;
02026 }
02027
02028 double EcalSelectiveReadoutValidation::getEeEventSize(double nReadXtals) const{
02029 double ruHeaderPayload = 0.;
02030 const unsigned firstEbDcc0 = nEeDccs/2;
02031 for(unsigned iDcc0 = 0; iDcc0 < nDccs_; ++iDcc0){
02032
02033 if(iDcc0== firstEbDcc0) iDcc0 += nEbDccs;
02034 ruHeaderPayload += getRuCount(iDcc0)*8.;
02035 }
02036 return getDccOverhead(EE)*nEeDccs + nReadXtals*getBytesPerCrystal()
02037 + ruHeaderPayload;
02038 }
02039
02040 void EcalSelectiveReadoutValidation::normalizeHists(double eventCount){
02041 MonitorElement* mes[] = { meChOcc_, meTtf_, meTp_, meFullRoRu_,
02042 meZs1Ru_, meForcedRu_, meLiTtf_, meMiTtf_, meHiTtf_,
02043
02044
02045 };
02046
02047 double scale = 1./eventCount;
02048 stringstream buf;
02049 for(unsigned i = 0; i < sizeof(mes)/sizeof(mes[0]); ++i){
02050 if(mes[i] == 0) continue;
02051 TH1* h = mes[i]->getTH1();
02052 if(dynamic_cast<TH2*>(h)){
02053 h->GetZaxis()->SetTitle("Frequency");
02054 } else{
02055 h->GetYaxis()->SetTitle("<Count>");
02056 }
02057 buf << "Normalising " << h->GetName() << ". Factor: " << scale << "\n";
02058 h->Scale(scale);
02059
02060
02061 h->SetBit(TH1::kIsAverage);
02062 }
02063 edm::LogInfo("EcalSrValid") << buf.str();
02064 }
02065
02066
02067
02068 int
02069 EcalSelectiveReadoutValidation::dccZsFIR(const EcalDataFrame& frame,
02070 const std::vector<int>& firWeights,
02071 int firstFIRSample,
02072 bool* saturated){
02073 const int nFIRTaps = 6;
02074
02075 const vector<int>& w = firWeights;
02076
02077
02078 int acc = 0;
02079 bool gain12saturated = false;
02080 const int gain12 = 0x01;
02081 const int lastFIRSample = firstFIRSample + nFIRTaps - 1;
02082
02083 int iWeight = 0;
02084 for(int iSample=firstFIRSample-1;
02085 iSample<lastFIRSample; ++iSample, ++iWeight){
02086 if(iSample>=0 && iSample < frame.size()){
02087 EcalMGPASample sample(frame[iSample]);
02088 if(sample.gainId()!=gain12) gain12saturated = true;
02089 LogTrace("DccFir") << (iSample>=firstFIRSample?"+":"") << sample.adc()
02090 << "*(" << w[iWeight] << ")";
02091 acc+=sample.adc()*w[iWeight];
02092 } else{
02093 edm::LogWarning("DccFir") << __FILE__ << ":" << __LINE__ <<
02094 ": Not enough samples in data frame or 'ecalDccZs1stSample' module "
02095 "parameter is not valid...";
02096 }
02097 }
02098 LogTrace("DccFir") << "\n";
02099
02100
02101
02102 acc = (acc>=0)?(acc >> 8):-(-acc >> 8);
02103
02104
02105
02106
02107 LogTrace("DccFir") << "acc: " << acc << "\n"
02108 << "saturated: " << (gain12saturated?"yes":"no") << "\n";
02109
02110 if(saturated){
02111 *saturated = gain12saturated;
02112 }
02113
02114 return gain12saturated?numeric_limits<int>::max():acc;
02115 }
02116
02117 std::vector<int>
02118 EcalSelectiveReadoutValidation::getFIRWeights(const std::vector<double>&
02119 normalizedWeights){
02120 const int nFIRTaps = 6;
02121 vector<int> firWeights(nFIRTaps, 0);
02122 const static int maxWeight = 0xEFF;
02123 for(unsigned i=0; i < min((size_t)nFIRTaps,normalizedWeights.size()); ++i){
02124 firWeights[i] = lround(normalizedWeights[i] * (1<<10));
02125 if(abs(firWeights[i])>maxWeight){
02126 firWeights[i] = firWeights[i]<0?-maxWeight:maxWeight;
02127 }
02128 }
02129 return firWeights;
02130 }
02131
02132 void
02133 EcalSelectiveReadoutValidation::configFirWeights(vector<double> weightsForZsFIR){
02134 bool notNormalized = false;
02135 bool notInt = false;
02136 for(unsigned i=0; i < weightsForZsFIR.size(); ++i){
02137 if(weightsForZsFIR[i] > 1.) notNormalized = true;
02138 if((int)weightsForZsFIR[i]!=weightsForZsFIR[i]) notInt = true;
02139 }
02140 if(notInt && notNormalized){
02141 throw cms::Exception("InvalidParameter")
02142 << "weigtsForZsFIR paramater values are not valid: they "
02143 << "must either be integer and uses the hardware representation "
02144 << "of the weights or less or equal than 1 and used the normalized "
02145 << "representation.";
02146 }
02147 LogInfo log("DccFir");
02148 if(notNormalized){
02149 firWeights_ = vector<int>(weightsForZsFIR.size());
02150 for(unsigned i = 0; i< weightsForZsFIR.size(); ++i){
02151 firWeights_[i] = (int)weightsForZsFIR[i];
02152 }
02153 } else{
02154 firWeights_ = getFIRWeights(weightsForZsFIR);
02155 }
02156
02157 log << "Input weights for FIR: ";
02158 for(unsigned i = 0; i < weightsForZsFIR.size(); ++i){
02159 log << weightsForZsFIR[i] << "\t";
02160 }
02161
02162 double s2 = 0.;
02163 log << "\nActual FIR weights: ";
02164 for(unsigned i = 0; i < firWeights_.size(); ++i){
02165 log << firWeights_[i] << "\t";
02166 s2 += firWeights_[i]*firWeights_[i];
02167 }
02168
02169 s2 = sqrt(s2);
02170 log << "\nNormalized FIR weights after hw representation rounding: ";
02171 for(unsigned i = 0; i < firWeights_.size(); ++i){
02172 log << firWeights_[i] / (double)(1<<10) << "\t";
02173 }
02174
02175 log <<"\nFirst FIR sample: " << firstFIRSample_;
02176 }
02177
02178 void EcalSelectiveReadoutValidation::initAsciiFile(){
02179 if(logSrpAlgoErrors_){
02180 srpAlgoErrorLog_.open(srpAlgoErrorLogFileName_.c_str(), ios::out | ios::trunc);
02181 if(!srpAlgoErrorLog_.good()){
02182 throw cms::Exception("Output")
02183 << "Failed to open the log file '"
02184 << srpAlgoErrorLogFileName_
02185 << "' for SRP algorithm result check.\n";
02186 }
02187 }
02188
02189 if(logSrApplicationErrors_){
02190 srApplicationErrorLog_.open(srApplicationErrorLogFileName_.c_str(), ios::out | ios::trunc);
02191 if(!srApplicationErrorLog_.good()){
02192 throw cms::Exception("Output")
02193 << "Failed to open the log file '"
02194 << srApplicationErrorLogFileName_
02195 << "' for Selective Readout decision application check.\n";
02196 }
02197 }
02198 }
02199
02200
02201
02202
02203 template<class T>
02204 void EcalSelectiveReadoutValidation::compareSrfColl(const edm::Event& event, T& srfFromData, T& computedSrf){
02205 typedef typename T::const_iterator SrFlagCollectionConstIt;
02206 typedef typename T::key_type MyRuDetIdType;
02207 SrFlagCollectionConstIt itSrfFromData = srfFromData.begin();
02208 SrFlagCollectionConstIt itComputedSr = computedSrf.begin();
02209
02210 {
02211 PgTiming t("collection comparison");
02212
02213
02214
02215 while(itSrfFromData != srfFromData.end()
02216 || itComputedSr != computedSrf.end()){
02217
02218 MyRuDetIdType inconsistentRu = 0;
02219 bool inconsistent = false;
02220 if(itComputedSr == computedSrf.end() ||
02221 (itSrfFromData != srfFromData.end()
02222 && itSrfFromData->id() < itComputedSr->id())){
02223
02224 pair<int, int> ch = dccCh(itSrfFromData->id());
02225 srpAlgoErrorLog_ << event.id() << ": " << itSrfFromData->id()
02226 << ", DCC " << ch.first << " ch " << ch.second
02227 << " found in data (SRF:" << itSrfFromData->flagName()
02228 << ") but not in the set of SRFs computed from the data TTF.\n";
02229 inconsistentRu = itSrfFromData->id();
02230 inconsistent = true;
02231 ++itSrfFromData;
02232 } else if(itSrfFromData==srfFromData.end() ||
02233 (itComputedSr != computedSrf.end()
02234 && itComputedSr->id() < itSrfFromData->id())){
02235
02236 pair<int, int> ch = dccCh(itComputedSr->id());
02237 if(logErrForDccs_[ch.first-minDccId_]){
02238 srpAlgoErrorLog_ << event.id() << ": " << itComputedSr->id()
02239 << ", DCC " << ch.first << " ch " << ch.second
02240 << " not found in data. Computed SRF: "
02241 << itComputedSr->flagName() << ".\n";
02242 inconsistentRu = itComputedSr->id();
02243 inconsistent = true;
02244 }
02245 ++itComputedSr;
02246 } else{
02247
02248 if(itComputedSr->value()!=itSrfFromData->value()){
02249
02250 pair<int, int> ch = dccCh(itSrfFromData->id());
02251 srpAlgoErrorLog_ << event.id() << ", "
02252 << itSrfFromData->id()
02253 << ", DCC " << ch.first << " ch " << ch.second
02254 << ", SRF inconsistency: "
02255 << "from data: " << itSrfFromData->flagName()
02256 << ", computed from TTF: "
02257 << itComputedSr->flagName()
02258 << "\n";
02259
02260 inconsistentRu = itComputedSr->id();
02261 inconsistent = true;
02262 }
02263 if(itComputedSr != computedSrf.end()) ++itComputedSr;
02264 if(itSrfFromData != srfFromData.end()) ++itSrfFromData;
02265 }
02266
02267 if(inconsistent) fill(meSRFlagsConsistency_, ruGraphX(inconsistentRu),
02268 ruGraphY(inconsistentRu));
02269 }
02270 }
02271 }
02272
02273
02274 int EcalSelectiveReadoutValidation::dccId(const EcalScDetId& detId) const{
02275 return elecMap_->getDCCandSC(detId).first;
02276 }
02277
02278 int EcalSelectiveReadoutValidation::dccId(const EcalTrigTowerDetId& detId) const{
02279 if(detId.ietaAbs()>17){
02280 throw cms::Exception("InvalidArgument")
02281 << "Argument of EcalSelectiveReadoutValidation::dccId(const EcalTrigTowerDetId&) "
02282 << "must be a barrel trigger tower Id\n";
02283 }
02284 return dccCh(detId).first;
02285
02286
02287
02288
02289
02290 }
02291
02292
02293
02294 void EcalSelectiveReadoutValidation::selectFedsForLog(){
02295 logErrForDccs_ = vector<bool>(nDccs_, false);
02296
02297 for(EBSrFlagCollection::const_iterator it = ebSrFlags_->begin();
02298 it != ebSrFlags_->end();
02299 ++it){
02300
02301
02302
02303
02304 int iDcc = dccId(it->id()) - minDccId_;
02305
02306
02307 logErrForDccs_.at(iDcc) = true;
02308 }
02309
02310 for(EESrFlagCollection::const_iterator it = eeSrFlags_->begin();
02311 it != eeSrFlags_->end();
02312 ++it){
02313 int iDcc = dccId(it->id()) - minDccId_;
02314
02315
02316 logErrForDccs_.at(iDcc) = true;
02317 }
02318
02319 stringstream buf;
02320 buf << "List of DCCs found in the first processed event: ";
02321 bool first = true;
02322 for(unsigned iDcc = 0; iDcc < nDccs_; ++iDcc){
02323 if(logErrForDccs_[iDcc]){
02324 buf << (first?"":", ") << (iDcc + minDccId_);
02325 first = false;
02326 }
02327 }
02328 buf << "\nOnly DCCs from this list will be considered for error logging\n";
02329 srpAlgoErrorLog_ << buf.str();
02330 srApplicationErrorLog_<< buf.str();
02331 LogInfo("EcalSrValid") << buf;
02332 }
02333
02334
02335 template<class T>
02336 void EcalSelectiveReadoutValidation::checkSrApplication(const edm::Event& event,
02337 T& srfs){
02338 typedef typename T::const_iterator SrFlagCollectionConstIt;
02339 typedef typename T::key_type MyRuDetIdType;
02340
02341 for(SrFlagCollectionConstIt itSrf = srfs.begin();
02342 itSrf != srfs.end(); ++itSrf){
02343 int flag = itSrf->value() & ~EcalSrFlag::SRF_FORCED_MASK;
02344 pair<int,int> ru = dccCh(itSrf->id());
02345
02346 if(flag == EcalSrFlag::SRF_FULL){
02347 if(nPerRu_[ru.first-minDccId_][ru.second-1]==getCrystalCount(ru.first, ru.second)){
02348 fill(meIncompleteFRORateMap_, ruGraphX(itSrf->id()),
02349 ruGraphY(itSrf->id()), 0);
02350 fill(meDroppedFRORateMap_,
02351 ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 0);
02352 } else if(nPerRu_[ru.first-minDccId_][ru.second-1]==0) {
02353 fill(meIncompleteFRORateMap_,
02354 ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 0);
02355 fill(meDroppedFRORateMap_,
02356 ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 1);
02357 fill(meDroppedFROMap_, ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 1);
02358 ++nDroppedFRO_;
02359 srApplicationErrorLog_ << event.id() << ": Flag of RU "
02360 << itSrf->id() << " (DCC " << ru.first
02361 << " ch " << ru.second << ") is 'Full readout' "
02362 << "while none of its channel was read out\n";
02363 } else{
02364 fill(meIncompleteFRORateMap_,
02365 ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 1);
02366 fill(meDroppedFRORateMap_,
02367 ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 0);
02368 fill(meIncompleteFROMap_,
02369 ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 1);
02370 ++nIncompleteFRO_;
02371 srApplicationErrorLog_ << event.id() << ": Flag of RU"
02372 << itSrf->id() << " (DCC " << ru.first
02373 << " ch " << ru.second << ") is 'Full readout' "
02374 << "while only "
02375 << nPerRu_[ru.first-minDccId_][ru.second-1]
02376 << " / " << getCrystalCount(ru.first, ru.second)
02377 << " channels were read out.\n";
02378 }
02379 }
02380
02381 if(flag == EcalSrFlag::SRF_ZS1 || flag == EcalSrFlag::SRF_ZS2){
02382 if(nPerRu_[ru.first-minDccId_][ru.second-1]
02383 ==getCrystalCount(ru.first, ru.second)){
02384
02385
02386 fill(meCompleteZSMap_, ruGraphX(itSrf->id()), ruGraphY(itSrf->id()));
02387 fill(meCompleteZSRateMap_,
02388 ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 1);
02389
02390
02391
02392 ++nCompleteZS_;
02393 } else{
02394 fill(meCompleteZSRateMap_,
02395 ruGraphX(itSrf->id()), ruGraphY(itSrf->id()), 0);
02396 }
02397 }
02398 }
02399 }
02400
02401 int EcalSelectiveReadoutValidation::getCrystalCount(int iDcc, int iDccCh){
02402 if(iDcc < minDccId_ || iDcc > maxDccId_){
02403 return 0;
02404 } else if (10 <= iDcc && iDcc <= 45) {
02405 return 25;
02406 } else {
02407 int iDccPhi;
02408 if(iDcc < 10) iDccPhi = iDcc;
02409 else iDccPhi = iDcc - 45;
02410 switch(iDccPhi*100+iDccCh){
02411 case 110:
02412 case 232:
02413 case 312:
02414 case 412:
02415 case 532:
02416 case 610:
02417 case 830:
02418 case 806:
02419
02420 return 20;
02421 case 134:
02422 case 634:
02423 case 827:
02424 case 803:
02425 return 10;
02426 case 330:
02427 case 430:
02428 return 20;
02429 case 203:
02430 case 503:
02431 case 721:
02432 case 921:
02433 return 21;
02434 default:
02435 return 25;
02436 }
02437 }
02438 }
1.7.3