Back to home page

Project CMSSW displayed by LXR

 
 

    


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

0001 #include "Validation/MuonME0Validation/interface/ME0DigisValidation.h"
0002 #include <TMath.h>
0003 
0004 ME0DigisValidation::ME0DigisValidation(const edm::ParameterSet &cfg) : ME0BaseValidation(cfg) {
0005   InputTagToken_ = consumes<edm::PSimHitContainer>(cfg.getParameter<edm::InputTag>("simInputLabel"));
0006   InputTagToken_Digi = consumes<ME0DigiPreRecoCollection>(cfg.getParameter<edm::InputTag>("digiInputLabel"));
0007   sigma_x_ = cfg.getParameter<double>("sigma_x");
0008   sigma_y_ = cfg.getParameter<double>("sigma_y");
0009 }
0010 
0011 void ME0DigisValidation::bookHistograms(DQMStore::IBooker &ibooker,
0012                                         edm::Run const &Run,
0013                                         edm::EventSetup const &iSetup) {
0014   LogDebug("MuonME0DigisValidation") << "Info: Loading Geometry information\n";
0015   ibooker.setCurrentFolder("MuonME0DigisV/ME0DigisTask");
0016 
0017   unsigned int nregion = 2;
0018 
0019   edm::LogInfo("MuonME0DigisValidation") << "+++ Info : # of region: " << nregion << std::endl;
0020 
0021   LogDebug("MuonME0DigisValidation") << "+++ Info : finish to get geometry information from ES.\n";
0022 
0023   num_evts = ibooker.book1D("num_evts", "Number of events; ; Number of events", 1, 0, 2);
0024 
0025   me0_strip_dg_x_local_tot =
0026       ibooker.book1D("me0_strip_dg_x_local_tot", "Local X; X_{local} [cm]; Entries", 60, -30.0, +30.0);
0027   me0_strip_dg_y_local_tot =
0028       ibooker.book1D("me0_strip_dg_y_local_tot", "Local Y; Y_{local} [cm]; Entries", 100, -50.0, +50.0);
0029   me0_strip_dg_time_tot = ibooker.book1D("me0_strip_dg_time_tot", "ToF; ToF [ns]; Entries", 400, -200, +200);
0030 
0031   me0_strip_dg_dx_local_tot_Muon =
0032       ibooker.book1D("me0_strip_dg_dx_local_tot", "Local DeltaX; #Delta X_{local} [cm]; Entries", 50, -0.1, +0.1);
0033   me0_strip_dg_dy_local_tot_Muon =
0034       ibooker.book1D("me0_strip_dg_dy_local_tot", "Local DeltaY; #Delta Y_{local} [cm]; Entries", 500, -10.0, +10.0);
0035   me0_strip_dg_dphi_global_tot_Muon = ibooker.book1D(
0036       "me0_strip_dg_dphi_global_tot", "Global DeltaPhi; #Delta #phi_{global} [rad]; Entries", 50, -0.01, +0.01);
0037   me0_strip_dg_dtime_tot_Muon =
0038       ibooker.book1D("me0_strip_dg_dtime_tot", "DeltaToF; #Delta ToF [ns]; Entries", 50, -5, +5);
0039 
0040   me0_strip_dg_dphi_vs_phi_global_tot_Muon = ibooker.book2D("me0_strip_dg_dphi_vs_phi_global_tot",
0041                                                             "Global DeltaPhi vs. Phi; #phi_{global} [rad]; #Delta "
0042                                                             "#phi_{global} [rad]",
0043                                                             72,
0044                                                             -M_PI,
0045                                                             +M_PI,
0046                                                             50,
0047                                                             -0.01,
0048                                                             +0.01);
0049 
0050   me0_strip_dg_den_eta_tot = ibooker.book1D("me0_strip_dg_den_eta_tot", "Denominator; #eta; Entries", 12, 1.8, 3.0);
0051   me0_strip_dg_num_eta_tot = ibooker.book1D("me0_strip_dg_num_eta_tot", "Numerator; #eta; Entries", 12, 1.8, 3.0);
0052 
0053   float bins[] = {62.3, 70.0, 77.7, 87.1, 96.4, 108.2, 119.9, 134.7, 149.5};
0054   int binnum = sizeof(bins) / sizeof(float) - 1;
0055 
0056   me0_strip_dg_bkg_rad_tot =
0057       ibooker.book1D("me0_strip_dg_bkg_radius_tot", "Total neutron background; Radius [cm]; Entries", binnum, bins);
0058   me0_strip_dg_bkgElePos_rad = ibooker.book1D(
0059       "me0_strip_dg_bkgElePos_radius", "Neutron background: electrons+positrons; Radius [cm]; Entries", binnum, bins);
0060   me0_strip_dg_bkgNeutral_rad = ibooker.book1D(
0061       "me0_strip_dg_bkgNeutral_radius", "Neutron background: gammas+neutrons; Radius [cm]; Entries", binnum, bins);
0062 
0063   me0_strip_exp_bkg_rad_tot = ibooker.book1D("me0_strip_exp_bkg_radius_tot",
0064                                              "Total expected neutron background; Radius [cm]; Hit Rate [Hz/cm^{2}]",
0065                                              binnum,
0066                                              bins);
0067   me0_strip_exp_bkgElePos_rad = ibooker.book1D("me0_strip_exp_bkgElePos_radius",
0068                                                "Expected neutron background: electrons+positrons; Radius "
0069                                                "[cm]; Hit Rate [Hz/cm^{2}]",
0070                                                binnum,
0071                                                bins);
0072   me0_strip_exp_bkgNeutral_rad = ibooker.book1D("me0_strip_exp_bkgNeutral_radius",
0073                                                 "Expected neutron background: gammas+neutrons; Radius "
0074                                                 "[cm]; Hit Rate [Hz/cm^{2}]",
0075                                                 binnum,
0076                                                 bins);
0077 
0078   std::vector<double> neuBkg, eleBkg;
0079   neuBkg.push_back(899644.0);
0080   neuBkg.push_back(-30841.0);
0081   neuBkg.push_back(441.28);
0082   neuBkg.push_back(-3.3405);
0083   neuBkg.push_back(0.0140588);
0084   neuBkg.push_back(-3.11473e-05);
0085   neuBkg.push_back(2.83736e-08);
0086   eleBkg.push_back(4.68590e+05);
0087   eleBkg.push_back(-1.63834e+04);
0088   eleBkg.push_back(2.35700e+02);
0089   eleBkg.push_back(-1.77706e+00);
0090   eleBkg.push_back(7.39960e-03);
0091   eleBkg.push_back(-1.61448e-05);
0092   eleBkg.push_back(1.44368e-08);
0093 
0094   for (int i = 1; i < me0_strip_exp_bkgNeutral_rad->getTH1F()->GetSize(); i++) {
0095     double pos = me0_strip_exp_bkgNeutral_rad->getTH1F()->GetBinCenter(i);
0096     double neutral = 0;
0097     double charged = 0;
0098 
0099     double pos_helper = 1.0;
0100     for (int j = 0; j < 7; ++j) {
0101       neutral += neuBkg[j] * pos_helper;
0102       charged += eleBkg[j] * pos_helper;
0103       pos_helper *= pos;
0104     }
0105 
0106     me0_strip_exp_bkgNeutral_rad->setBinContent(i, neutral);
0107     me0_strip_exp_bkgElePos_rad->setBinContent(i, charged);
0108     me0_strip_exp_bkg_rad_tot->setBinContent(i, neutral + charged);
0109   }
0110 
0111   for (unsigned int region_num = 0; region_num < nregion; region_num++) {
0112     me0_strip_dg_zr_tot[region_num] = BookHistZR(ibooker, "me0_strip_dg_tot", "Digi", region_num);
0113     me0_strip_dg_zr_tot_Muon[region_num] = BookHistZR(ibooker, "me0_strip_dg_tot_Muon", "Digi Muon", region_num);
0114     for (unsigned int layer_num = 0; layer_num < 6; layer_num++) {
0115       std::string hist_name_for_dx_local =
0116           std::string("me0_strip_dg_dx_local") + regionLabel[region_num] + "_l" + layerLabel[layer_num];
0117       std::string hist_name_for_dy_local =
0118           std::string("me0_strip_dg_dy_local") + regionLabel[region_num] + "_l" + layerLabel[layer_num];
0119       std::string hist_name_for_dphi_global =
0120           std::string("me0_strip_dg_dphi_global") + regionLabel[region_num] + "_l" + layerLabel[layer_num];
0121 
0122       std::string hist_name_for_den_eta =
0123           std::string("me0_strip_dg_den_eta") + regionLabel[region_num] + "_l" + layerLabel[layer_num];
0124       std::string hist_name_for_num_eta =
0125           std::string("me0_strip_dg_num_eta") + regionLabel[region_num] + "_l" + layerLabel[layer_num];
0126 
0127       std::string hist_label_for_dx_local = "Local DeltaX: region" + regionLabel[region_num] + " layer " +
0128                                             layerLabel[layer_num] + " ; #Delta X_{local} [cm]; Entries";
0129       std::string hist_label_for_dy_local = "Local DeltaY: region" + regionLabel[region_num] + " layer " +
0130                                             layerLabel[layer_num] + " ; #Delta Y_{local} [cm]; Entries";
0131       std::string hist_label_for_dphi_global = "Global DeltaPhi: region" + regionLabel[region_num] + " layer " +
0132                                                layerLabel[layer_num] + " ; #Delta #phi_{global} [rad]; Entries";
0133 
0134       std::string hist_label_for_den_eta =
0135           "Denominator: region" + regionLabel[region_num] + " layer " + layerLabel[layer_num] + " ; #eta; Entries";
0136       std::string hist_label_for_num_eta =
0137           "Numerator: region" + regionLabel[region_num] + " layer " + layerLabel[layer_num] + " ; #eta; Entries";
0138 
0139       me0_strip_dg_xy[region_num][layer_num] = BookHistXY(ibooker, "me0_strip_dg", "Digi", region_num, layer_num);
0140       me0_strip_dg_xy_Muon[region_num][layer_num] =
0141           BookHistXY(ibooker, "me0_strip_dg_Muon", "Digi Muon", region_num, layer_num);
0142 
0143       me0_strip_dg_dx_local_Muon[region_num][layer_num] =
0144           ibooker.book1D(hist_name_for_dx_local.c_str(), hist_label_for_dx_local.c_str(), 50, -0.1, +0.1);
0145       me0_strip_dg_dy_local_Muon[region_num][layer_num] =
0146           ibooker.book1D(hist_name_for_dy_local.c_str(), hist_label_for_dy_local.c_str(), 500, -10.0, +10.0);
0147       me0_strip_dg_dphi_global_Muon[region_num][layer_num] =
0148           ibooker.book1D(hist_name_for_dphi_global.c_str(), hist_label_for_dphi_global.c_str(), 50, -0.01, +0.01);
0149 
0150       me0_strip_dg_den_eta[region_num][layer_num] =
0151           ibooker.book1D(hist_name_for_den_eta, hist_label_for_den_eta, 12, 1.8, 3.0);
0152       me0_strip_dg_num_eta[region_num][layer_num] =
0153           ibooker.book1D(hist_name_for_num_eta, hist_label_for_num_eta, 12, 1.8, 3.0);
0154     }
0155   }
0156 }
0157 
0158 ME0DigisValidation::~ME0DigisValidation() {}
0159 
0160 void ME0DigisValidation::analyze(const edm::Event &e, const edm::EventSetup &iSetup) {
0161   const ME0Geometry *ME0Geometry_ = &iSetup.getData(geomToken_);
0162 
0163   edm::Handle<edm::PSimHitContainer> ME0Hits;
0164   e.getByToken(InputTagToken_, ME0Hits);
0165 
0166   edm::Handle<ME0DigiPreRecoCollection> ME0Digis;
0167   e.getByToken(InputTagToken_Digi, ME0Digis);
0168 
0169   if (!ME0Hits.isValid() || !ME0Digis.isValid()) {
0170     edm::LogError("ME0DigisValidation") << "Cannot get ME0Hits/ME0Digis by Token simInputTagToken";
0171     return;
0172   }
0173 
0174   num_evts->Fill(1);
0175   bool toBeCounted1 = true;
0176 
0177   for (ME0DigiPreRecoCollection::DigiRangeIterator cItr = ME0Digis->begin(); cItr != ME0Digis->end(); cItr++) {
0178     ME0DetId id = (*cItr).first;
0179 
0180     const GeomDet *gdet = ME0Geometry_->idToDet(id);
0181     if (gdet == nullptr) {
0182       edm::LogWarning("ME0DigisValidation") << "Getting DetId failed. Discard this gem strip hit. Maybe it comes "
0183                                                "from unmatched geometry.";
0184       continue;
0185     }
0186     const BoundPlane &surface = gdet->surface();
0187 
0188     int region = (int)id.region();
0189     int layer = (int)id.layer();
0190     int chamber = (int)id.chamber();
0191     int roll = (int)id.roll();
0192 
0193     ME0DigiPreRecoCollection::const_iterator digiItr;
0194     for (digiItr = (*cItr).second.first; digiItr != (*cItr).second.second; ++digiItr) {
0195       LocalPoint lp(digiItr->x(), digiItr->y(), 0);
0196 
0197       GlobalPoint gp = surface.toGlobal(lp);
0198 
0199       float g_r = (float)gp.perp();
0200       float g_x = (float)gp.x();
0201       float g_y = (float)gp.y();
0202       float g_z = (float)gp.z();
0203 
0204       int particleType = digiItr->pdgid();
0205       int isPrompt = digiItr->prompt();
0206 
0207       float timeOfFlight = digiItr->tof();
0208 
0209       me0_strip_dg_x_local_tot->Fill(lp.x());
0210       me0_strip_dg_y_local_tot->Fill(lp.y());
0211       me0_strip_dg_time_tot->Fill(timeOfFlight);
0212 
0213       // fill hist
0214       int region_num = 0;
0215       if (region == -1)
0216         region_num = 0;
0217       else if (region == 1)
0218         region_num = 1;
0219       int layer_num = layer - 1;
0220 
0221       bool toBeCounted2 = true;
0222 
0223       if (isPrompt == 1 && abs(particleType) == 13) {
0224         me0_strip_dg_zr_tot_Muon[region_num]->Fill(g_z, g_r);
0225         me0_strip_dg_xy_Muon[region_num][layer_num]->Fill(g_x, g_y);
0226 
0227         for (auto hits = ME0Hits->begin(); hits != ME0Hits->end(); hits++) {
0228           int particleType_sh = hits->particleType();
0229           int evtId_sh = hits->eventId().event();
0230           int bx_sh = hits->eventId().bunchCrossing();
0231           int procType_sh = hits->processType();
0232 
0233           if (!(abs(particleType_sh) == 13 && evtId_sh == 0 && bx_sh == 0 && procType_sh == 0))
0234             continue;
0235 
0236           const ME0DetId id(hits->detUnitId());
0237           int region_sh = id.region();
0238           int layer_sh = id.layer();
0239           int chamber_sh = id.chamber();
0240           int roll_sh = id.roll();
0241 
0242           int region_sh_num = 0;
0243           if (region_sh == -1)
0244             region_sh_num = 0;
0245           else if (region_sh == 1)
0246             region_sh_num = 1;
0247           int layer_sh_num = layer_sh - 1;
0248 
0249           LocalPoint lp_sh = hits->localPosition();
0250           GlobalPoint gp_sh = ME0Geometry_->idToDet(id)->surface().toGlobal(lp_sh);
0251 
0252           if (toBeCounted1) {
0253             me0_strip_dg_den_eta[region_sh_num][layer_sh_num]->Fill(fabs(gp_sh.eta()));
0254             me0_strip_dg_den_eta_tot->Fill(fabs(gp_sh.eta()));
0255           }
0256 
0257           if (!(region == region_sh && layer == layer_sh && chamber == chamber_sh && roll == roll_sh))
0258             continue;
0259 
0260           float dx_loc = lp_sh.x() - lp.x();
0261           float dy_loc = lp_sh.y() - lp.y();
0262           float dphi_glob = gp_sh.phi() - gp.phi();
0263           float deta_glob = gp_sh.eta() - gp.eta();
0264 
0265           if (!(fabs(dphi_glob) < 3 * sigma_x_ && fabs(deta_glob) < 3 * sigma_y_))
0266             continue;
0267 
0268           float timeOfFlight_sh = hits->tof();
0269           const LocalPoint centralLP(0., 0., 0.);
0270           const GlobalPoint centralGP(ME0Geometry_->idToDet(id)->surface().toGlobal(centralLP));
0271           float centralTOF(centralGP.mag() / 29.98);  // speed of light
0272           float timeOfFlight_sh_corr = timeOfFlight_sh - centralTOF;
0273 
0274           me0_strip_dg_dx_local_Muon[region_num][layer_num]->Fill(dx_loc);
0275           me0_strip_dg_dy_local_Muon[region_num][layer_num]->Fill(dy_loc);
0276           me0_strip_dg_dphi_global_Muon[region_num][layer_num]->Fill(dphi_glob);
0277 
0278           me0_strip_dg_dx_local_tot_Muon->Fill(dx_loc);
0279           me0_strip_dg_dy_local_tot_Muon->Fill(dy_loc);
0280           me0_strip_dg_dphi_global_tot_Muon->Fill(dphi_glob);
0281 
0282           me0_strip_dg_dphi_vs_phi_global_tot_Muon->Fill(gp_sh.phi(), dphi_glob);
0283           me0_strip_dg_dtime_tot_Muon->Fill(timeOfFlight - timeOfFlight_sh_corr);
0284 
0285           if (toBeCounted2) {
0286             me0_strip_dg_num_eta[region_num][layer_num]->Fill(fabs(gp_sh.eta()));
0287             me0_strip_dg_num_eta_tot->Fill(fabs(gp_sh.eta()));
0288           }
0289           toBeCounted2 = false;
0290 
0291         }  // loop SH
0292 
0293         toBeCounted1 = false;
0294 
0295       } else {
0296         me0_strip_dg_zr_tot[region_num]->Fill(g_z, g_r);
0297         me0_strip_dg_xy[region_num][layer_num]->Fill(g_x, g_y);
0298       }
0299 
0300       if ((abs(particleType) == 11 || abs(particleType) == 22 || abs(particleType) == 2112) && isPrompt == 0)
0301         me0_strip_dg_bkg_rad_tot->Fill(fabs(gp.perp()));
0302       if ((abs(particleType) == 11) && isPrompt == 0)
0303         me0_strip_dg_bkgElePos_rad->Fill(fabs(gp.perp()));
0304       if ((abs(particleType) == 22 || abs(particleType) == 2112) && isPrompt == 0)
0305         me0_strip_dg_bkgNeutral_rad->Fill(fabs(gp.perp()));
0306 
0307     }  // loop DG
0308   }
0309 }