Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-03-24 05:59:48 

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