ATLAS Offline Software
DiMuMon.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #include <sstream>
6 #include "GaudiKernel/PhysicalConstants.h"
8 
9 #include "DiMuMon.h"
10 
12 
13 #include <cmath>
14 
15 #include "TF1.h"
16 #include "TCanvas.h"
17 #include "TStyle.h"
18 
19 #include "RooRealVar.h"
20 #include "RooDataHist.h"
21 #include "RooPlot.h"
22 #include "RooBreitWigner.h"
23 #include "RooCBShape.h"
24 #include "RooFFTConvPdf.h"
25 #include "RooGlobalFunc.h"
26 #include "RooArgList.h"
27 #include "RooFitResult.h"
28 
29 
30 //using namespace Analysis;
31 //using namespace Rec;
32 
33 DiMuMon::DiMuMon( const std::string & type, const std::string & name, const IInterface* parent )
35  , m_triggerChainName("NoTrig")
36 {
37  declareProperty( "resonName", m_resonName = "Zmumu" );
38  declareProperty( "triggerChainName", m_triggerChainName = "NoTrig" );
39  declareProperty( "setDebug", m_setDebug = false );
40  declareProperty( "minInvmass", m_minInvmass = 60.);
41  declareProperty( "maxInvmass", m_maxInvmass = 120.);
42  declareProperty( "regions", m_regions );
43  declareProperty( "nVarBins", m_nVarBins = 10);
44  declareProperty( "nPtBins", m_nPtBins = 100);
45  declareProperty( "nEtaBins", m_nEtaBins = 32);
46  declareProperty( "nPhiBins", m_nPhiBins = 36);
47  declareProperty( "nMassBins", m_nMassBins = 60);
48  declareProperty( "varsVSmean", m_varsVSmean);
49  declareProperty( "varsVSwidth", m_varsVSwidth);
50  declareProperty( "varsDistr", m_varsDistr);
51  declareProperty( "doSaveFits", m_doSaveFits = false );
52  declareProperty( "doFits", m_doFits = false );
53  declareProperty( "ptCut", m_ptCut = 0.);//GeV
54  declareProperty( "etaCut", m_etaCut = 2.5);
55 }
56 
57 
59 {
60 }
61 
63 
65  ATH_CHECK( m_muonCollection.initialize() );
66 
67  if (m_regions.empty()) {
68  m_regions.emplace_back("All");
69  }
70  // m_variables[] = {"eta","etaAll","etaPos","etaNeg","phi","phiAll","phiPos","phiNeg","pt","ptAll","ptPos","ptNeg","etaDiff","etaSumm","phiDiff","phiSumm","crtDiff"};
71  m_varLabels["eta"] = "Dimuon #eta";
72  m_varLabels["etaAll"] = "#eta (all #mu)";
73  m_varLabels["etaPos"] = "#eta(+)";
74  m_varLabels["etaNeg"] = "#eta(-)";
75 
76  m_varLabels["phi"] = "Dimuon #phi";
77  m_varLabels["phiAll"] = "#phi (all #mu)";
78  m_varLabels["phiPos"] = "#phi(+)";
79  m_varLabels["phiNeg"] = "#phi(-)";
80 
81  m_varLabels["pt"] = "Dimuon p_{T}";
82  m_varLabels["ptAll"] = "p_{T} (all #mu)";
83  m_varLabels["ptPos"] = "p_{T}(+) [GeV]";
84  m_varLabels["ptNeg"] = "p_{T}(-) [GeV]";
85 
86  m_varLabels["etaDiff"] = "#eta(+)-#eta(-)";
87  m_varLabels["etaSumm"] = "#eta(+)+#eta(-)";
88  m_varLabels["phiDiff"] = "#phi(+)-#phi(-)";
89  m_varLabels["phiSumm"] = "#phi(+)+#phi(-)";
90  m_varLabels["crtDiff"] = "1/p_{T}(+) - 1/p_{T}(-) [GeV^{-1}]";
91 
92  //resonance independent
93  // for eta these are filled as the histograms are declared due to the dependence between region and eta
94  m_varRanges["phi"] = std::make_pair(-Gaudi::Units::pi,Gaudi::Units::pi);
95  m_varRanges["phiAll"] = std::make_pair(-Gaudi::Units::pi,Gaudi::Units::pi);
96  m_varRanges["phiPos"] = std::make_pair(-Gaudi::Units::pi,Gaudi::Units::pi);
97  m_varRanges["phiNeg"] = std::make_pair(-Gaudi::Units::pi,Gaudi::Units::pi);
98  m_varRanges["etaSumm"] = std::make_pair(-5.,5.);
99 
100  //resonance dependent
101  double ptMax = 0.;
102  if (m_resonName=="Jpsi" || m_resonName=="Upsi"){
103  m_varRanges["eta"] = std::make_pair(-3.,3.);
104  m_varRanges["phiDiff"] = std::make_pair(0.,1.);
105  m_varRanges["etaDiff"] = std::make_pair(-1.,1.);
106  m_varRanges["crtDiff"] = std::make_pair(-0.15,0.15);
107  m_varRanges["phiSumm"] = std::make_pair(-6.,6.);
108  if (m_ptCut == 0.) m_ptCut = 4.0;
109  ptMax = 18.;
110  } else if (m_resonName=="Zmumu") {
111  m_varRanges["eta"] = std::make_pair(-5.,5.);
112  m_varRanges["phiDiff"] = std::make_pair(0.,Gaudi::Units::pi);
113  m_varRanges["etaDiff"] = std::make_pair(-3.,3.);
114  m_varRanges["crtDiff"] = std::make_pair(-0.03,0.03);
115  m_varRanges["phiSumm"] = std::make_pair(-3.5,3.5);
116  if (m_ptCut == 0.) m_ptCut = 20.0;
117  ptMax = 100.;
118  }
119  m_varRanges["pt"] = std::make_pair(0.,ptMax);
120  m_varRanges["ptAll"] = std::make_pair(m_ptCut,ptMax);
121  m_varRanges["ptPos"] = std::make_pair(m_ptCut,ptMax);
122  m_varRanges["ptNeg"] = std::make_pair(m_ptCut,ptMax);
123 
124 
125  m_coneSize = 0.4;
126  m_isolationCut = 0.2;
127 
128  return StatusCode::SUCCESS;
129 
130 }
131 
133 {
134 
135 
136  MonGroup dimuMonObj_shift( this, "DiMuMon/"+m_resonName+"/" + m_triggerChainName, run );
137  MonGroup dimuMonObj_expert( this, "DiMuMon/"+m_resonName+"/" + m_triggerChainName + "_detail", run );
138 
140  // book histograms that are only made in the online environment...
141  }
142 
144  // book histograms that are only relevant for cosmics data...
145  }
146 
147 
148 
149 
150  if( newRunFlag() ) {
151 
152 
153  m_chi2 = new TH1F("chi2","chi2",100,0.,10.);
154  RegisterHisto(dimuMonObj_expert,m_chi2);
155 
156  const int nCuts = 8;
157  m_stat = new TH1F("statistics","Statistics",nCuts, 0., double(nCuts));
158  TString binLabels[] ={"MuonIdTrk","TrkIsoPt40/Pt<.2","CombMuon","MCPgoodTrk","MuPtOK","eta<2.5","OppChargePair","InvMassOK"};
159  for(int ilabel=0;ilabel<nCuts;ilabel++){
160  m_stat->GetXaxis()->SetBinLabel(ilabel+1,binLabels[ilabel]);
161  }
162  RegisterHisto(dimuMonObj_expert,m_stat);
163 
164  //create histos for each region
165  for (const std::string& reg : m_regions) {
166 
167  //declare region-dependent eta range
168  if (reg == "BB") {
169  m_varRanges["etaAll"] = std::make_pair(-1.05,1.05);
170  m_varRanges["etaPos"] = std::make_pair(-1.05,1.05);
171  m_varRanges["etaNeg"] = std::make_pair(-1.05,1.05);
172  } else if (reg == "EAEA") {
173  m_varRanges["etaAll"] = std::make_pair(1.05,2.5);
174  m_varRanges["etaPos"] = std::make_pair(1.05,2.5);
175  m_varRanges["etaNeg"] = std::make_pair(1.05,2.5);
176  } else if (reg == "ECEC") {
177  m_varRanges["etaAll"] = std::make_pair(-2.5,-1.05);
178  m_varRanges["etaPos"] = std::make_pair(-2.5,-1.05);
179  m_varRanges["etaNeg"] = std::make_pair(-2.5,-1.05);
180  } else {
181  m_varRanges["etaAll"] = std::make_pair(-2.5,2.5);
182  m_varRanges["etaPos"] = std::make_pair(-2.5,2.5);
183  m_varRanges["etaNeg"] = std::make_pair(-2.5,2.5);
184  }
185 
186  //mass plots
187  TString hname = m_resonName + "_invmass_" + reg;
188  TString htitle = hname + "; Invmass[GeV/c^{2}]";
190  RegisterHisto(dimuMonObj_shift,m_invmass[reg]);
191 
192  //for each var vs mass plot
193  for (const std::string& varM : m_varsVSmean) {
194  hname = m_resonName + "_2DinvmassVS" + varM + "_" + reg;
195  htitle = hname + ";" + m_varLabels[varM] + ";Invmass[GeV/c^{2}]";
197  RegisterHisto(dimuMonObj_shift,m_2DinvmassVSx[reg][varM]);
198 
199  hname = m_resonName + "_invmassVS" + varM + "_" + reg;
200  htitle = hname + ";" + m_varLabels[varM] + ";Invmass[GeV/c^{2}]";
201  m_invmassVSx[reg][varM] = new TH1F(hname, htitle, m_nVarBins,m_varRanges[varM].first, m_varRanges[varM].second);
202  RegisterHisto(dimuMonObj_shift,m_invmassVSx[reg][varM]);
203  }
204 
205  //for each var vs width plot
206  for (const std::string& varW : m_varsVSwidth) {
207  //book the corresponding 2D histo if it hasn't already been booked
208  if (m_2DinvmassVSx[reg].find(varW)==m_2DinvmassVSx[reg].end()){
209  hname = m_resonName + "_2DinvmassVS" + varW + "_" + reg;
210  htitle = hname + ";" + m_varLabels[varW] + ";Invmass[GeV/c^{2}]";
212  RegisterHisto(dimuMonObj_shift,m_2DinvmassVSx[reg][varW]);
213  }
214 
215  hname = m_resonName + "_widthVS" + varW + "_" + reg;
216  htitle = hname + ";" + m_varLabels[varW] + ";Width[GeV/c^{2}]";
217  m_widthVSx[reg][varW] = new TH1F(hname,htitle ,m_nVarBins,m_varRanges[varW].first, m_varRanges[varW].second);
218  RegisterHisto(dimuMonObj_expert,m_widthVSx[reg][varW]);
219 
220  }
221 
222  //for the each variable's distribution
223  for (const std::string& varD: m_varsDistr) {
224  hname = m_resonName + "_" + varD + "_" + reg;
225  htitle = hname + ";" + m_varLabels[varD];
226  if (varD == "eta" || varD == "etaAll" || varD == "etaPos" || varD == "etaNeg" ){
227  m_xDistr[reg][varD] = new TH1F(hname, htitle, m_nEtaBins, m_varRanges[varD].first, m_varRanges[varD].second);
228  } else if (varD == "pt" || varD == "ptAll" || varD == "ptPos" || varD == "ptNeg" ){
229  m_xDistr[reg][varD] = new TH1F(hname, htitle, m_nPtBins, m_varRanges[varD].first, m_varRanges[varD].second);
230  } else if (varD == "phi" || varD == "phiAll" || varD == "phiPos" || varD == "phiNeg" ){
231  m_xDistr[reg][varD] = new TH1F(hname, htitle, m_nPhiBins, m_varRanges[varD].first, m_varRanges[varD].second);
232  } else {
233  m_xDistr[reg][varD] = new TH1F(hname, htitle, m_nVarBins, m_varRanges[varD].first, m_varRanges[varD].second);
234  }
235  RegisterHisto(dimuMonObj_expert,m_xDistr[reg][varD]);
236  }
237  }
238  }
239 
240  return StatusCode::SUCCESS;
241 }
242 
243 
245 {
246 
247  const double muonMass = 105.66*Gaudi::Units::MeV;
248  //retrieve all muons
250  if(!muons.isValid()){
251  ATH_MSG_WARNING("Could not retrieve muon container");
252  return StatusCode::FAILURE;
253  } else ATH_MSG_DEBUG("Muon container successfully retrieved.");
254 
255  //make a new container
257 
258  //pick out the good muon tracks and store in the new container
259  for(const auto* muon : *muons ) {
260 
261  const xAOD::TrackParticle *idTrk = muon->trackParticle(xAOD::Muon::InnerDetectorTrackParticle);
262  if (!idTrk) continue;
263  m_stat->Fill("MuonIdTrk",1);
264  double idTrkPt(0),ptSum(0);
265  float iso_pt40(0);
266  if( !muon->isolation(iso_pt40, xAOD::Iso::ptcone40) ) {
267  return StatusCode::FAILURE;
268  }
269  else {
270  idTrkPt = idTrk->pt();
271  ptSum = xAOD::Iso::ptcone40;
272  }
273  if(m_resonName=="Zmumu" && ptSum/idTrkPt > m_isolationCut) continue;
274  m_stat->Fill("TrkIsoPt40/Pt<.2",1);
275  if (!muon->trackParticle(xAOD::Muon::CombinedTrackParticle)) continue;
276  m_stat->Fill("CombMuon",1);
277 
278  if (!trackQuality(idTrk)) continue;
279  m_stat->Fill("MCPgoodTrk",1);
280 
281  if (idTrkPt < m_ptCut*1000) continue;
282  m_stat->Fill("MuPtOK",1);
283 
284  double idTrkEta = idTrk->eta();
285  if (std::abs(idTrkEta)>2.5) continue;
286  m_stat->Fill("eta<2.5",1);
287 
288  goodMuons.push_back(muon);
289  }
290 
291  //pair up the tracks of the good muons and fill histograms
292  int nMuons = goodMuons.size();
293 
294  if (nMuons>1){
295  xAOD::MuonContainer::const_iterator mu1 = goodMuons.begin();
296  xAOD::MuonContainer::const_iterator muEnd = goodMuons.end();
297  for (; mu1!=muEnd;++mu1){
298  const xAOD::TrackParticle *id1 = (*mu1)->trackParticle(xAOD::Muon::InnerDetectorTrackParticle);
300  for (; mu2!=muEnd; ++mu2){
301  const xAOD::TrackParticle *id2 = (*mu2)->trackParticle(xAOD::Muon::InnerDetectorTrackParticle);
302  //consider only opposite sign muons
303  double q1 = id1->charge();
304  double q2 = id2->charge();
305  if (q1*q2>0) continue;
306  m_stat->Fill("OppChargePair",1);
307  const xAOD::TrackParticle *idPos = id1;
308  const xAOD::TrackParticle *idNeg = id2;
309  if (q1<0){
310  idPos = id2;
311  idNeg = id1;
312  }
313 
314  //cut on the pair invariant mass
315  double invmass = getInvmass(id1,id2,muonMass);
316  if (invmass<m_minInvmass || invmass>m_maxInvmass) continue;
317  m_stat->Fill("InvMassOK",1);
318 
319  //compute variables
320  m_varValues["eta"] = getEta(idPos,idNeg);
321  double etaPos = idPos->eta();
322  m_varValues["etaPos"] = etaPos;
323  double etaNeg = idNeg->eta();
324  m_varValues["etaNeg"] = etaNeg;
325  m_varValues["phi"] = getPhi(idPos,idNeg);
326  double phiPos = idPos->phi();
327  m_varValues["phiPos"] = phiPos;
328  double phiNeg = idNeg->phi();
329  m_varValues["phiNeg"] = phiNeg;
330  m_varValues["pt"] = getPt(idPos,idNeg);
331  double ptPos = idPos->pt()/Gaudi::Units::GeV;
332  m_varValues["ptPos"] = ptPos;
333  double ptNeg = idNeg->pt()/Gaudi::Units::GeV;
334  m_varValues["ptNeg"] = ptNeg;
335 
336  m_varValues["crtDiff"] = getCrtDiff(idPos,idNeg);
337  m_varValues["etaDiff"] = etaPos - etaNeg;
338  double phiDiff = std::abs(phiPos - phiNeg);
339  if (phiDiff>Gaudi::Units::pi) phiDiff = 2*(Gaudi::Units::pi) - phiDiff;
340  m_varValues["phiDiff"] = phiDiff;
341  m_varValues["etaSumm"] = etaPos + etaNeg;
342  m_varValues["phiSumm"] = phiPos + phiNeg;
343 
344  //determine which region muons are in
345  std::string region = "";
346  if ((std::abs(etaPos)<1.05 || std::abs(etaPos)==1.05) && (std::abs(etaNeg)<1.05 || std::abs(etaNeg)==1.05)) region="BB";
347  else if ((etaPos>1.05 && etaPos<2.5) && (etaNeg>1.05 && etaNeg<2.5)) region="EAEA";
348  else if ((etaPos<-1.05 && etaPos>-2.5) && (etaNeg<-1.05 && etaNeg>-2.5)) region="ECEC";
349 
350  //do we care about hese muons?
351  bool fillAll = m_invmass.find("All")!=m_invmass.end();
352  bool fillReg = false;
353  if (region!=""){
354  if (m_invmass.find(region)!=m_invmass.end()) fillReg=true;
355  }
356  if (!fillAll && !fillReg) continue;
357 
358 
359  //fill invmass histos
360  if (fillAll) m_invmass["All"]->Fill(invmass);
361  if (fillReg) m_invmass[region]->Fill(invmass);
362 
363  //fill 2D histos
364  //first, retrieve the overlap of the width and mean variable lists which by construction is the list of vars in the 2D histos map
365  std::map<std::string, TH2F*>* allVars = nullptr;
366  if (fillAll){
367  allVars = &m_2DinvmassVSx["All"];
368  } else { //fillReg=true
369  allVars = &m_2DinvmassVSx[region];
370  }
371 
372  for (const std::pair<const std::string, TH2F*>& p : *allVars) {
373  if (p.first!="etaAll" && p.first!="phiAll" && p.first!="ptAll"){
374  if (fillAll) m_2DinvmassVSx["All"][p.first]->Fill(m_varValues[p.first],invmass);
375  if (fillReg) m_2DinvmassVSx[region][p.first]->Fill(m_varValues[p.first],invmass);
376  } else if (p.first=="etaAll"){
377  if (fillAll){
378  m_2DinvmassVSx["All"]["etaAll"]->Fill(etaPos,invmass);
379  m_2DinvmassVSx["All"]["etaAll"]->Fill(etaNeg,invmass);
380  }
381  if (fillReg){
382  m_2DinvmassVSx[region]["etaAll"]->Fill(etaPos,invmass);
383  m_2DinvmassVSx[region]["etaAll"]->Fill(etaNeg,invmass);
384  }
385  } else if (p.first=="phiAll"){
386  if (fillAll){
387  m_2DinvmassVSx["All"]["phiAll"]->Fill(phiPos,invmass);
388  m_2DinvmassVSx["All"]["phiAll"]->Fill(phiNeg,invmass);
389  }
390  if (fillReg){
391  m_2DinvmassVSx[region]["phiAll"]->Fill(phiPos,invmass);
392  m_2DinvmassVSx[region]["phiAll"]->Fill(phiNeg,invmass);
393  }
394  } else if (p.first=="ptAll"){
395  if (fillAll){
396  m_2DinvmassVSx["All"]["ptAll"]->Fill(ptPos,invmass);
397  m_2DinvmassVSx["All"]["ptAll"]->Fill(ptNeg,invmass);
398  }
399  if (fillReg){
400  m_2DinvmassVSx[region]["ptAll"]->Fill(ptPos,invmass);
401  m_2DinvmassVSx[region]["ptAll"]->Fill(ptNeg,invmass);
402  }
403  }
404  }
405 
406  //fill var distributions
407  //here we already know the list of variables so no need for gymnastics
408  for (const std::string& varD : m_varsDistr) {
409  if (varD!="etaAll" && varD!="phiAll" && varD!="ptAll"){
410  if (fillAll) m_xDistr["All"][varD]->Fill(m_varValues[varD]);
411  if (fillReg) m_xDistr[region][varD]->Fill(m_varValues[varD]);
412  } else if (varD=="etaAll"){
413  if (fillAll){
414  m_xDistr["All"]["etaAll"]->Fill(etaPos);
415  m_xDistr["All"]["etaAll"]->Fill(etaNeg);
416  }
417  if (fillReg){
418  m_xDistr[region]["etaAll"]->Fill(etaPos);
419  m_xDistr[region]["etaAll"]->Fill(etaNeg);
420  }
421  } else if (varD=="phiAll"){
422  if (fillAll){
423  m_xDistr["All"]["phiAll"]->Fill(phiPos);
424  m_xDistr["All"]["phiAll"]->Fill(phiNeg);
425  }
426  if (fillReg){
427  m_xDistr[region]["phiAll"]->Fill(phiPos);
428  m_xDistr[region]["phiAll"]->Fill(phiNeg);
429  }
430  } else if (varD=="ptAll"){
431  if (fillAll){
432  m_xDistr["All"]["ptAll"]->Fill(ptPos);
433  m_xDistr["All"]["ptAll"]->Fill(ptNeg);
434  }
435  if (fillReg){
436  m_xDistr[region]["ptAll"]->Fill(ptPos);
437  m_xDistr[region]["ptAll"]->Fill(ptNeg);
438  }
439  }
440  }
441 
442 
443  }// mu2 loop
444  }// mu1 loop
445  }// do we have more than 1 good muon?
446  return StatusCode::SUCCESS;
447 }
448 
449 
451 {
452 
453 
454  if(endOfRunFlag() && m_doFits) {
455  for (const std::string& reg : m_regions) {
456  for (const std::pair<const std::string, TH2F*>& p : m_2DinvmassVSx[reg]) {
457  int mode = 999.; //0 = both, 1 = mean only, 2 = wdth only
458  std::vector<TH1F*> hout;//1D histograms for fit results
459  if (m_invmassVSx[reg].find(p.first)!=m_invmassVSx[reg].end()){
460  if (m_widthVSx[reg].find(p.first)!=m_widthVSx[reg].end()) {
461  mode =0;
462  hout.push_back(m_invmassVSx[reg][p.first]);
463  hout.push_back(m_widthVSx[reg][p.first]);
464  } else {
465  mode = 1;
466  hout.push_back(m_invmassVSx[reg][p.first]);
467  }
468  } else {
469  mode = 2;
470  hout.push_back(m_widthVSx[reg][p.first]);
471  }
472  iterativeGausFit(m_2DinvmassVSx[reg][p.first],hout,mode);
473  } //variables
474  } //regions
475  } //isEndOfRun
476  return StatusCode::SUCCESS;
477 }
478 
479 
480 void DiMuMon::iterativeGausFit (TH2F* hin, const std::vector<TH1F*>& hout, int mode){
481  // a canvas may be needed when implmenting this into the post-processing file
482  TString hname = hin->GetName();
483  TString psName = hname + m_triggerChainName + ".ps";
484  TCanvas* ctemp = new TCanvas("ctemp","ctemp",500,500);
485  if (m_doSaveFits) ctemp->Print(psName+"[");
486  hin->SetMarkerSize(1.2);
487  hin->Draw();
488  if (m_doSaveFits) ctemp->Print(psName);
489  int nbins=hin->GetNbinsX();
490  for (int i=0; i<nbins;i++){
491  std::ostringstream o; o<<i;
492  TString projName = hname + o.str();
493  TH1D* htemp = (TH1D*) (hin->ProjectionY(projName,i+1,i+1)->Clone());
494  htemp->SetTitle(projName);
495  htemp->Sumw2();
496  if (htemp->GetEntries()>50){
497  double mean = 999., meanErr = 999., sigma = 999., sigmaErr = 999., chi2=0;
498  if (m_resonName == "Jpsi" || m_resonName == "Upsi"){
499  mean = htemp->GetMean();
500  sigma= htemp->GetRMS();
501  TF1* fn = new TF1("fn","gaus",mean-2*sigma,mean+2*sigma);
502  fn->SetParameters(float(htemp->GetEntries())/10.,mean,sigma);
503  htemp->Fit("fn","RMLQN");
504  mean = fn->GetParameter(1);
505  sigma= fn->GetParameter(2);
506  fn->SetRange(mean-1.2*sigma,mean+1.2*sigma);
507  fn->SetParameters(float(htemp->GetEntries())/10.,mean,sigma);
508  //gStyle->SetOptStat(1); // not thread-safe
509  if (m_doSaveFits) {
510  htemp->Fit("fn","RML");
511  ctemp->Print(psName);
512  } else htemp->Fit("fn","RMLQN");
513  double frange = 2.4*sigma;
514  double hrange = htemp->GetXaxis()->GetXmax()-htemp->GetXaxis()->GetXmin();
515  double ndf = frange/hrange*(htemp->GetNbinsX()) - 3;//substract the fit parameters
516  chi2 = (fn->GetChisquare())/ndf;
517  //fill results
518  mean = fn->GetParameter(1);
519  meanErr = fn->GetParError(1);
520  sigma = fn->GetParameter(2);
521  sigmaErr = fn->GetParError(2);
522  delete fn;
523  } else {
524  //fit Z peak with a convolution of BreitWigner and Crystal Ball fns, fit by Louise, implementation by Jike taken from IDPerfMon
525  RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
526  RooRealVar m("mass", "dimuon invariant mass", 91.2, 71., 111., "GeV");
527  RooDataHist *data = nullptr;
528  data = new RooDataHist("data", "data", m, htemp);
529  RooRealVar bwm0("bw_#mu", "bw_#mu", 91.2, 85.2, 97.2) ;
530  RooRealVar bwsg("bw_#sigma","bw_#sigma", 2.4952) ;
531  RooBreitWigner bw("bw","bw", m, bwm0, bwsg);
532 
533  RooRealVar cbm0("cb_#mu", "cb_#mu", 0 ) ;
534  RooRealVar cbsg("cb_#sigma", "cb_#sigma", 3., 1., 10.) ;
535  RooRealVar cbal("cb_#alpha", "cb_#alpha", 2.0) ;
536  RooRealVar cbn( "cb_n", "cb_n", 1., 0.05, 3.) ;
537  RooCBShape cb( "cb", "cb", m, cbm0, cbsg, cbal, cbn);
538 
539  m.setBins(5000);
540  RooFFTConvPdf bxc("bxc", "BW (X) CB", m, bw, cb) ;
541  bxc.fitTo(*data, RooFit::PrintLevel(-1), RooFit::PrintEvalErrors(-1));
542  RooPlot* frame = m.frame();
543  data->plotOn(frame, RooFit::MarkerSize(0.9));
544  bxc.paramOn(frame, RooFit::Format("NELU", RooFit::AutoPrecision(2)), RooFit::Layout(0.1,0.4,0.9));
545  bxc.plotOn (frame, RooFit::LineColor(kBlue));
546  if (m_doSaveFits) {
547  frame->Draw();
548  ctemp->Print(psName);
549  }
550  mean = bwm0.getVal();
551  meanErr = bwm0.getError();
552  sigma = cbsg.getVal();
553  sigmaErr = cbsg.getError();
554  chi2 = frame->chiSquare();
555  delete data;
556  }
557  //fill results
558  m_chi2->Fill(chi2);
559  if (mode == 0){//plot both invmass and width vs variable
560  hout.at(0)->SetBinContent(i+1,mean);
561  hout.at(0)->SetBinError(i+1,meanErr);
562  hout.at(1)->SetBinContent(i+1,sigma);
563  hout.at(1)->SetBinError(i+1,sigmaErr);
564  } else if (mode==1){//plot only invmass vs variable
565  hout.at(0)->SetBinContent(i+1,mean);
566  hout.at(0)->SetBinError(i+1,meanErr);
567  } else if (mode==2){//plot only width vs variable
568  hout.at(0)->SetBinContent(i+1,sigma);
569  hout.at(0)->SetBinError(i+1,sigmaErr);
570  }
571 
572  }// more than 50 events
573 
574  delete htemp;
575  }
576  hout.at(0)->Sumw2();
577  if (hout.size()==2) hout.at(1)->Sumw2();
578  if (m_doSaveFits) ctemp->Print(psName+"]");
579  delete ctemp;
580 }
581 
582 template<class T>
584 
585  StatusCode sc = mon.regHist(histo);
586  if (sc.isFailure() ) {
587  ATH_MSG_WARNING( "Cannot book histogram:" );
588  }
589 }
590 
592  int countPass = 0;
593  // https://twiki.cern.ch/twiki/bin/view/AtlasProtected/MCPAnalysisGuidelinesEPS2011
594  if (idTrk) {
595  uint8_t dummy(-1);
596  bool eBLhits = idTrk->summaryValue( dummy, xAOD::expectInnermostPixelLayerHit )? dummy :false;
597  int nBLhits = idTrk->summaryValue( dummy, xAOD::numberOfInnermostPixelLayerHits )? dummy :-1;
598  int nhitsPIX = idTrk->summaryValue( dummy, xAOD::numberOfPixelHits )? dummy :-1;
599  int nPIXDS = idTrk->summaryValue( dummy, xAOD::numberOfPixelDeadSensors )? dummy :-1;
600  int nhitsSCT = idTrk->summaryValue( dummy, xAOD::numberOfSCTHits )? dummy :-1;
601  int nSCTDS = idTrk->summaryValue( dummy, xAOD::numberOfSCTDeadSensors )? dummy :-1;
602  int nPIXH = idTrk->summaryValue( dummy, xAOD::numberOfPixelHoles )? dummy :-1;
603  int nSCTH = idTrk->summaryValue( dummy, xAOD::numberOfSCTHoles )? dummy :-1;
604 
605  if (nBLhits>0 || !(eBLhits)) countPass+=1;
606  if ((nhitsPIX+nPIXDS)>1) countPass+=1;
607  if ((nhitsSCT+nSCTDS)>5) countPass+=1;
608  if ((nPIXH+nSCTH)<3) countPass+=1;
609  int nTRTout = idTrk->summaryValue( dummy, xAOD::numberOfTRTOutliers )? dummy :-1;
610  int nTRThits = idTrk->summaryValue( dummy, xAOD::numberOfTRTHits )? dummy :-1;
611  int n = nTRTout + nTRThits;
612  if (std::abs(idTrk->eta())<1.9){
613  if (n>5 && nTRTout<(0.9*n)) countPass+=1;
614  } else {
615  if (n>5){
616  if (nTRTout<(0.9*n)) countPass+=1;
617  } else countPass+=1;
618  }
619  }
620  else{
621  ATH_MSG_WARNING("Trying to check trackquality but no xAOD::TrackParticle available!");
622  }
623  return countPass==5;
624 }
625 
626 //methods from IDPerformanceMoniroting by Weina
627 double DiMuMon::getInvmass(const xAOD::TrackParticle* id1, const xAOD::TrackParticle* id2, double Mass) const {
628  TLorentzVector particle1,particle2,v;
629  particle1.SetPtEtaPhiE(id1->pt(),id1->eta(),id1->phi(),sqrt(pow(Mass,2)+pow(id1->p4().Px(),2)+pow(id1->p4().Py(),2)+pow(id1->p4().Pz(),2)));
630  particle2.SetPtEtaPhiE(id2->pt(),id2->eta(),id2->phi(),sqrt(pow(Mass,2)+pow(id2->p4().Px(),2)+pow(id2->p4().Py(),2)+pow(id2->p4().Pz(),2)));
631  v=particle1+particle2;
632  double invmass = v.Mag()/Gaudi::Units::GeV;
633  return invmass;
634 }
635 
636 double DiMuMon::getPt(const xAOD::TrackParticle* id1, const xAOD::TrackParticle* id2 ) const {
637  double transmom;
638  double px = id1->p4().Px()+id2->p4().Px();
639  double py = id1->p4().Py()+id2->p4().Py();
640  transmom=sqrt(px*px+py*py);
641  return transmom/Gaudi::Units::GeV; //Gev
642 }
643 
644 double DiMuMon::getEta(const xAOD::TrackParticle* id1, const xAOD::TrackParticle* id2 ) const {
645  double eta;
646  double px=id1->p4().Px()+id2->p4().Px();
647  double py = id1->p4().Py()+id2->p4().Py();
648  double pz = id1->p4().Pz()+id2->p4().Pz();
649  double p=sqrt(px*px+py*py+pz*pz);
650  eta=0.5*log((p+pz)/(p-pz));
651  return eta;
652 }
653 
654 double DiMuMon::getPhi(const xAOD::TrackParticle* id1, const xAOD::TrackParticle* id2 ) const {
655  double px=id1->p4().Px()+id2->p4().Px();
656  double py = id1->p4().Py()+id2->p4().Py();
657  double p=sqrt(px*px+py*py);
658  double phi=acos(px/p);
659  if (py<0.) phi=-1.*phi;
660  return phi;
661 }
662 //curvature difference unit 1/GeV
664  double qoverpt1=id1->charge()/id1->pt();
665  double qoverpt2=id2->charge()/id2->pt();
666  double asym;
667  asym=std::abs(qoverpt1)-std::abs(qoverpt2);
668  return asym*1000;
669 }
670 
671 
covarianceTool.ndf
ndf
Definition: covarianceTool.py:678
xAOD::TrackParticle_v1::pt
virtual double pt() const override final
The transverse momentum ( ) of the particle.
Definition: TrackParticle_v1.cxx:73
xAOD::numberOfPixelHoles
@ numberOfPixelHoles
number of pixel layers on track with absence of hits [unit8_t].
Definition: TrackingPrimitives.h:261
python.SystemOfUnits.second
int second
Definition: SystemOfUnits.py:120
DiMuMon::procHistograms
virtual StatusCode procHistograms()
An inheriting class should either override this function or finalHists().
Definition: DiMuMon.cxx:450
xAOD::muon
@ muon
Definition: TrackingPrimitives.h:195
data
char data[hepevt_bytes_allocation_ATLAS]
Definition: HepEvt.cxx:11
DiMuMon::m_varRanges
std::map< std::string, std::pair< double, double > > m_varRanges
Definition: DiMuMon.h:77
DataModel_detail::const_iterator
Const iterator class for DataVector/DataList.
Definition: DVLIterator.h:82
DiMuMon::fillHistograms
virtual StatusCode fillHistograms()
An inheriting class should either override this function or fillHists().
Definition: DiMuMon.cxx:244
pdg_comparison.sigma
sigma
Definition: pdg_comparison.py:324
mean
void mean(std::vector< double > &bins, std::vector< double > &values, const std::vector< std::string > &files, const std::string &histname, const std::string &tplotname, const std::string &label="")
Definition: dependence.cxx:254
test_pyathena.px
px
Definition: test_pyathena.py:18
python.SystemOfUnits.m
int m
Definition: SystemOfUnits.py:91
python.PerfMonSerializer.p
def p
Definition: PerfMonSerializer.py:743
phi
Scalar phi() const
phi method
Definition: AmgMatrixBasePlugin.h:64
python.Constants.FATAL
int FATAL
Definition: Control/AthenaCommon/python/Constants.py:19
DiMuMon::getEta
double getEta(const xAOD::TrackParticle *id1, const xAOD::TrackParticle *id2) const
Definition: DiMuMon.cxx:644
xAOD::uint8_t
uint8_t
Definition: Muon_v1.cxx:575
find
std::string find(const std::string &s)
return a remapped string
Definition: hcg.cxx:135
DiMuMon::m_varsVSmean
std::vector< std::string > m_varsVSmean
Definition: DiMuMon.h:74
DiMuMon::m_varsVSwidth
std::vector< std::string > m_varsVSwidth
Definition: DiMuMon.h:75
ManagedMonitorToolBase
Provides functionality for users to implement and save histograms, ntuples, and summary data,...
Definition: ManagedMonitorToolBase.h:74
SG::VIEW_ELEMENTS
@ VIEW_ELEMENTS
this data object is a view, it does not own its elmts
Definition: OwnershipPolicy.h:18
TH2F
Definition: rootspy.cxx:420
ConstDataVector::end
iterator end() noexcept
Return an iterator pointing past the end of the collection.
xAOD::TrackParticle_v1::charge
float charge() const
Returns the charge.
Definition: TrackParticle_v1.cxx:150
dqt_zlumi_pandas.hname
string hname
Definition: dqt_zlumi_pandas.py:272
eta
Scalar eta() const
pseudorapidity method
Definition: AmgMatrixBasePlugin.h:79
ConstDataVector.h
DataVector adapter that acts like it holds const pointers.
SG::ReadHandle
Definition: StoreGate/StoreGate/ReadHandle.h:70
AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T > &t)
Definition: AthCommonDataStore.h:145
xAOD::TrackParticle_v1::eta
virtual double eta() const override final
The pseudorapidity ( ) of the particle.
Definition: TrackParticle_v1.cxx:77
StateLessPT_NewConfig.Format
Format
Definition: StateLessPT_NewConfig.py:146
DiMuMon::m_chi2
TH1F * m_chi2
Definition: DiMuMon.h:81
xAOD::TrackParticle_v1::summaryValue
bool summaryValue(uint8_t &value, const SummaryType &information) const
Accessor for TrackSummary values.
Definition: TrackParticle_v1.cxx:736
AthenaMonManager::cosmics
@ cosmics
Definition: AthenaMonManager.h:58
conifer::pow
constexpr int pow(int x)
Definition: conifer.h:20
DiMuMon.h
python.SystemOfUnits.MeV
int MeV
Definition: SystemOfUnits.py:154
TH1D
Definition: rootspy.cxx:342
DiMuMon::getPt
double getPt(const xAOD::TrackParticle *id1, const xAOD::TrackParticle *id2) const
Definition: DiMuMon.cxx:636
DiMuMon::m_resonName
std::string m_resonName
Definition: DiMuMon.h:64
xAOD::numberOfPixelHits
@ numberOfPixelHits
these are the pixel hits, including the b-layer [unit8_t].
Definition: TrackingPrimitives.h:259
DiMuMon::m_isolationCut
double m_isolationCut
Definition: DiMuMon.h:58
xAOD::expectInnermostPixelLayerHit
@ expectInnermostPixelLayerHit
Do we expect a 0th-layer barrel hit for this track?
Definition: TrackingPrimitives.h:236
xAOD::numberOfTRTHits
@ numberOfTRTHits
number of TRT hits [unit8_t].
Definition: TrackingPrimitives.h:275
DiMuMon::m_nVarBins
int m_nVarBins
Definition: DiMuMon.h:52
DiMuMon::~DiMuMon
virtual ~DiMuMon()
Definition: DiMuMon.cxx:58
DiMuMon::m_nEtaBins
int m_nEtaBins
Definition: DiMuMon.h:54
python.TrigEgammaMonitorHelper.TH2F
def TH2F(name, title, nxbins, bins_par2, bins_par3, bins_par4, bins_par5=None, bins_par6=None, path='', **kwargs)
Definition: TrigEgammaMonitorHelper.py:45
python.DomainsRegistry.reg
reg
globals -----------------------------------------------------------------—
Definition: DomainsRegistry.py:343
pi
#define pi
Definition: TileMuonFitter.cxx:65
AthenaPoolTestRead.sc
sc
Definition: AthenaPoolTestRead.py:27
mergePhysValFiles.end
end
Definition: DataQuality/DataQualityUtils/scripts/mergePhysValFiles.py:93
DiMuMon::m_xDistr
std::map< std::string, std::map< std::string, TH1F * > > m_xDistr
Definition: DiMuMon.h:87
AthenaMonManager::environment
static Environment_t environment()
Returns the running environment of the monitoring application to help ManagedMonitorToolBase objects ...
Definition: AthenaMonManager.cxx:321
DiMuMon::m_invmass
std::map< std::string, TH1F * > m_invmass
Definition: DiMuMon.h:83
SCT_CalibAlgs::nbins
@ nbins
Definition: SCT_CalibNumbers.h:10
instance
std::map< std::string, double > instance
Definition: Run_To_Get_Tags.h:8
DiMuMon::m_etaCut
double m_etaCut
Definition: DiMuMon.h:62
xAOD::TrackParticle_v1::p4
virtual FourMom_t p4() const override final
The full 4-momentum of the particle.
Definition: TrackParticle_v1.cxx:129
DiMuMon::m_setDebug
bool m_setDebug
Definition: DiMuMon.h:70
ManagedMonitorToolBase::MonGroup
A container of information describing a monitoring object.
Definition: ManagedMonitorToolBase.h:138
DiMuMon::m_stat
TH1F * m_stat
Definition: DiMuMon.h:82
python.getCurrentFolderTag.fn
fn
Definition: getCurrentFolderTag.py:65
AthenaMonManager::dataType
static DataType_t dataType()
Returns the data type that the monitoring application is running over to help ManagedMonitorToolBase ...
Definition: AthenaMonManager.cxx:338
ManagedMonitorToolBase::initialize
virtual StatusCode initialize()
Definition: ManagedMonitorToolBase.cxx:669
id2
HWIdentifier id2
Definition: LArRodBlockPhysicsV0.cxx:564
lumiFormat.i
int i
Definition: lumiFormat.py:92
DiMuMon::m_2DinvmassVSx
std::map< std::string, std::map< std::string, TH2F * > > m_2DinvmassVSx
Definition: DiMuMon.h:84
DiMuMon::m_invmassVSx
std::map< std::string, std::map< std::string, TH1F * > > m_invmassVSx
Definition: DiMuMon.h:85
beamspotman.n
n
Definition: beamspotman.py:731
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
AthenaMonManager::online
@ online
Definition: AthenaMonManager.h:49
chi2
double chi2(TH1 *h0, TH1 *h1)
Definition: comparitor.cxx:522
Amg::pz
@ pz
Definition: GeoPrimitives.h:40
test_pyathena.parent
parent
Definition: test_pyathena.py:15
python.xAODType.dummy
dummy
Definition: xAODType.py:4
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
Preparation.mode
mode
Definition: Preparation.py:95
jet::CompScaleVar::Mass
@ Mass
Definition: UncertaintyEnum.h:98
run
Definition: run.py:1
xAOD::numberOfSCTHoles
@ numberOfSCTHoles
number of SCT holes [unit8_t].
Definition: TrackingPrimitives.h:270
DiMuMon::m_coneSize
double m_coneSize
Definition: DiMuMon.h:57
DiMuMon::m_minInvmass
double m_minInvmass
Definition: DiMuMon.h:59
DiMuMon::m_nPtBins
int m_nPtBins
Definition: DiMuMon.h:53
DiMuMon::m_muonCollection
SG::ReadHandleKey< xAOD::MuonContainer > m_muonCollection
Definition: DiMuMon.h:66
DiMuMon::m_widthVSx
std::map< std::string, std::map< std::string, TH1F * > > m_widthVSx
Definition: DiMuMon.h:86
DiMuMon::m_maxInvmass
double m_maxInvmass
Definition: DiMuMon.h:60
Amg::py
@ py
Definition: GeoPrimitives.h:39
DiMuMon::m_doFits
bool m_doFits
Definition: DiMuMon.h:68
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:195
DiMuMon::m_triggerChainName
std::string m_triggerChainName
Definition: DiMuMon.h:65
ConstDataVector::push_back
value_type push_back(value_type pElem)
Add an element to the end of the collection.
DiMuMon::getPhi
double getPhi(const xAOD::TrackParticle *id1, const xAOD::TrackParticle *id2) const
Definition: DiMuMon.cxx:654
DiMuMon::m_varValues
std::map< std::string, double > m_varValues
Definition: DiMuMon.h:78
DiMuMon::getCrtDiff
double getCrtDiff(const xAOD::TrackParticle *id1, const xAOD::TrackParticle *id2) const
Definition: DiMuMon.cxx:663
DiMuMon::m_doSaveFits
bool m_doSaveFits
Definition: DiMuMon.h:69
DiMuMon::bookHistograms
virtual StatusCode bookHistograms()
An inheriting class should either override this function or bookHists().
Definition: DiMuMon.cxx:132
xAOD::Iso::ptcone40
@ ptcone40
Definition: IsolationType.h:42
xAOD::numberOfTRTOutliers
@ numberOfTRTOutliers
number of TRT outliers [unit8_t].
Definition: TrackingPrimitives.h:276
python.PyAthena.v
v
Definition: PyAthena.py:157
DiMuMon::trackQuality
bool trackQuality(const xAOD::TrackParticle *idTrk)
Definition: DiMuMon.cxx:591
DiMuMon::m_regions
std::vector< std::string > m_regions
Definition: DiMuMon.h:73
DiMuMon::initialize
virtual StatusCode initialize()
Definition: DiMuMon.cxx:62
DiMuMon::iterativeGausFit
void iterativeGausFit(TH2F *hin, const std::vector< TH1F * > &hout, int mode)
Definition: DiMuMon.cxx:480
DiMuMon::DiMuMon
DiMuMon(const std::string &type, const std::string &name, const IInterface *parent)
Definition: DiMuMon.cxx:33
StateLessPT_NewConfig.Layout
Layout
Definition: StateLessPT_NewConfig.py:159
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
xAOD::numberOfSCTDeadSensors
@ numberOfSCTDeadSensors
number of dead SCT sensors crossed [unit8_t].
Definition: TrackingPrimitives.h:273
ConstDataVector
DataVector adapter that acts like it holds const pointers.
Definition: ConstDataVector.h:76
ManagedMonitorToolBase::endOfRunFlag
bool endOfRunFlag() const
Definition: ManagedMonitorToolBase.h:859
python.CaloScaleNoiseConfig.type
type
Definition: CaloScaleNoiseConfig.py:78
DeMoScan.first
bool first
Definition: DeMoScan.py:534
plotBeamSpotMon.mon
mon
Definition: plotBeamSpotMon.py:67
DiMuMon::RegisterHisto
void RegisterHisto(MonGroup &mon, T *histo)
Definition: DiMuMon.cxx:583
DiMuMon::m_nPhiBins
int m_nPhiBins
Definition: DiMuMon.h:55
python.CaloCondTools.log
log
Definition: CaloCondTools.py:20
xAOD::numberOfSCTHits
@ numberOfSCTHits
number of hits in SCT [unit8_t].
Definition: TrackingPrimitives.h:268
DiMuMon::m_ptCut
double m_ptCut
Definition: DiMuMon.h:61
xAOD::numberOfPixelDeadSensors
@ numberOfPixelDeadSensors
number of dead pixel sensors crossed [unit8_t].
Definition: TrackingPrimitives.h:266
python.TrigEgammaMonitorHelper.TH1F
def TH1F(name, title, nxbins, bins_par2, bins_par3=None, path='', **kwargs)
Definition: TrigEgammaMonitorHelper.py:24
DiMuMon::m_varsDistr
std::vector< std::string > m_varsDistr
Definition: DiMuMon.h:76
xAOD::TrackParticle_v1
Class describing a TrackParticle.
Definition: TrackParticle_v1.h:43
plotBeamSpotCompare.histo
histo
Definition: plotBeamSpotCompare.py:415
checker_macros.h
Define macros for attributes used to control the static checker.
GeV
#define GeV
Definition: CaloTransverseBalanceVecMon.cxx:30
ConstDataVector::begin
iterator begin() noexcept
Return an iterator pointing at the beginning of the collection.
DiMuMon::m_nMassBins
int m_nMassBins
Definition: DiMuMon.h:51
DiMuMon::getInvmass
double getInvmass(const xAOD::TrackParticle *track1, const xAOD::TrackParticle *track2, double Mass) const
Definition: DiMuMon.cxx:627
TSU::T
unsigned long long T
Definition: L1TopoDataTypes.h:35
xAOD::TrackParticle_v1::phi
virtual double phi() const override final
The azimuthal angle ( ) of the particle (has range to .)
xAOD::numberOfInnermostPixelLayerHits
@ numberOfInnermostPixelLayerHits
these are the hits in the 0th pixel barrel layer
Definition: TrackingPrimitives.h:237
ManagedMonitorToolBase::newRunFlag
bool newRunFlag() const
Definition: ManagedMonitorToolBase.h:854
DiMuMon::m_varLabels
std::map< std::string, std::string > m_varLabels
Definition: DiMuMon.h:79