Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:32:31 

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   // energy in ECAL
0203   double eEcalCone = 0.;
0204 
0205   // MC info
0206   double phi_MC = -999999.;  // phi of initial particle from HepMC
0207   double eta_MC = -999999.;  // eta of initial particle from HepMC
0208 
0209   // HCAL energy around MC eta-phi at all depths;
0210   double partR = 0.3;
0211 
0212   if (imc != 0) {
0213     const edm::Handle<edm::HepMCProduct> &evtMC = ev.getHandle(tok_evt_);  // generator in late 310_preX
0214     if (!evtMC.isValid()) {
0215       edm::LogVerbatim("HcalRecHitsValidation") << "no HepMCProduct found";
0216     } else {
0217 #ifdef EDM_ML_DEBUG
0218       edm::LogVerbatim("HcalRecHitsValidation") << "*** source HepMCProduct found";
0219 #endif
0220     }
0221 
0222     // MC particle with highest pt is taken as a direction reference
0223     double maxPt = -99999.;
0224     const HepMC::GenEvent *myGenEvent = evtMC->GetEvent();
0225     for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0226          ++p) {
0227       double phip = (*p)->momentum().phi();
0228       double etap = (*p)->momentum().eta();
0229       //    phi_MC = phip;
0230       //    eta_MC = etap;
0231       double pt = (*p)->momentum().perp();
0232       if (pt > maxPt) {
0233         maxPt = pt;
0234         phi_MC = phip;
0235         eta_MC = etap;
0236       }
0237     }
0238 #ifdef EDM_ML_DEBUG
0239     edm::LogVerbatim("HcalRecHitsValidation") << "*** Max pT = " << maxPt;
0240 #endif
0241   }
0242 
0243 #ifdef EDM_ML_DEBUG
0244   edm::LogVerbatim("HcalRecHitsValidation") << "*** 2";
0245 #endif
0246 
0247   geometry_ = &c.getData(tok_Geom_);
0248 
0249   // Fill working vectors of HCAL RecHits quantities (all of these are drawn)
0250   fillRecHitsTmp(subdet_, ev);
0251 
0252 #ifdef EDM_ML_DEBUG
0253   edm::LogVerbatim("HcalRecHitsValidation") << "*** 3";
0254 #endif
0255 
0256   //===========================================================================
0257   // IN ALL other CASES : ieta-iphi maps
0258   //===========================================================================
0259 
0260   // ECAL
0261   if (ecalselector_ == "yes" && (subdet_ == 1 || subdet_ == 2 || subdet_ == 5)) {
0262     const edm::Handle<EBRecHitCollection> &rhitEB = ev.getHandle(tok_EB_);
0263 
0264     EcalRecHitCollection::const_iterator RecHit;
0265     EcalRecHitCollection::const_iterator RecHitEnd;
0266 
0267     if (rhitEB.isValid()) {
0268       RecHit = rhitEB.product()->begin();
0269       RecHitEnd = rhitEB.product()->end();
0270 
0271       for (; RecHit != RecHitEnd; ++RecHit) {
0272         EBDetId EBid = EBDetId(RecHit->id());
0273 
0274         auto cellGeometry = geometry_->getSubdetectorGeometry(EBid)->getGeometry(EBid);
0275         double eta = cellGeometry->getPosition().eta();
0276         double phi = cellGeometry->getPosition().phi();
0277         double en = RecHit->energy();
0278 
0279         double r = dR(eta_MC, phi_MC, eta, phi);
0280         if (r < partR) {
0281           eEcalCone += en;
0282         }
0283       }
0284     }
0285 
0286     const edm::Handle<EERecHitCollection> &rhitEE = ev.getHandle(tok_EE_);
0287 
0288     if (rhitEE.isValid()) {
0289       RecHit = rhitEE.product()->begin();
0290       RecHitEnd = rhitEE.product()->end();
0291 
0292       for (; RecHit != RecHitEnd; ++RecHit) {
0293         EEDetId EEid = EEDetId(RecHit->id());
0294 
0295         auto cellGeometry = geometry_->getSubdetectorGeometry(EEid)->getGeometry(EEid);
0296         double eta = cellGeometry->getPosition().eta();
0297         double phi = cellGeometry->getPosition().phi();
0298         double en = RecHit->energy();
0299 
0300         double r = dR(eta_MC, phi_MC, eta, phi);
0301         if (r < partR) {
0302           eEcalCone += en;
0303         }
0304       }
0305     }
0306   }  // end of ECAL selection
0307 
0308 #ifdef EDM_ML_DEBUG
0309   edm::LogVerbatim("HcalRecHitsValidation") << "*** 4";
0310 #endif
0311 
0312   //===========================================================================
0313   // SUBSYSTEMS,
0314   //===========================================================================
0315 
0316   if ((subdet_ != 6) && (subdet_ != 0)) {
0317 #ifdef EDM_ML_DEBUG
0318     edm::LogVerbatim("HcalRecHitsValidation") << "*** 6";
0319 #endif
0320 
0321     double HcalCone = 0.;
0322 
0323     int ietaMax = 9999;
0324     double etaMax = 9999.;
0325 
0326     //   CYCLE over cells ====================================================
0327 
0328     for (unsigned int i = 0; i < cen.size(); i++) {
0329       int sub = csub[i];
0330       int depth = cdepth[i];
0331       double eta = ceta[i];
0332       double phi = cphi[i];
0333       double en = cen[i];
0334       double t = ctime[i];
0335       int ieta = cieta[i];
0336       double chi2 = cchi2[i];
0337 
0338       double chi2_log10 = 9.99;  // initial value above histos limits , keep it if chi2 <= 0.
0339       if (chi2 > 0.)
0340         chi2_log10 = log10(chi2);
0341 
0342       double r = dR(eta_MC, phi_MC, eta, phi);
0343       if (r < partR) {
0344         if (sub == 1)
0345           eHcalConeHB += en;
0346         if (sub == 2)
0347           eHcalConeHE += en;
0348         if (sub == 3)
0349           eHcalConeHO += en;
0350         if (sub == 4) {
0351           eHcalConeHF += en;
0352           if (depth == 1)
0353             eHcalConeHFL += en;
0354           else
0355             eHcalConeHFS += en;
0356         }
0357 
0358         HcalCone += en;
0359         // alternative: ietamax -> closest to MC eta  !!!
0360         float eta_diff = fabs(eta_MC - eta);
0361         if (eta_diff < etaMax) {
0362           etaMax = eta_diff;
0363           ietaMax = ieta;
0364         }
0365       }
0366 
0367       // The energy and overall timing histos are drawn while
0368       // the ones split by depth are not
0369       if (sub == 1 && (subdet_ == 1 || subdet_ == 5)) {
0370         meRecHitsM2Chi2HB->Fill(chi2_log10);
0371         meLog10Chi2profileHB->Fill(en, chi2_log10);
0372 
0373         meRecHitsEnergyHB->Fill(en);
0374 
0375         meTEprofileHB_Low->Fill(en, t);
0376         meTEprofileHB->Fill(en, t);
0377         meTEprofileHB_High->Fill(en, t);
0378       }
0379       if (sub == 2 && (subdet_ == 2 || subdet_ == 5)) {
0380         meRecHitsM2Chi2HE->Fill(chi2_log10);
0381         meLog10Chi2profileHE->Fill(en, chi2_log10);
0382 
0383         meRecHitsEnergyHE->Fill(en);
0384 
0385         meTEprofileHE_Low->Fill(en, t);
0386         meTEprofileHE->Fill(en, t);
0387       }
0388       if (sub == 4 && (subdet_ == 4 || subdet_ == 5)) {
0389         meRecHitsEnergyHF->Fill(en);
0390 
0391         meTEprofileHF_Low->Fill(en, t);
0392         meTEprofileHF->Fill(en, t);
0393       }
0394       if (sub == 3 && (subdet_ == 3 || subdet_ == 5)) {
0395         meRecHitsEnergyHO->Fill(en);
0396 
0397         meTEprofileHO->Fill(en, t);
0398         meTEprofileHO_High->Fill(en, t);
0399       }
0400     }
0401 
0402     if (imc != 0) {
0403       meEnConeEtaProfile->Fill(double(ietaMax), HcalCone);  //
0404       meEnConeEtaProfile_E->Fill(double(ietaMax), eEcalCone);
0405       meEnConeEtaProfile_EH->Fill(double(ietaMax), HcalCone + eEcalCone);
0406     }
0407 
0408 #ifdef EDM_ML_DEBUG
0409     edm::LogVerbatim("HcalRecHitsValidation") << "*** 7";
0410 #endif
0411   }
0412 
0413   // SimHits vs. RecHits
0414   if (subdet_ > 0 && subdet_ < 6 && imc != 0) {  // not noise
0415 
0416     const edm::Handle<PCaloHitContainer> &hcalHits = ev.getHandle(tok_hh_);
0417     if (hcalHits.isValid()) {
0418       const PCaloHitContainer *SimHitResult = hcalHits.product();
0419 
0420       double enSimHitsHB = 0.;
0421       double enSimHitsHE = 0.;
0422       double enSimHitsHO = 0.;
0423       double enSimHitsHF = 0.;
0424       double enSimHitsHFL = 0.;
0425       double enSimHitsHFS = 0.;
0426       // sum of SimHits in the cone
0427 
0428       for (std::vector<PCaloHit>::const_iterator SimHits = SimHitResult->begin(); SimHits != SimHitResult->end();
0429            ++SimHits) {
0430         int sub, depth;
0431         HcalDetId cell;
0432 
0433         if (testNumber_)
0434           cell = HcalHitRelabeller::relabel(SimHits->id(), hcons);
0435         else
0436           cell = HcalDetId(SimHits->id());
0437 
0438         sub = cell.subdet();
0439         depth = cell.depth();
0440 
0441         if (sub != subdet_ && subdet_ != 5)
0442           continue;  // If we are not looking at all of the subdetectors and the
0443                      // simhit doesn't come from the specific subdetector of
0444                      // interest, then we won't do any thing with it
0445 
0446         const HcalGeometry *cellGeometry =
0447             dynamic_cast<const HcalGeometry *>(geometry_->getSubdetectorGeometry(DetId::Hcal, cell.subdet()));
0448         double etaS = cellGeometry->getPosition(cell).eta();
0449         double phiS = cellGeometry->getPosition(cell).phi();
0450         double en = SimHits->energy();
0451 
0452         double r = dR(eta_MC, phi_MC, etaS, phiS);
0453 
0454         if (r < partR) {  // just energy in the small cone
0455 
0456           if (sub == static_cast<int>(HcalBarrel))
0457             enSimHitsHB += en;
0458           if (sub == static_cast<int>(HcalEndcap))
0459             enSimHitsHE += en;
0460           if (sub == static_cast<int>(HcalOuter))
0461             enSimHitsHO += en;
0462           if (sub == static_cast<int>(HcalForward)) {
0463             enSimHitsHF += en;
0464             if (depth == 1)
0465               enSimHitsHFL += en;
0466             else
0467               enSimHitsHFS += en;
0468           }
0469         }
0470       }
0471 
0472       // Now some histos with SimHits
0473 
0474       if (subdet_ == 4 || subdet_ == 5) {
0475         meRecHitSimHitHF->Fill(enSimHitsHF, eHcalConeHF);
0476         meRecHitSimHitProfileHF->Fill(enSimHitsHF, eHcalConeHF);
0477 
0478         meRecHitSimHitHFL->Fill(enSimHitsHFL, eHcalConeHFL);
0479         meRecHitSimHitProfileHFL->Fill(enSimHitsHFL, eHcalConeHFL);
0480         meRecHitSimHitHFS->Fill(enSimHitsHFS, eHcalConeHFS);
0481         meRecHitSimHitProfileHFS->Fill(enSimHitsHFS, eHcalConeHFS);
0482       }
0483       if (subdet_ == 1 || subdet_ == 5) {
0484         meRecHitSimHitHB->Fill(enSimHitsHB, eHcalConeHB);
0485         meRecHitSimHitProfileHB->Fill(enSimHitsHB, eHcalConeHB);
0486       }
0487       if (subdet_ == 2 || subdet_ == 5) {
0488         meRecHitSimHitHE->Fill(enSimHitsHE, eHcalConeHE);
0489         meRecHitSimHitProfileHE->Fill(enSimHitsHE, eHcalConeHE);
0490       }
0491       if (subdet_ == 3 || subdet_ == 5) {
0492         meRecHitSimHitHO->Fill(enSimHitsHO, eHcalConeHO);
0493         meRecHitSimHitProfileHO->Fill(enSimHitsHO, eHcalConeHO);
0494       }
0495     }
0496   }
0497 
0498   nevtot++;
0499 }
0500 
0501 ///////////////////////////////////////////////////////////////////////////////
0502 void HcalRecHitsValidation::fillRecHitsTmp(int subdet_, edm::Event const &ev) {
0503   using namespace edm;
0504 
0505   // initialize data vectors
0506   csub.clear();
0507   cen.clear();
0508   ceta.clear();
0509   cphi.clear();
0510   ctime.clear();
0511   cieta.clear();
0512   ciphi.clear();
0513   cdepth.clear();
0514   cz.clear();
0515   cchi2.clear();
0516 
0517   if (subdet_ == 1 || subdet_ == 2 || subdet_ == 5 || subdet_ == 6 || subdet_ == 0) {
0518     // HBHE
0519     const edm::Handle<HBHERecHitCollection> &hbhecoll = ev.getHandle(tok_hbhe_);
0520     if (hbhecoll.isValid()) {
0521       for (HBHERecHitCollection::const_iterator j = hbhecoll->begin(); j != hbhecoll->end(); j++) {
0522         HcalDetId cell(j->id());
0523         const HcalGeometry *cellGeometry =
0524             dynamic_cast<const HcalGeometry *>(geometry_->getSubdetectorGeometry(DetId::Hcal, cell.subdet()));
0525         double eta = cellGeometry->getPosition(cell).eta();
0526         double phi = cellGeometry->getPosition(cell).phi();
0527         double zc = cellGeometry->getPosition(cell).z();
0528         int sub = cell.subdet();
0529         int depth = cell.depth();
0530         int inteta = cell.ieta();
0531         int intphi = cell.iphi() - 1;
0532         double en = j->energy();
0533         double t = j->time();
0534         double chi2 = j->chi2();
0535 
0536         if ((iz > 0 && eta > 0.) || (iz < 0 && eta < 0.) || iz == 0) {
0537           csub.push_back(sub);
0538           cen.push_back(en);
0539           ceta.push_back(eta);
0540           cphi.push_back(phi);
0541           ctime.push_back(t);
0542           cieta.push_back(inteta);
0543           ciphi.push_back(intphi);
0544           cdepth.push_back(depth);
0545           cz.push_back(zc);
0546           cchi2.push_back(chi2);
0547         }
0548       }
0549     }
0550   }
0551 
0552   if (subdet_ == 4 || subdet_ == 5 || subdet_ == 6 || subdet_ == 0) {
0553     // HF
0554     const edm::Handle<HFRecHitCollection> &hfcoll = ev.getHandle(tok_hf_);
0555     if (hfcoll.isValid()) {
0556       for (HFRecHitCollection::const_iterator j = hfcoll->begin(); j != hfcoll->end(); j++) {
0557         HcalDetId cell(j->id());
0558         auto cellGeometry = geometry_->getSubdetectorGeometry(cell)->getGeometry(cell);
0559 
0560         double eta = cellGeometry->getPosition().eta();
0561         double phi = cellGeometry->getPosition().phi();
0562         double zc = cellGeometry->getPosition().z();
0563         int sub = cell.subdet();
0564         int depth = cell.depth();
0565         int inteta = cell.ieta();
0566         int intphi = cell.iphi() - 1;
0567         double en = j->energy();
0568         double t = j->time();
0569 
0570         if ((iz > 0 && eta > 0.) || (iz < 0 && eta < 0.) || iz == 0) {
0571           csub.push_back(sub);
0572           cen.push_back(en);
0573           ceta.push_back(eta);
0574           cphi.push_back(phi);
0575           ctime.push_back(t);
0576           cieta.push_back(inteta);
0577           ciphi.push_back(intphi);
0578           cdepth.push_back(depth);
0579           cz.push_back(zc);
0580           cchi2.push_back(0.);
0581         }
0582       }
0583     }
0584   }
0585 
0586   // HO
0587   if (subdet_ == 3 || subdet_ == 5 || subdet_ == 6 || subdet_ == 0) {
0588     const edm::Handle<HORecHitCollection> &hocoll = ev.getHandle(tok_ho_);
0589     if (hocoll.isValid()) {
0590       for (HORecHitCollection::const_iterator j = hocoll->begin(); j != hocoll->end(); j++) {
0591         HcalDetId cell(j->id());
0592         auto cellGeometry = geometry_->getSubdetectorGeometry(cell)->getGeometry(cell);
0593 
0594         double eta = cellGeometry->getPosition().eta();
0595         double phi = cellGeometry->getPosition().phi();
0596         double zc = cellGeometry->getPosition().z();
0597         int sub = cell.subdet();
0598         int depth = cell.depth();
0599         int inteta = cell.ieta();
0600         int intphi = cell.iphi() - 1;
0601         double t = j->time();
0602         double en = j->energy();
0603 
0604         if ((iz > 0 && eta > 0.) || (iz < 0 && eta < 0.) || iz == 0) {
0605           csub.push_back(sub);
0606           cen.push_back(en);
0607           ceta.push_back(eta);
0608           cphi.push_back(phi);
0609           ctime.push_back(t);
0610           cieta.push_back(inteta);
0611           ciphi.push_back(intphi);
0612           cdepth.push_back(depth);
0613           cz.push_back(zc);
0614           cchi2.push_back(0.);
0615         }
0616       }
0617     }
0618   }
0619 }
0620 
0621 double HcalRecHitsValidation::dR(double eta1, double phi1, double eta2, double phi2) {
0622   double PI = 3.1415926535898;
0623   double deltaphi = phi1 - phi2;
0624   if (phi2 > phi1) {
0625     deltaphi = phi2 - phi1;
0626   }
0627   if (deltaphi > PI) {
0628     deltaphi = 2. * PI - deltaphi;
0629   }
0630   double deltaeta = eta2 - eta1;
0631   double tmp = sqrt(deltaeta * deltaeta + deltaphi * deltaphi);
0632   return tmp;
0633 }
0634 
0635 double HcalRecHitsValidation::phi12(double phi1, double en1, double phi2, double en2) {
0636   // weighted mean value of phi1 and phi2
0637 
0638   double tmp;
0639   double PI = 3.1415926535898;
0640   double a1 = phi1;
0641   double a2 = phi2;
0642 
0643   if (a1 > 0.5 * PI && a2 < 0.)
0644     a2 += 2 * PI;
0645   if (a2 > 0.5 * PI && a1 < 0.)
0646     a1 += 2 * PI;
0647   tmp = (a1 * en1 + a2 * en2) / (en1 + en2);
0648   if (tmp > PI)
0649     tmp -= 2. * PI;
0650 
0651   return tmp;
0652 }
0653 
0654 double HcalRecHitsValidation::dPhiWsign(double phi1, double phi2) {
0655   // clockwise      phi2 w.r.t phi1 means "+" phi distance
0656   // anti-clockwise phi2 w.r.t phi1 means "-" phi distance
0657 
0658   double PI = 3.1415926535898;
0659   double a1 = phi1;
0660   double a2 = phi2;
0661   double tmp = a2 - a1;
0662   if (a1 * a2 < 0.) {
0663     if (a1 > 0.5 * PI)
0664       tmp += 2. * PI;
0665     if (a2 > 0.5 * PI)
0666       tmp -= 2. * PI;
0667   }
0668   return tmp;
0669 }
0670 
0671 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!