Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-01-13 23:40:10 

0001 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0002 #include "Geometry/HcalCommonData/interface/HcalHitRelabeller.h"
0003 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0004 #include "Validation/HcalRecHits/interface/HcalRecHitsValidation.h"
0005 
0006 //#define EDM_ML_DEBUG
0007 
0008 HcalRecHitsValidation::HcalRecHitsValidation(edm::ParameterSet const &conf)
0009     : topFolderName_(conf.getParameter<std::string>("TopFolderName")),
0010       outputFile_(conf.getUntrackedParameter<std::string>("outputFile", "myfile.root")),
0011       hcalselector_(conf.getUntrackedParameter<std::string>("hcalselector", "all")),
0012       ecalselector_(conf.getUntrackedParameter<std::string>("ecalselector", "yes")),
0013       sign_(conf.getUntrackedParameter<std::string>("sign", "*")),
0014       mc_(conf.getUntrackedParameter<std::string>("mc", "yes")),
0015       testNumber_(conf.getParameter<bool>("TestNumber")),
0016       EBRecHitCollectionLabel_(conf.getParameter<edm::InputTag>("EBRecHitCollectionLabel")),
0017       EERecHitCollectionLabel_(conf.getParameter<edm::InputTag>("EERecHitCollectionLabel")),
0018       tok_evt_(consumes<edm::HepMCProduct>(edm::InputTag("generatorSmeared"))),
0019       tok_EB_(consumes<EBRecHitCollection>(EBRecHitCollectionLabel_)),
0020       tok_EE_(consumes<EERecHitCollection>(EERecHitCollectionLabel_)),
0021       tok_hh_(consumes<edm::PCaloHitContainer>(conf.getUntrackedParameter<edm::InputTag>("SimHitCollectionLabel"))),
0022       tok_hbhe_(consumes<HBHERecHitCollection>(conf.getUntrackedParameter<edm::InputTag>("HBHERecHitCollectionLabel"))),
0023       tok_hf_(consumes<HFRecHitCollection>(conf.getUntrackedParameter<edm::InputTag>("HFRecHitCollectionLabel"))),
0024       tok_ho_(consumes<HORecHitCollection>(conf.getUntrackedParameter<edm::InputTag>("HORecHitCollectionLabel"))),
0025       tok_HRNDC_(esConsumes<HcalDDDRecConstants, HcalRecNumberingRecord>()),
0026       tok_Geom_(esConsumes<CaloGeometry, CaloGeometryRecord>()) {
0027   // DQM ROOT output
0028   if (!outputFile_.empty()) {
0029     edm::LogVerbatim("OutputInfo") << " Hcal RecHit Task histograms will be saved to '" << outputFile_.c_str() << "'";
0030   } else {
0031     edm::LogVerbatim("OutputInfo") << " Hcal RecHit Task histograms will NOT be saved";
0032   }
0033 
0034   nevtot = 0;
0035 
0036   // Collections
0037 
0038   // register for data access
0039 
0040   subdet_ = 5;
0041   if (hcalselector_ == "noise")
0042     subdet_ = 0;
0043   if (hcalselector_ == "HB")
0044     subdet_ = 1;
0045   if (hcalselector_ == "HE")
0046     subdet_ = 2;
0047   if (hcalselector_ == "HO")
0048     subdet_ = 3;
0049   if (hcalselector_ == "HF")
0050     subdet_ = 4;
0051   if (hcalselector_ == "all")
0052     subdet_ = 5;
0053   if (hcalselector_ == "ZS")
0054     subdet_ = 6;
0055 
0056   iz = 1;
0057   if (sign_ == "-")
0058     iz = -1;
0059   if (sign_ == "*")
0060     iz = 0;
0061 
0062   imc = 1;
0063   if (mc_ == "no")
0064     imc = 0;
0065 }
0066 
0067 void HcalRecHitsValidation::bookHistograms(DQMStore::IBooker &ib, edm::Run const &run, edm::EventSetup const &es) {
0068   Char_t histo[200];
0069 
0070   ib.setCurrentFolder(topFolderName_);
0071 
0072   //======================= Now various cases one by one ===================
0073 
0074   // Histograms drawn for single pion scan
0075   if (subdet_ != 0 && imc != 0) {  // just not for noise
0076     sprintf(histo, "HcalRecHitTask_En_rechits_cone_profile_vs_ieta_all_depths");
0077     meEnConeEtaProfile = ib.bookProfile(histo, histo, 83, -41.5, 41.5, -100., 2000., " ");
0078 
0079     sprintf(histo, "HcalRecHitTask_En_rechits_cone_profile_vs_ieta_all_depths_E");
0080     meEnConeEtaProfile_E = ib.bookProfile(histo, histo, 83, -41.5, 41.5, -100., 2000., " ");
0081 
0082     sprintf(histo, "HcalRecHitTask_En_rechits_cone_profile_vs_ieta_all_depths_EH");
0083     meEnConeEtaProfile_EH = ib.bookProfile(histo, histo, 83, -41.5, 41.5, -100., 2000., " ");
0084   }
0085 
0086   // ************** HB **********************************
0087   if (subdet_ == 1 || subdet_ == 5) {
0088     sprintf(histo, "HcalRecHitTask_M2Log10Chi2_of_rechits_HB");  //   Chi2
0089     meRecHitsM2Chi2HB = ib.book1D(histo, histo, 120, -2., 10.);
0090 
0091     sprintf(histo, "HcalRecHitTask_Log10Chi2_vs_energy_profile_HB");
0092     meLog10Chi2profileHB = ib.bookProfile(histo, histo, 300, -5., 295., -2., 9.9, " ");
0093 
0094     sprintf(histo, "HcalRecHitTask_energy_of_rechits_HB");
0095     meRecHitsEnergyHB = ib.book1D(histo, histo, 2010, -10., 2000.);
0096 
0097     sprintf(histo, "HcalRecHitTask_timing_vs_energy_profile_HB");
0098     meTEprofileHB = ib.bookProfile(histo, histo, 150, -5., 295., -48., 92., " ");
0099 
0100     sprintf(histo, "HcalRecHitTask_timing_vs_energy_profile_Low_HB");
0101     meTEprofileHB_Low = ib.bookProfile(histo, histo, 150, -5., 295., -48., 92., " ");
0102 
0103     sprintf(histo, "HcalRecHitTask_timing_vs_energy_profile_High_HB");
0104     meTEprofileHB_High = ib.bookProfile(histo, histo, 150, -5., 295., 48., 92., " ");
0105 
0106     if (imc != 0) {
0107       sprintf(histo, "HcalRecHitTask_energy_rechits_vs_simhits_HB");
0108       meRecHitSimHitHB = ib.book2D(histo, histo, 120, 0., 1.2, 300, 0., 150.);
0109       sprintf(histo, "HcalRecHitTask_energy_rechits_vs_simhits_profile_HB");
0110       meRecHitSimHitProfileHB = ib.bookProfile(histo, histo, 120, 0., 1.2, 0., 500., " ");
0111     }
0112   }
0113 
0114   // ********************** HE ************************************
0115   if (subdet_ == 2 || subdet_ == 5) {
0116     sprintf(histo, "HcalRecHitTask_M2Log10Chi2_of_rechits_HE");  //   Chi2
0117     meRecHitsM2Chi2HE = ib.book1D(histo, histo, 120, -2., 10.);
0118 
0119     sprintf(histo, "HcalRecHitTask_Log10Chi2_vs_energy_profile_HE");
0120     meLog10Chi2profileHE = ib.bookProfile(histo, histo, 1000, -5., 995., -2., 9.9, " ");
0121 
0122     sprintf(histo, "HcalRecHitTask_energy_of_rechits_HE");
0123     meRecHitsEnergyHE = ib.book1D(histo, histo, 2010, -10., 2000.);
0124 
0125     sprintf(histo, "HcalRecHitTask_timing_vs_energy_profile_Low_HE");
0126     meTEprofileHE_Low = ib.bookProfile(histo, histo, 80, -5., 75., -48., 92., " ");
0127 
0128     sprintf(histo, "HcalRecHitTask_timing_vs_energy_profile_HE");
0129     meTEprofileHE = ib.bookProfile(histo, histo, 200, -5., 2995., -48., 92., " ");
0130 
0131     if (imc != 0) {
0132       sprintf(histo, "HcalRecHitTask_energy_rechits_vs_simhits_HE");
0133       meRecHitSimHitHE = ib.book2D(histo, histo, 120, 0., 0.6, 300, 0., 150.);
0134       sprintf(histo, "HcalRecHitTask_energy_rechits_vs_simhits_profile_HE");
0135       meRecHitSimHitProfileHE = ib.bookProfile(histo, histo, 120, 0., 0.6, 0., 500., " ");
0136     }
0137   }
0138 
0139   // ************** HO ****************************************
0140   if (subdet_ == 3 || subdet_ == 5) {
0141     sprintf(histo, "HcalRecHitTask_energy_of_rechits_HO");
0142     meRecHitsEnergyHO = ib.book1D(histo, histo, 2010, -10., 2000.);
0143 
0144     sprintf(histo, "HcalRecHitTask_timing_vs_energy_profile_HO");
0145     meTEprofileHO = ib.bookProfile(histo, histo, 60, -5., 55., -48., 92., " ");
0146 
0147     sprintf(histo, "HcalRecHitTask_timing_vs_energy_profile_High_HO");
0148     meTEprofileHO_High = ib.bookProfile(histo, histo, 100, -5., 995., -48., 92., " ");
0149 
0150     if (imc != 0) {
0151       sprintf(histo, "HcalRecHitTask_energy_rechits_vs_simhits_HO");
0152       meRecHitSimHitHO = ib.book2D(histo, histo, 150, 0., 1.5, 350, 0., 350.);
0153       sprintf(histo, "HcalRecHitTask_energy_rechits_vs_simhits_profile_HO");
0154       meRecHitSimHitProfileHO = ib.bookProfile(histo, histo, 150, 0., 1.5, 0., 500., " ");
0155     }
0156   }
0157 
0158   // ********************** HF ************************************
0159   if (subdet_ == 4 || subdet_ == 5) {
0160     sprintf(histo, "HcalRecHitTask_energy_of_rechits_HF");
0161     meRecHitsEnergyHF = ib.book1D(histo, histo, 2010, -10., 2000.);
0162 
0163     sprintf(histo, "HcalRecHitTask_timing_vs_energy_profile_Low_HF");
0164     meTEprofileHF_Low = ib.bookProfile(histo, histo, 100, -5., 195., -48., 92., " ");
0165 
0166     sprintf(histo, "HcalRecHitTask_timing_vs_energy_profile_HF");
0167     meTEprofileHF = ib.bookProfile(histo, histo, 200, -5., 995., -48., 92., " ");
0168 
0169     if (imc != 0) {
0170       sprintf(histo, "HcalRecHitTask_energy_rechits_vs_simhits_HF");
0171       meRecHitSimHitHF = ib.book2D(histo, histo, 50, 0., 50., 150, 0., 150.);
0172       sprintf(histo, "HcalRecHitTask_energy_rechits_vs_simhits_HFL");
0173       meRecHitSimHitHFL = ib.book2D(histo, histo, 50, 0., 50., 150, 0., 150.);
0174       sprintf(histo, "HcalRecHitTask_energy_rechits_vs_simhits_HFS");
0175       meRecHitSimHitHFS = ib.book2D(histo, histo, 50, 0., 50., 150, 0., 150.);
0176       sprintf(histo, "HcalRecHitTask_energy_rechits_vs_simhits_profile_HF");
0177       meRecHitSimHitProfileHF = ib.bookProfile(histo, histo, 50, 0., 50., 0., 500., " ");
0178       sprintf(histo, "HcalRecHitTask_energy_rechits_vs_simhits_profile_HFL");
0179       meRecHitSimHitProfileHFL = ib.bookProfile(histo, histo, 50, 0., 50., 0., 500., " ");
0180       sprintf(histo, "HcalRecHitTask_energy_rechits_vs_simhits_profile_HFS");
0181       meRecHitSimHitProfileHFS = ib.bookProfile(histo, histo, 50, 0., 50., 0., 500., " ");
0182     }
0183   }
0184 }
0185 
0186 void HcalRecHitsValidation::analyze(edm::Event const &ev, edm::EventSetup const &c) {
0187   using namespace edm;
0188 
0189   const HcalDDDRecConstants *hcons = &c.getData(tok_HRNDC_);
0190 
0191   // cuts for each subdet_ector mimiking  "Scheme B"
0192   //  double cutHB = 0.9, cutHE = 1.4, cutHO = 1.1, cutHFL = 1.2, cutHFS = 1.8;
0193 
0194   // energy in HCAL
0195   double eHcalConeHB = 0.;
0196   double eHcalConeHE = 0.;
0197   double eHcalConeHO = 0.;
0198   double eHcalConeHF = 0.;
0199   double eHcalConeHFL = 0.;
0200   double eHcalConeHFS = 0.;
0201 
0202   // Total numbet of RecHits in HCAL, in the cone, above 1 GeV theshold
0203   int nrechits = 0;
0204   int nrechitsCone = 0;
0205   int nrechitsThresh = 0;
0206 
0207   // energy in ECAL
0208   double eEcalCone = 0.;
0209   int numrechitsEcal = 0;
0210 
0211   // MC info
0212   double phi_MC = -999999.;  // phi of initial particle from HepMC
0213   double eta_MC = -999999.;  // eta of initial particle from HepMC
0214 
0215   // HCAL energy around MC eta-phi at all depths;
0216   double partR = 0.3;
0217 
0218   if (imc != 0) {
0219     const edm::Handle<edm::HepMCProduct> &evtMC = ev.getHandle(tok_evt_);  // generator in late 310_preX
0220     if (!evtMC.isValid()) {
0221       edm::LogVerbatim("HcalRecHitsValidation") << "no HepMCProduct found";
0222     } else {
0223 #ifdef EDM_ML_DEBUG
0224       edm::LogVerbatim("HcalRecHitsValidation") << "*** source HepMCProduct found";
0225 #endif
0226     }
0227 
0228     // MC particle with highest pt is taken as a direction reference
0229     double maxPt = -99999.;
0230     int npart = 0;
0231     const HepMC::GenEvent *myGenEvent = evtMC->GetEvent();
0232     for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0233          ++p) {
0234       double phip = (*p)->momentum().phi();
0235       double etap = (*p)->momentum().eta();
0236       //    phi_MC = phip;
0237       //    eta_MC = etap;
0238       double pt = (*p)->momentum().perp();
0239       if (pt > maxPt) {
0240         npart++;
0241         maxPt = pt;
0242         phi_MC = phip;
0243         eta_MC = etap;
0244       }
0245     }
0246 #ifdef EDM_ML_DEBUG
0247     edm::LogVerbatim("HcalRecHitsValidation") << "*** Max pT = " << maxPt;
0248 #endif
0249   }
0250 
0251 #ifdef EDM_ML_DEBUG
0252   edm::LogVerbatim("HcalRecHitsValidation") << "*** 2";
0253 #endif
0254 
0255   geometry_ = &c.getData(tok_Geom_);
0256 
0257   // Fill working vectors of HCAL RecHits quantities (all of these are drawn)
0258   fillRecHitsTmp(subdet_, ev);
0259 
0260 #ifdef EDM_ML_DEBUG
0261   edm::LogVerbatim("HcalRecHitsValidation") << "*** 3";
0262 #endif
0263 
0264   //===========================================================================
0265   // IN ALL other CASES : ieta-iphi maps
0266   //===========================================================================
0267 
0268   // ECAL
0269   if (ecalselector_ == "yes" && (subdet_ == 1 || subdet_ == 2 || subdet_ == 5)) {
0270     const edm::Handle<EBRecHitCollection> &rhitEB = ev.getHandle(tok_EB_);
0271 
0272     EcalRecHitCollection::const_iterator RecHit;
0273     EcalRecHitCollection::const_iterator RecHitEnd;
0274 
0275     if (rhitEB.isValid()) {
0276       RecHit = rhitEB.product()->begin();
0277       RecHitEnd = rhitEB.product()->end();
0278 
0279       for (; RecHit != RecHitEnd; ++RecHit) {
0280         EBDetId EBid = EBDetId(RecHit->id());
0281 
0282         auto cellGeometry = geometry_->getSubdetectorGeometry(EBid)->getGeometry(EBid);
0283         double eta = cellGeometry->getPosition().eta();
0284         double phi = cellGeometry->getPosition().phi();
0285         double en = RecHit->energy();
0286 
0287         double r = dR(eta_MC, phi_MC, eta, phi);
0288         if (r < partR) {
0289           eEcalCone += en;
0290           numrechitsEcal++;
0291         }
0292       }
0293     }
0294 
0295     const edm::Handle<EERecHitCollection> &rhitEE = ev.getHandle(tok_EE_);
0296 
0297     if (rhitEE.isValid()) {
0298       RecHit = rhitEE.product()->begin();
0299       RecHitEnd = rhitEE.product()->end();
0300 
0301       for (; RecHit != RecHitEnd; ++RecHit) {
0302         EEDetId EEid = EEDetId(RecHit->id());
0303 
0304         auto cellGeometry = geometry_->getSubdetectorGeometry(EEid)->getGeometry(EEid);
0305         double eta = cellGeometry->getPosition().eta();
0306         double phi = cellGeometry->getPosition().phi();
0307         double en = RecHit->energy();
0308 
0309         double r = dR(eta_MC, phi_MC, eta, phi);
0310         if (r < partR) {
0311           eEcalCone += en;
0312           numrechitsEcal++;
0313         }
0314       }
0315     }
0316   }  // end of ECAL selection
0317 
0318 #ifdef EDM_ML_DEBUG
0319   edm::LogVerbatim("HcalRecHitsValidation") << "*** 4";
0320 #endif
0321 
0322   //===========================================================================
0323   // SUBSYSTEMS,
0324   //===========================================================================
0325 
0326   if ((subdet_ != 6) && (subdet_ != 0)) {
0327 #ifdef EDM_ML_DEBUG
0328     edm::LogVerbatim("HcalRecHitsValidation") << "*** 6";
0329 #endif
0330 
0331     double HcalCone = 0.;
0332 
0333     int ietaMax = 9999;
0334     double etaMax = 9999.;
0335 
0336     //   CYCLE over cells ====================================================
0337 
0338     for (unsigned int i = 0; i < cen.size(); i++) {
0339       int sub = csub[i];
0340       int depth = cdepth[i];
0341       double eta = ceta[i];
0342       double phi = cphi[i];
0343       double en = cen[i];
0344       double t = ctime[i];
0345       int ieta = cieta[i];
0346       double chi2 = cchi2[i];
0347 
0348       double chi2_log10 = 9.99;  // initial value above histos limits , keep it if chi2 <= 0.
0349       if (chi2 > 0.)
0350         chi2_log10 = log10(chi2);
0351 
0352       nrechits++;
0353       if (en > 1.)
0354         nrechitsThresh++;
0355 
0356       double r = dR(eta_MC, phi_MC, eta, phi);
0357       if (r < partR) {
0358         if (sub == 1)
0359           eHcalConeHB += en;
0360         if (sub == 2)
0361           eHcalConeHE += en;
0362         if (sub == 3)
0363           eHcalConeHO += en;
0364         if (sub == 4) {
0365           eHcalConeHF += en;
0366           if (depth == 1)
0367             eHcalConeHFL += en;
0368           else
0369             eHcalConeHFS += en;
0370         }
0371         nrechitsCone++;
0372 
0373         HcalCone += en;
0374         // alternative: ietamax -> closest to MC eta  !!!
0375         float eta_diff = fabs(eta_MC - eta);
0376         if (eta_diff < etaMax) {
0377           etaMax = eta_diff;
0378           ietaMax = ieta;
0379         }
0380       }
0381 
0382       // The energy and overall timing histos are drawn while
0383       // the ones split by depth are not
0384       if (sub == 1 && (subdet_ == 1 || subdet_ == 5)) {
0385         meRecHitsM2Chi2HB->Fill(chi2_log10);
0386         meLog10Chi2profileHB->Fill(en, chi2_log10);
0387 
0388         meRecHitsEnergyHB->Fill(en);
0389 
0390         meTEprofileHB_Low->Fill(en, t);
0391         meTEprofileHB->Fill(en, t);
0392         meTEprofileHB_High->Fill(en, t);
0393       }
0394       if (sub == 2 && (subdet_ == 2 || subdet_ == 5)) {
0395         meRecHitsM2Chi2HE->Fill(chi2_log10);
0396         meLog10Chi2profileHE->Fill(en, chi2_log10);
0397 
0398         meRecHitsEnergyHE->Fill(en);
0399 
0400         meTEprofileHE_Low->Fill(en, t);
0401         meTEprofileHE->Fill(en, t);
0402       }
0403       if (sub == 4 && (subdet_ == 4 || subdet_ == 5)) {
0404         meRecHitsEnergyHF->Fill(en);
0405 
0406         meTEprofileHF_Low->Fill(en, t);
0407         meTEprofileHF->Fill(en, t);
0408       }
0409       if (sub == 3 && (subdet_ == 3 || subdet_ == 5)) {
0410         meRecHitsEnergyHO->Fill(en);
0411 
0412         meTEprofileHO->Fill(en, t);
0413         meTEprofileHO_High->Fill(en, t);
0414       }
0415     }
0416 
0417     if (imc != 0) {
0418       meEnConeEtaProfile->Fill(double(ietaMax), HcalCone);  //
0419       meEnConeEtaProfile_E->Fill(double(ietaMax), eEcalCone);
0420       meEnConeEtaProfile_EH->Fill(double(ietaMax), HcalCone + eEcalCone);
0421     }
0422 
0423 #ifdef EDM_ML_DEBUG
0424     edm::LogVerbatim("HcalRecHitsValidation") << "*** 7";
0425 #endif
0426   }
0427 
0428   // SimHits vs. RecHits
0429   if (subdet_ > 0 && subdet_ < 6 && imc != 0) {  // not noise
0430 
0431     const edm::Handle<PCaloHitContainer> &hcalHits = ev.getHandle(tok_hh_);
0432     if (hcalHits.isValid()) {
0433       const PCaloHitContainer *SimHitResult = hcalHits.product();
0434 
0435       double enSimHitsHB = 0.;
0436       double enSimHitsHE = 0.;
0437       double enSimHitsHO = 0.;
0438       double enSimHitsHF = 0.;
0439       double enSimHitsHFL = 0.;
0440       double enSimHitsHFS = 0.;
0441       // sum of SimHits in the cone
0442 
0443       for (std::vector<PCaloHit>::const_iterator SimHits = SimHitResult->begin(); SimHits != SimHitResult->end();
0444            ++SimHits) {
0445         int sub, depth;
0446         HcalDetId cell;
0447 
0448         if (testNumber_)
0449           cell = HcalHitRelabeller::relabel(SimHits->id(), hcons);
0450         else
0451           cell = HcalDetId(SimHits->id());
0452 
0453         sub = cell.subdet();
0454         depth = cell.depth();
0455 
0456         if (sub != subdet_ && subdet_ != 5)
0457           continue;  // If we are not looking at all of the subdetectors and the
0458                      // simhit doesn't come from the specific subdetector of
0459                      // interest, then we won't do any thing with it
0460 
0461         const HcalGeometry *cellGeometry =
0462             dynamic_cast<const HcalGeometry *>(geometry_->getSubdetectorGeometry(DetId::Hcal, cell.subdet()));
0463         double etaS = cellGeometry->getPosition(cell).eta();
0464         double phiS = cellGeometry->getPosition(cell).phi();
0465         double en = SimHits->energy();
0466 
0467         double r = dR(eta_MC, phi_MC, etaS, phiS);
0468 
0469         if (r < partR) {  // just energy in the small cone
0470 
0471           if (sub == static_cast<int>(HcalBarrel))
0472             enSimHitsHB += en;
0473           if (sub == static_cast<int>(HcalEndcap))
0474             enSimHitsHE += en;
0475           if (sub == static_cast<int>(HcalOuter))
0476             enSimHitsHO += en;
0477           if (sub == static_cast<int>(HcalForward)) {
0478             enSimHitsHF += en;
0479             if (depth == 1)
0480               enSimHitsHFL += en;
0481             else
0482               enSimHitsHFS += en;
0483           }
0484         }
0485       }
0486 
0487       // Now some histos with SimHits
0488 
0489       if (subdet_ == 4 || subdet_ == 5) {
0490         meRecHitSimHitHF->Fill(enSimHitsHF, eHcalConeHF);
0491         meRecHitSimHitProfileHF->Fill(enSimHitsHF, eHcalConeHF);
0492 
0493         meRecHitSimHitHFL->Fill(enSimHitsHFL, eHcalConeHFL);
0494         meRecHitSimHitProfileHFL->Fill(enSimHitsHFL, eHcalConeHFL);
0495         meRecHitSimHitHFS->Fill(enSimHitsHFS, eHcalConeHFS);
0496         meRecHitSimHitProfileHFS->Fill(enSimHitsHFS, eHcalConeHFS);
0497       }
0498       if (subdet_ == 1 || subdet_ == 5) {
0499         meRecHitSimHitHB->Fill(enSimHitsHB, eHcalConeHB);
0500         meRecHitSimHitProfileHB->Fill(enSimHitsHB, eHcalConeHB);
0501       }
0502       if (subdet_ == 2 || subdet_ == 5) {
0503         meRecHitSimHitHE->Fill(enSimHitsHE, eHcalConeHE);
0504         meRecHitSimHitProfileHE->Fill(enSimHitsHE, eHcalConeHE);
0505       }
0506       if (subdet_ == 3 || subdet_ == 5) {
0507         meRecHitSimHitHO->Fill(enSimHitsHO, eHcalConeHO);
0508         meRecHitSimHitProfileHO->Fill(enSimHitsHO, eHcalConeHO);
0509       }
0510     }
0511   }
0512 
0513   nevtot++;
0514 }
0515 
0516 ///////////////////////////////////////////////////////////////////////////////
0517 void HcalRecHitsValidation::fillRecHitsTmp(int subdet_, edm::Event const &ev) {
0518   using namespace edm;
0519 
0520   // initialize data vectors
0521   csub.clear();
0522   cen.clear();
0523   ceta.clear();
0524   cphi.clear();
0525   ctime.clear();
0526   cieta.clear();
0527   ciphi.clear();
0528   cdepth.clear();
0529   cz.clear();
0530   cchi2.clear();
0531 
0532   if (subdet_ == 1 || subdet_ == 2 || subdet_ == 5 || subdet_ == 6 || subdet_ == 0) {
0533     // HBHE
0534     const edm::Handle<HBHERecHitCollection> &hbhecoll = ev.getHandle(tok_hbhe_);
0535     if (hbhecoll.isValid()) {
0536       for (HBHERecHitCollection::const_iterator j = hbhecoll->begin(); j != hbhecoll->end(); j++) {
0537         HcalDetId cell(j->id());
0538         const HcalGeometry *cellGeometry =
0539             dynamic_cast<const HcalGeometry *>(geometry_->getSubdetectorGeometry(DetId::Hcal, cell.subdet()));
0540         double eta = cellGeometry->getPosition(cell).eta();
0541         double phi = cellGeometry->getPosition(cell).phi();
0542         double zc = cellGeometry->getPosition(cell).z();
0543         int sub = cell.subdet();
0544         int depth = cell.depth();
0545         int inteta = cell.ieta();
0546         int intphi = cell.iphi() - 1;
0547         double en = j->energy();
0548         double t = j->time();
0549         double chi2 = j->chi2();
0550 
0551         if ((iz > 0 && eta > 0.) || (iz < 0 && eta < 0.) || iz == 0) {
0552           csub.push_back(sub);
0553           cen.push_back(en);
0554           ceta.push_back(eta);
0555           cphi.push_back(phi);
0556           ctime.push_back(t);
0557           cieta.push_back(inteta);
0558           ciphi.push_back(intphi);
0559           cdepth.push_back(depth);
0560           cz.push_back(zc);
0561           cchi2.push_back(chi2);
0562         }
0563       }
0564     }
0565   }
0566 
0567   if (subdet_ == 4 || subdet_ == 5 || subdet_ == 6 || subdet_ == 0) {
0568     // HF
0569     const edm::Handle<HFRecHitCollection> &hfcoll = ev.getHandle(tok_hf_);
0570     if (hfcoll.isValid()) {
0571       for (HFRecHitCollection::const_iterator j = hfcoll->begin(); j != hfcoll->end(); j++) {
0572         HcalDetId cell(j->id());
0573         auto cellGeometry = geometry_->getSubdetectorGeometry(cell)->getGeometry(cell);
0574 
0575         double eta = cellGeometry->getPosition().eta();
0576         double phi = cellGeometry->getPosition().phi();
0577         double zc = cellGeometry->getPosition().z();
0578         int sub = cell.subdet();
0579         int depth = cell.depth();
0580         int inteta = cell.ieta();
0581         int intphi = cell.iphi() - 1;
0582         double en = j->energy();
0583         double t = j->time();
0584 
0585         if ((iz > 0 && eta > 0.) || (iz < 0 && eta < 0.) || iz == 0) {
0586           csub.push_back(sub);
0587           cen.push_back(en);
0588           ceta.push_back(eta);
0589           cphi.push_back(phi);
0590           ctime.push_back(t);
0591           cieta.push_back(inteta);
0592           ciphi.push_back(intphi);
0593           cdepth.push_back(depth);
0594           cz.push_back(zc);
0595           cchi2.push_back(0.);
0596         }
0597       }
0598     }
0599   }
0600 
0601   // HO
0602   if (subdet_ == 3 || subdet_ == 5 || subdet_ == 6 || subdet_ == 0) {
0603     const edm::Handle<HORecHitCollection> &hocoll = ev.getHandle(tok_ho_);
0604     if (hocoll.isValid()) {
0605       for (HORecHitCollection::const_iterator j = hocoll->begin(); j != hocoll->end(); j++) {
0606         HcalDetId cell(j->id());
0607         auto cellGeometry = geometry_->getSubdetectorGeometry(cell)->getGeometry(cell);
0608 
0609         double eta = cellGeometry->getPosition().eta();
0610         double phi = cellGeometry->getPosition().phi();
0611         double zc = cellGeometry->getPosition().z();
0612         int sub = cell.subdet();
0613         int depth = cell.depth();
0614         int inteta = cell.ieta();
0615         int intphi = cell.iphi() - 1;
0616         double t = j->time();
0617         double en = j->energy();
0618 
0619         if ((iz > 0 && eta > 0.) || (iz < 0 && eta < 0.) || iz == 0) {
0620           csub.push_back(sub);
0621           cen.push_back(en);
0622           ceta.push_back(eta);
0623           cphi.push_back(phi);
0624           ctime.push_back(t);
0625           cieta.push_back(inteta);
0626           ciphi.push_back(intphi);
0627           cdepth.push_back(depth);
0628           cz.push_back(zc);
0629           cchi2.push_back(0.);
0630         }
0631       }
0632     }
0633   }
0634 }
0635 
0636 double HcalRecHitsValidation::dR(double eta1, double phi1, double eta2, double phi2) {
0637   double PI = 3.1415926535898;
0638   double deltaphi = phi1 - phi2;
0639   if (phi2 > phi1) {
0640     deltaphi = phi2 - phi1;
0641   }
0642   if (deltaphi > PI) {
0643     deltaphi = 2. * PI - deltaphi;
0644   }
0645   double deltaeta = eta2 - eta1;
0646   double tmp = sqrt(deltaeta * deltaeta + deltaphi * deltaphi);
0647   return tmp;
0648 }
0649 
0650 double HcalRecHitsValidation::phi12(double phi1, double en1, double phi2, double en2) {
0651   // weighted mean value of phi1 and phi2
0652 
0653   double tmp;
0654   double PI = 3.1415926535898;
0655   double a1 = phi1;
0656   double a2 = phi2;
0657 
0658   if (a1 > 0.5 * PI && a2 < 0.)
0659     a2 += 2 * PI;
0660   if (a2 > 0.5 * PI && a1 < 0.)
0661     a1 += 2 * PI;
0662   tmp = (a1 * en1 + a2 * en2) / (en1 + en2);
0663   if (tmp > PI)
0664     tmp -= 2. * PI;
0665 
0666   return tmp;
0667 }
0668 
0669 double HcalRecHitsValidation::dPhiWsign(double phi1, double phi2) {
0670   // clockwise      phi2 w.r.t phi1 means "+" phi distance
0671   // anti-clockwise phi2 w.r.t phi1 means "-" phi distance
0672 
0673   double PI = 3.1415926535898;
0674   double a1 = phi1;
0675   double a2 = phi2;
0676   double tmp = a2 - a1;
0677   if (a1 * a2 < 0.) {
0678     if (a1 > 0.5 * PI)
0679       tmp += 2. * PI;
0680     if (a2 > 0.5 * PI)
0681       tmp -= 2. * PI;
0682   }
0683   return tmp;
0684 }
0685 
0686 DEFINE_FWK_MODULE(HcalRecHitsValidation);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!