OfflineValidation/plugins/DiMuonValidation.cc

0001 // system includes
0002 #include <memory>
0003 #include <cmath>
0004 #include <fmt/printf.h>
0005
0006 // user include files
0007 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0008 #include "DataFormats/Candidate/interface/Particle.h"
0009 #include "DataFormats/TrackReco/interface/Track.h"
0010 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0011 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/Framework/interface/Frameworkfwd.h"
0014 #include "FWCore/Framework/interface/MakerMacros.h"
0015 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0017 #include "FWCore/ServiceRegistry/interface/Service.h"
0018 #include "FWCore/Utilities/interface/InputTag.h"
0019 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0020 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0021 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0022
0023 // ROOT includes
0024 #include "TH1D.h"
0025 #include "TH2D.h"
0026 #include "TH3D.h"
0027
0028 namespace DiMuonValid {
0029   using LV = reco::Particle::LorentzVector;
0030 }
0031 //
0032 // class declaration
0033 //
0034 class DiMuonValidation : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0035 public:
0036   explicit DiMuonValidation(const edm::ParameterSet& pset)
0037       : eBeam_(pset.getParameter<double>("eBeam")),
0038         compressionSettings_(pset.getUntrackedParameter<int>("compressionSettings", -1)),
0039         TkTag_(pset.getParameter<std::string>("TkTag")),
0040         pair_mass_min_(pset.getParameter<double>("Pair_mass_min")),
0041         pair_mass_max_(pset.getParameter<double>("Pair_mass_max")),
0042         pair_mass_nbins_(pset.getParameter<int>("Pair_mass_nbins")),
0043         pair_etaminpos_(pset.getParameter<double>("Pair_etaminpos")),
0044         pair_etamaxpos_(pset.getParameter<double>("Pair_etamaxpos")),
0045         pair_etaminneg_(pset.getParameter<double>("Pair_etaminneg")),
0046         pair_etamaxneg_(pset.getParameter<double>("Pair_etamaxneg")),
0047         variable_CosThetaCS_xmin_(pset.getParameter<double>("Variable_CosThetaCS_xmin")),
0048         variable_CosThetaCS_xmax_(pset.getParameter<double>("Variable_CosThetaCS_xmax")),
0049         variable_CosThetaCS_nbins_(pset.getParameter<int>("Variable_CosThetaCS_nbins")),
0050         variable_DeltaEta_xmin_(pset.getParameter<double>("Variable_DeltaEta_xmin")),
0051         variable_DeltaEta_xmax_(pset.getParameter<double>("Variable_DeltaEta_xmax")),
0052         variable_DeltaEta_nbins_(pset.getParameter<int>("Variable_DeltaEta_nbins")),
0053         variable_EtaMinus_xmin_(pset.getParameter<double>("Variable_EtaMinus_xmin")),
0054         variable_EtaMinus_xmax_(pset.getParameter<double>("Variable_EtaMinus_xmax")),
0055         variable_EtaMinus_nbins_(pset.getParameter<int>("Variable_EtaMinus_nbins")),
0056         variable_EtaPlus_xmin_(pset.getParameter<double>("Variable_EtaPlus_xmin")),
0057         variable_EtaPlus_xmax_(pset.getParameter<double>("Variable_EtaPlus_xmax")),
0058         variable_EtaPlus_nbins_(pset.getParameter<int>("Variable_EtaPlus_nbins")),
0059         variable_PhiCS_xmin_(pset.getParameter<double>("Variable_PhiCS_xmin")),
0060         variable_PhiCS_xmax_(pset.getParameter<double>("Variable_PhiCS_xmax")),
0061         variable_PhiCS_nbins_(pset.getParameter<int>("Variable_PhiCS_nbins")),
0062         variable_PhiMinus_xmin_(pset.getParameter<double>("Variable_PhiMinus_xmin")),
0063         variable_PhiMinus_xmax_(pset.getParameter<double>("Variable_PhiMinus_xmax")),
0064         variable_PhiMinus_nbins_(pset.getParameter<int>("Variable_PhiMinus_nbins")),
0065         variable_PhiPlus_xmin_(pset.getParameter<double>("Variable_PhiPlus_xmin")),
0066         variable_PhiPlus_xmax_(pset.getParameter<double>("Variable_PhiPlus_xmax")),
0067         variable_PhiPlus_nbins_(pset.getParameter<int>("Variable_PhiPlus_nbins")),
0068         variable_PairPt_xmin_(pset.getParameter<double>("Variable_PairPt_xmin")),
0069         variable_PairPt_xmax_(pset.getParameter<double>("Variable_PairPt_xmax")),
0070         variable_PairPt_nbins_(pset.getParameter<int>("Variable_PairPt_nbins")) {
0071     usesResource(TFileService::kSharedResource);
0072     theTrackCollectionToken_ = consumes<reco::TrackCollection>(TkTag_);
0073
0074     variables_min_[0] = variable_CosThetaCS_xmin_;
0075     variables_min_[1] = variable_DeltaEta_xmin_;
0076     variables_min_[2] = variable_EtaMinus_xmin_;
0077     variables_min_[3] = variable_EtaPlus_xmin_;
0078     variables_min_[4] = variable_PhiCS_xmin_;
0079     variables_min_[5] = variable_PhiMinus_xmin_;
0080     variables_min_[6] = variable_PhiPlus_xmin_;
0081     variables_min_[7] = variable_PairPt_xmin_;
0082
0083     variables_max_[0] = variable_CosThetaCS_xmax_;
0084     variables_max_[1] = variable_DeltaEta_xmax_;
0085     variables_max_[2] = variable_EtaMinus_xmax_;
0086     variables_max_[3] = variable_EtaPlus_xmax_;
0087     variables_max_[4] = variable_PhiCS_xmax_;
0088     variables_max_[5] = variable_PhiMinus_xmax_;
0089     variables_max_[6] = variable_PhiPlus_xmax_;
0090     variables_max_[7] = variable_PairPt_xmax_;
0091
0092     variables_bins_number_[0] = variable_CosThetaCS_nbins_;
0093     variables_bins_number_[1] = variable_DeltaEta_nbins_;
0094     variables_bins_number_[2] = variable_EtaMinus_nbins_;
0095     variables_bins_number_[3] = variable_EtaPlus_nbins_;
0096     variables_bins_number_[4] = variable_PhiCS_nbins_;
0097     variables_bins_number_[5] = variable_PhiMinus_nbins_;
0098     variables_bins_number_[6] = variable_PhiPlus_nbins_;
0099     variables_bins_number_[7] = variable_PairPt_nbins_;
0100   }
0101
0102   ~DiMuonValidation() override = default;
0103
0104   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0105   static constexpr int varNumber_ = 8;
0106   static constexpr double mu_mass2_ = 0.105658 * 0.105658;  //The invariant mass of muon 105.658MeV
0107
0108 private:
0109   void beginJob() override;
0110   void analyze(const edm::Event&, const edm::EventSetup&) override;
0111
0112   // ----------member data ---------------------------
0113   float eBeam_;
0114   int compressionSettings_;
0115   std::string TkTag_;
0116
0117   double pair_mass_min_;
0118   double pair_mass_max_;
0119   int pair_mass_nbins_;
0120   double pair_etaminpos_;
0121   double pair_etamaxpos_;
0122   double pair_etaminneg_;
0123   double pair_etamaxneg_;
0124
0125   double variable_CosThetaCS_xmin_;
0126   double variable_CosThetaCS_xmax_;
0127   int variable_CosThetaCS_nbins_;
0128
0129   double variable_DeltaEta_xmin_;
0130   double variable_DeltaEta_xmax_;
0131   int variable_DeltaEta_nbins_;
0132
0133   double variable_EtaMinus_xmin_;
0134   double variable_EtaMinus_xmax_;
0135   int variable_EtaMinus_nbins_;
0136
0137   double variable_EtaPlus_xmin_;
0138   double variable_EtaPlus_xmax_;
0139   int variable_EtaPlus_nbins_;
0140
0141   double variable_PhiCS_xmin_;
0142   double variable_PhiCS_xmax_;
0143   int variable_PhiCS_nbins_;
0144
0145   double variable_PhiMinus_xmin_;
0146   double variable_PhiMinus_xmax_;
0147   int variable_PhiMinus_nbins_;
0148
0149   double variable_PhiPlus_xmin_;
0150   double variable_PhiPlus_xmax_;
0151   int variable_PhiPlus_nbins_;
0152
0153   double variable_PairPt_xmin_;
0154   double variable_PairPt_xmax_;
0155   int variable_PairPt_nbins_;
0156
0157   edm::EDGetTokenT<reco::TrackCollection> theTrackCollectionToken_;
0158   TH2D* th2d_mass_variables_[varNumber_];  // actual histograms
0159   std::string variables_name_[varNumber_] = {
0160       "CosThetaCS", "DeltaEta", "EtaMinus", "EtaPlus", "PhiCS", "PhiMinus", "PhiPlus", "Pt"};
0161
0162   int variables_bins_number_[varNumber_];  // = {20, 20, 12, 12, 20, 16, 16, 100};
0163   double variables_min_[varNumber_];       // = {-1, -4.8, -2.4, -2.4, -M_PI / 2, -M_PI, -M_PI, 0};
0164   double variables_max_[varNumber_];       // = {+1, +4.8, +2.4, +2.4, +M_PI / 2, +M_PI, +M_PI, 100};
0165 };
0166
0167 //
0168 // member functions
0169 //
0170
0171 // ------------ method called for each event  ------------
0172 void DiMuonValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0173   using namespace edm;
0174
0175   const reco::TrackCollection& tC = iEvent.get(theTrackCollectionToken_);
0176
0177   DiMuonValid::LV LV_mother(0., 0., 0., 0.);
0178   //for (reco::TrackCollection::const_iterator track1 = tC.begin(); track1 != tC.end(); track1++) {
0179   for (const auto& track1 : tC) {
0180     DiMuonValid::LV LV_track1(track1.px(),
0181                               track1.py(),
0182                               track1.pz(),
0183                               sqrt((track1.p() * track1.p()) + mu_mass2_));  //old 106
0184
0185     for (const auto& track2 : tC) {
0186       if (&track1 == &track2) {
0187         continue;
0188       }  // discard the same track
0189
0190       if (track1.charge() == track2.charge()) {
0191         continue;
0192       }  // only reconstruct opposite charge pair
0193
0194       DiMuonValid::LV LV_track2(track2.px(), track2.py(), track2.pz(), sqrt((track2.p() * track2.p()) + mu_mass2_));
0195
0196       LV_mother = LV_track1 + LV_track2;
0197       double mother_mass = LV_mother.M();
0198       double mother_pt = LV_mother.Pt();
0199
0200       int charge1 = track1.charge();
0201       double etaMu1 = track1.eta();
0202       double phiMu1 = track1.phi();
0203       double ptMu1 = track1.pt();
0204
0205       int charge2 = track2.charge();
0206       double etaMu2 = track2.eta();
0207       double phiMu2 = track2.phi();
0208       double ptMu2 = track2.pt();
0209
0210       if (charge1 < 0) {  // use Mu+ for charge1, Mu- for charge2
0211         std::swap(charge1, charge2);
0212         std::swap(etaMu1, etaMu2);
0213         std::swap(phiMu1, phiMu2);
0214         std::swap(ptMu1, ptMu2);
0215       }
0216       //eta cut
0217       if (etaMu1 < pair_etaminpos_ || etaMu1 > pair_etamaxpos_ || etaMu2 < pair_etaminneg_ ||
0218           etaMu2 > pair_etamaxneg_) {
0219         continue;
0220       }
0221
0222       double delta_eta = etaMu1 - etaMu2;
0223
0224       double muplus = 1.0 / sqrt(2.0) * (LV_track1.E() + LV_track1.Z());
0225       double muminus = 1.0 / sqrt(2.0) * (LV_track1.E() - LV_track1.Z());
0226       double mubarplus = 1.0 / sqrt(2.0) * (LV_track2.E() + LV_track2.Z());
0227       double mubarminus = 1.0 / sqrt(2.0) * (LV_track2.E() - LV_track2.Z());
0228       //double costheta = 2.0 / Q.mag() / sqrt(pow(Q.mag(), 2) + pow(Q.Pt(), 2)) * (muplus * mubarminus - muminus * mubarplus);
0229       double costhetaCS = 2.0 / LV_mother.mag() / sqrt(pow(LV_mother.mag(), 2) + pow(LV_mother.Pt(), 2)) *
0230                           (muplus * mubarminus - muminus * mubarplus);
0231
0232       DiMuonValid::LV Pbeam(0., 0., eBeam_, eBeam_);
0233       auto R = Pbeam.Vect().Cross(LV_mother.Vect());
0234       auto Runit = R.Unit();
0235       auto Qt = LV_mother.Vect();
0236       Qt.SetZ(0);
0237       auto Qtunit = Qt.Unit();
0238       DiMuonValid::LV D(LV_track1 - LV_track2);
0239       auto Dt = D.Vect();
0240       Dt.SetZ(0);
0241       double tanphi =
0242           sqrt(pow(LV_mother.mag(), 2) + pow(LV_mother.Pt(), 2)) / LV_mother.mag() * Dt.Dot(Runit) / Dt.Dot(Qtunit);
0243       double phiCS = atan(tanphi);
0244
0245       if (mother_mass > pair_mass_min_ && mother_mass < pair_mass_max_) {
0246         th2d_mass_variables_[0]->Fill(mother_mass, costhetaCS, 1);
0247         th2d_mass_variables_[1]->Fill(mother_mass, delta_eta, 1);
0248         th2d_mass_variables_[2]->Fill(mother_mass, etaMu2, 1);
0249         th2d_mass_variables_[3]->Fill(mother_mass, etaMu1, 1);
0250         th2d_mass_variables_[4]->Fill(mother_mass, phiCS, 1);
0251         th2d_mass_variables_[5]->Fill(mother_mass, phiMu2, 1);
0252         th2d_mass_variables_[6]->Fill(mother_mass, phiMu1, 1);
0253         th2d_mass_variables_[7]->Fill(mother_mass, mother_pt, 1);
0254       }
0255     }
0256   }
0257 }
0258
0259 // ------------ method called once each job just before starting event loop  ------------
0260 void DiMuonValidation::beginJob() {
0261   edm::Service<TFileService> fs;
0262   if (compressionSettings_ > 0) {
0263     fs->file().SetCompressionSettings(compressionSettings_);
0264   }
0265
0266   for (int i = 0; i < varNumber_; i++) {
0267     std::string th2d_name = fmt::sprintf("th2d_mass_%s", variables_name_[i].c_str());
0268     th2d_mass_variables_[i] = fs->make<TH2D>(th2d_name.c_str(),
0269                                              th2d_name.c_str(),
0270                                              pair_mass_nbins_,
0271                                              pair_mass_min_,
0272                                              pair_mass_max_,
0273                                              variables_bins_number_[i],
0274                                              variables_min_[i],
0275                                              variables_max_[i]);
0276   }
0277 }
0278
0279 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0280 void DiMuonValidation::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0281   edm::ParameterSetDescription desc;
0282   desc.setComment("Validates alignment payloads by evaluating bias in Z->mm mass distributions");
0283   desc.addUntracked<int>("compressionSettings", -1);
0284
0285   desc.add<double>("eBeam", 3500.)->setComment("beam energy in GeV");
0286   desc.add<std::string>("TkTag", "ALCARECOTkAlZMuMu");
0287
0288   desc.add<double>("Pair_mass_min", 60);
0289   desc.add<double>("Pair_mass_max", 120);
0290   desc.add<int>("Pair_mass_nbins", 120);
0291   desc.add<double>("Pair_etaminpos", 60);
0292   desc.add<double>("Pair_etamaxpos", 60);
0293   desc.add<double>("Pair_etaminneg", 60);
0294   desc.add<double>("Pair_etamaxneg", 60);
0295
0296   desc.add<double>("Variable_CosThetaCS_xmin", -1.);
0297   desc.add<double>("Variable_CosThetaCS_xmax", 1.);
0298   desc.add<int>("Variable_CosThetaCS_nbins", 20);
0299
0300   desc.add<double>("Variable_DeltaEta_xmin", -4.8);
0301   desc.add<double>("Variable_DeltaEta_xmax", 4.8);
0302   desc.add<int>("Variable_DeltaEta_nbins", 20);
0303
0304   desc.add<double>("Variable_EtaMinus_xmin", -2.4);
0305   desc.add<double>("Variable_EtaMinus_xmax", 2.4);
0306   desc.add<int>("Variable_EtaMinus_nbins", 12);
0307
0308   desc.add<double>("Variable_EtaPlus_xmin", -2.4);
0309   desc.add<double>("Variable_EtaPlus_xmax", 2.4);
0310   desc.add<int>("Variable_EtaPlus_nbins", 12);
0311
0312   desc.add<double>("Variable_PhiCS_xmin", -M_PI / 2);
0313   desc.add<double>("Variable_PhiCS_xmax", M_PI / 2);
0314   desc.add<int>("Variable_PhiCS_nbins", 20);
0315
0316   desc.add<double>("Variable_PhiMinus_xmin", -M_PI);
0317   desc.add<double>("Variable_PhiMinus_xmax", M_PI);
0318   desc.add<int>("Variable_PhiMinus_nbins", 16);
0319
0320   desc.add<double>("Variable_PhiPlus_xmin", -M_PI);
0321   desc.add<double>("Variable_PhiPlus_xmax", M_PI);
0322   desc.add<int>("Variable_PhiPlus_nbins", 16);
0323
0324   desc.add<double>("Variable_PairPt_xmin", 0.);
0325   desc.add<double>("Variable_PairPt_xmax", 100.);
0326   desc.add<int>("Variable_PairPt_nbins", 100);
0327
0328   descriptions.addWithDefaultLabel(desc);
0329 }
0330
0331 //define this as a plug-in
0332 DEFINE_FWK_MODULE(DiMuonValidation);