CMS 3D CMS Logo

DiMuonValidation.cc
Go to the documentation of this file.
1 // system includes
2 #include <memory>
3 #include <cmath>
4 #include <fmt/printf.h>
5 
6 // user include files
22 
23 // ROOT includes
24 #include "TH1D.h"
25 #include "TH2D.h"
26 #include "TH3D.h"
27 
28 namespace DiMuonValid {
30 }
31 //
32 // class declaration
33 //
34 class DiMuonValidation : public edm::one::EDAnalyzer<edm::one::SharedResources> {
35 public:
37  : eBeam_(pset.getParameter<double>("eBeam")),
38  compressionSettings_(pset.getUntrackedParameter<int>("compressionSettings", -1)),
39  TkTag_(pset.getParameter<std::string>("TkTag")),
40  pair_mass_min_(pset.getParameter<double>("Pair_mass_min")),
41  pair_mass_max_(pset.getParameter<double>("Pair_mass_max")),
42  pair_mass_nbins_(pset.getParameter<int>("Pair_mass_nbins")),
43  pair_etaminpos_(pset.getParameter<double>("Pair_etaminpos")),
44  pair_etamaxpos_(pset.getParameter<double>("Pair_etamaxpos")),
45  pair_etaminneg_(pset.getParameter<double>("Pair_etaminneg")),
46  pair_etamaxneg_(pset.getParameter<double>("Pair_etamaxneg")),
47  variable_CosThetaCS_xmin_(pset.getParameter<double>("Variable_CosThetaCS_xmin")),
48  variable_CosThetaCS_xmax_(pset.getParameter<double>("Variable_CosThetaCS_xmax")),
49  variable_CosThetaCS_nbins_(pset.getParameter<int>("Variable_CosThetaCS_nbins")),
50  variable_DeltaEta_xmin_(pset.getParameter<double>("Variable_DeltaEta_xmin")),
51  variable_DeltaEta_xmax_(pset.getParameter<double>("Variable_DeltaEta_xmax")),
52  variable_DeltaEta_nbins_(pset.getParameter<int>("Variable_DeltaEta_nbins")),
53  variable_EtaMinus_xmin_(pset.getParameter<double>("Variable_EtaMinus_xmin")),
54  variable_EtaMinus_xmax_(pset.getParameter<double>("Variable_EtaMinus_xmax")),
55  variable_EtaMinus_nbins_(pset.getParameter<int>("Variable_EtaMinus_nbins")),
56  variable_EtaPlus_xmin_(pset.getParameter<double>("Variable_EtaPlus_xmin")),
57  variable_EtaPlus_xmax_(pset.getParameter<double>("Variable_EtaPlus_xmax")),
58  variable_EtaPlus_nbins_(pset.getParameter<int>("Variable_EtaPlus_nbins")),
59  variable_PhiCS_xmin_(pset.getParameter<double>("Variable_PhiCS_xmin")),
60  variable_PhiCS_xmax_(pset.getParameter<double>("Variable_PhiCS_xmax")),
61  variable_PhiCS_nbins_(pset.getParameter<int>("Variable_PhiCS_nbins")),
62  variable_PhiMinus_xmin_(pset.getParameter<double>("Variable_PhiMinus_xmin")),
63  variable_PhiMinus_xmax_(pset.getParameter<double>("Variable_PhiMinus_xmax")),
64  variable_PhiMinus_nbins_(pset.getParameter<int>("Variable_PhiMinus_nbins")),
65  variable_PhiPlus_xmin_(pset.getParameter<double>("Variable_PhiPlus_xmin")),
66  variable_PhiPlus_xmax_(pset.getParameter<double>("Variable_PhiPlus_xmax")),
67  variable_PhiPlus_nbins_(pset.getParameter<int>("Variable_PhiPlus_nbins")),
68  variable_PairPt_xmin_(pset.getParameter<double>("Variable_PairPt_xmin")),
69  variable_PairPt_xmax_(pset.getParameter<double>("Variable_PairPt_xmax")),
70  variable_PairPt_nbins_(pset.getParameter<int>("Variable_PairPt_nbins")) {
71  usesResource(TFileService::kSharedResource);
72  theTrackCollectionToken_ = consumes<reco::TrackCollection>(TkTag_);
73 
82 
91 
100  }
101 
102  ~DiMuonValidation() override = default;
103 
104  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
105  static constexpr int varNumber_ = 8;
106  static constexpr double mu_mass2_ = 0.105658 * 0.105658; //The invariant mass of muon 105.658MeV
107 
108 private:
109  void beginJob() override;
110  void analyze(const edm::Event&, const edm::EventSetup&) override;
111 
112  // ----------member data ---------------------------
113  float eBeam_;
116 
124 
128 
132 
136 
140 
144 
148 
152 
156 
158 
159  TH1F* th1f_mass;
160  TH2D* th2d_mass_variables_[varNumber_]; // actual histograms
162  "CosThetaCS", "DeltaEta", "EtaMinus", "EtaPlus", "PhiCS", "PhiMinus", "PhiPlus", "Pt"};
163 
164  int variables_bins_number_[varNumber_]; // = {20, 20, 12, 12, 20, 16, 16, 100};
165  double variables_min_[varNumber_]; // = {-1, -4.8, -2.4, -2.4, -M_PI / 2, -M_PI, -M_PI, 0};
166  double variables_max_[varNumber_]; // = {+1, +4.8, +2.4, +2.4, +M_PI / 2, +M_PI, +M_PI, 100};
167 };
168 
169 //
170 // member functions
171 //
172 
173 // ------------ method called for each event ------------
175  using namespace edm;
176 
178 
179  DiMuonValid::LV LV_mother(0., 0., 0., 0.);
180  for (reco::TrackCollection::const_iterator track1 = tC.begin(); track1 != tC.end(); track1++) {
181  DiMuonValid::LV LV_track1(track1->px(),
182  track1->py(),
183  track1->pz(),
184  sqrt((track1->p() * track1->p()) + mu_mass2_)); //old 106
185 
186  for (reco::TrackCollection::const_iterator track2 = track1 + 1; track2 != tC.end(); track2++) {
187  if (track1->charge() == track2->charge()) {
188  continue;
189  } // only reconstruct opposite charge pair
190 
191  DiMuonValid::LV LV_track2(
192  track2->px(), track2->py(), track2->pz(), sqrt((track2->p() * track2->p()) + mu_mass2_));
193 
194  LV_mother = LV_track1 + LV_track2;
195  double mother_mass = LV_mother.M();
196  th1f_mass->Fill(mother_mass);
197 
198  double mother_pt = LV_mother.Pt();
199 
200  int charge1 = track1->charge();
201  double etaMu1 = track1->eta();
202  double phiMu1 = track1->phi();
203  double ptMu1 = track1->pt();
204 
205  int charge2 = track2->charge();
206  double etaMu2 = track2->eta();
207  double phiMu2 = track2->phi();
208  double ptMu2 = track2->pt();
209 
210  if (charge1 < 0) { // use Mu+ for charge1, Mu- for charge2
211  std::swap(charge1, charge2);
212  std::swap(etaMu1, etaMu2);
213  std::swap(phiMu1, phiMu2);
214  std::swap(ptMu1, ptMu2);
215  }
216  //eta cut
217  if (etaMu1 < pair_etaminpos_ || etaMu1 > pair_etamaxpos_ || etaMu2 < pair_etaminneg_ ||
218  etaMu2 > pair_etamaxneg_) {
219  continue;
220  }
221 
222  double delta_eta = etaMu1 - etaMu2;
223 
224  double muplus = 1.0 / sqrt(2.0) * (LV_track1.E() + LV_track1.Z());
225  double muminus = 1.0 / sqrt(2.0) * (LV_track1.E() - LV_track1.Z());
226  double mubarplus = 1.0 / sqrt(2.0) * (LV_track2.E() + LV_track2.Z());
227  double mubarminus = 1.0 / sqrt(2.0) * (LV_track2.E() - LV_track2.Z());
228  //double costheta = 2.0 / Q.mag() / sqrt(pow(Q.mag(), 2) + pow(Q.Pt(), 2)) * (muplus * mubarminus - muminus * mubarplus);
229  double costhetaCS = 2.0 / LV_mother.mag() / sqrt(pow(LV_mother.mag(), 2) + pow(LV_mother.Pt(), 2)) *
230  (muplus * mubarminus - muminus * mubarplus);
231 
232  DiMuonValid::LV Pbeam(0., 0., eBeam_, eBeam_);
233  auto R = Pbeam.Vect().Cross(LV_mother.Vect());
234  auto Runit = R.Unit();
235  auto Qt = LV_mother.Vect();
236  Qt.SetZ(0);
237  auto Qtunit = Qt.Unit();
238  DiMuonValid::LV D(LV_track1 - LV_track2);
239  auto Dt = D.Vect();
240  Dt.SetZ(0);
241  double tanphi =
242  sqrt(pow(LV_mother.mag(), 2) + pow(LV_mother.Pt(), 2)) / LV_mother.mag() * Dt.Dot(Runit) / Dt.Dot(Qtunit);
243  double phiCS = atan(tanphi);
244 
245  if (mother_mass > pair_mass_min_ && mother_mass < pair_mass_max_) {
246  th2d_mass_variables_[0]->Fill(mother_mass, costhetaCS, 1);
247  th2d_mass_variables_[1]->Fill(mother_mass, delta_eta, 1);
248  th2d_mass_variables_[2]->Fill(mother_mass, etaMu2, 1);
249  th2d_mass_variables_[3]->Fill(mother_mass, etaMu1, 1);
250  th2d_mass_variables_[4]->Fill(mother_mass, phiCS, 1);
251  th2d_mass_variables_[5]->Fill(mother_mass, phiMu2, 1);
252  th2d_mass_variables_[6]->Fill(mother_mass, phiMu1, 1);
253  th2d_mass_variables_[7]->Fill(mother_mass, mother_pt, 1);
254  }
255  }
256  }
257 }
258 
259 // ------------ method called once each job just before starting event loop ------------
262  if (compressionSettings_ > 0) {
263  fs->file().SetCompressionSettings(compressionSettings_);
264  }
265 
266  th1f_mass = fs->make<TH1F>("hMass", "mass;m_{#mu#mu} [GeV];events", 200, 0., 200.);
267 
268  for (int i = 0; i < varNumber_; i++) {
269  std::string th2d_name = fmt::sprintf("th2d_mass_%s", variables_name_[i].c_str());
270  th2d_mass_variables_[i] = fs->make<TH2D>(th2d_name.c_str(),
271  th2d_name.c_str(),
276  variables_min_[i],
277  variables_max_[i]);
278  }
279 }
280 
281 // ------------ method fills 'descriptions' with the allowed parameters for the module ------------
284  desc.setComment("Validates alignment payloads by evaluating bias in Z->mm mass distributions");
285  desc.addUntracked<int>("compressionSettings", -1);
286 
287  desc.add<double>("eBeam", 3500.)->setComment("beam energy in GeV");
288  desc.add<std::string>("TkTag", "ALCARECOTkAlZMuMu");
289 
290  desc.add<double>("Pair_mass_min", 60);
291  desc.add<double>("Pair_mass_max", 120);
292  desc.add<int>("Pair_mass_nbins", 120);
293  desc.add<double>("Pair_etaminpos", -2.4);
294  desc.add<double>("Pair_etamaxpos", 2.4);
295  desc.add<double>("Pair_etaminneg", -2.4);
296  desc.add<double>("Pair_etamaxneg", 2.4);
297 
298  desc.add<double>("Variable_CosThetaCS_xmin", -1.);
299  desc.add<double>("Variable_CosThetaCS_xmax", 1.);
300  desc.add<int>("Variable_CosThetaCS_nbins", 20);
301 
302  desc.add<double>("Variable_DeltaEta_xmin", -4.8);
303  desc.add<double>("Variable_DeltaEta_xmax", 4.8);
304  desc.add<int>("Variable_DeltaEta_nbins", 20);
305 
306  desc.add<double>("Variable_EtaMinus_xmin", -2.4);
307  desc.add<double>("Variable_EtaMinus_xmax", 2.4);
308  desc.add<int>("Variable_EtaMinus_nbins", 12);
309 
310  desc.add<double>("Variable_EtaPlus_xmin", -2.4);
311  desc.add<double>("Variable_EtaPlus_xmax", 2.4);
312  desc.add<int>("Variable_EtaPlus_nbins", 12);
313 
314  desc.add<double>("Variable_PhiCS_xmin", -M_PI / 2);
315  desc.add<double>("Variable_PhiCS_xmax", M_PI / 2);
316  desc.add<int>("Variable_PhiCS_nbins", 20);
317 
318  desc.add<double>("Variable_PhiMinus_xmin", -M_PI);
319  desc.add<double>("Variable_PhiMinus_xmax", M_PI);
320  desc.add<int>("Variable_PhiMinus_nbins", 16);
321 
322  desc.add<double>("Variable_PhiPlus_xmin", -M_PI);
323  desc.add<double>("Variable_PhiPlus_xmax", M_PI);
324  desc.add<int>("Variable_PhiPlus_nbins", 16);
325 
326  desc.add<double>("Variable_PairPt_xmin", 0.);
327  desc.add<double>("Variable_PairPt_xmax", 100.);
328  desc.add<int>("Variable_PairPt_nbins", 100);
329 
330  descriptions.addWithDefaultLabel(desc);
331 }
332 
333 //define this as a plug-in
static const std::string kSharedResource
Definition: TFileService.h:76
int variables_bins_number_[varNumber_]
void addWithDefaultLabel(ParameterSetDescription const &psetDescription)
DiMuonValidation(const edm::ParameterSet &pset)
static constexpr double mu_mass2_
TH2D * th2d_mass_variables_[varNumber_]
double variable_CosThetaCS_xmax_
std::vector< Track > TrackCollection
collection of Tracks
Definition: TrackFwd.h:14
double variables_min_[varNumber_]
double variable_CosThetaCS_xmin_
int iEvent
Definition: GenABIO.cc:224
void swap(edm::DataFrameContainer &lhs, edm::DataFrameContainer &rhs)
T sqrt(T t)
Definition: SSEVec.h:19
std::string variables_name_[varNumber_]
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
#define M_PI
void beginJob() override
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
DecomposeProduct< arg, typename Div::arg > D
Definition: Factorize.h:141
double variables_max_[varNumber_]
HLT enums.
void analyze(const edm::Event &, const edm::EventSetup &) override
static constexpr int varNumber_
reco::Particle::LorentzVector LV
edm::EDGetTokenT< reco::TrackCollection > theTrackCollectionToken_
~DiMuonValidation() override=default
math::XYZTLorentzVector LorentzVector
Lorentz vector.
Definition: Particle.h:21
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29