ATLAS Offline Software
MonitoringFile_DiMuPostProcess.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 "DataQualityUtils/MonitoringFile.h"
6 
7 #include <cmath>
8 #include <vector>
9 #include <sstream>
10 
11 #include <TCanvas.h>
12 #include <TF1.h>
13 #include <TFile.h>
14 #include <TH1.h>
15 #include <TH2.h>
16 #include <TKey.h>
17 #include <TMath.h>
18 #include <TProfile.h>
19 #include <TTree.h>
20 #include <TROOT.h>
21 
22 #include "RooRealVar.h"
23 #include "RooDataHist.h"
24 #include "RooArgList.h"
25 #include "RooPlot.h"
26 #include "RooBreitWigner.h"
27 #include "RooCBShape.h"
28 #include "RooFFTConvPdf.h"
29 #include "RooGlobalFunc.h"
30 #include "RooFitResult.h"
31 
32 namespace dqutils {
33  void
34  MonitoringFile::
35  fitMergedFile_DiMuMonManager(const std::string& inFilename, bool isIncremental) { //adapted from
36  // MonitoringFile_IDPerfPostProcess.cxx
37  if (isIncremental == true) return;
38 
39  TFile* f = TFile::Open(inFilename.c_str(), "UPDATE");
40  if (f == 0 || !f->IsOpen()) {
41  //std::cout<<"Failed to load or open file!"<<endl;
42  delete f;
43  return;
44  }
45  if (f->GetSize() < 1000.) {
46  //std::cout<<"File is empty!"<<endl;
47  delete f;
48  return;
49  }
50  std::string run_dir;
51  TIter next_run(f->GetListOfKeys());
52  TKey* key_run(0);
53  while ((key_run = dynamic_cast<TKey*>(next_run())) != 0) {
54  TObject* obj_run = key_run->ReadObj();
55  TDirectory* tdir_run = dynamic_cast<TDirectory*>(obj_run);
56  if (tdir_run != 0) {
57  std::string tdir_run_name(tdir_run->GetName());
58  if (tdir_run_name.find("run") != std::string::npos) {
59  run_dir = tdir_run_name;
60  TIter next_perf(tdir_run->GetListOfKeys());
61  TKey* key_perf(0);
62  while ((key_perf = dynamic_cast<TKey*>(next_perf())) != 0) {
63  std::string perf_name(key_perf->GetName());
64  if (perf_name.find("DiMuMon") != std::string::npos) {
65  TObject* obj_perf = key_perf->ReadObj();
66  TDirectory* tdir_perf = dynamic_cast<TDirectory*>(obj_perf);
67  if (tdir_perf != 0) {
68  run_dir = run_dir + '/';
69  TIter next_module(tdir_perf->GetListOfKeys());
70  TKey* key_module(0);
71  while ((key_module = dynamic_cast<TKey*>(next_module())) != 0) {
72  std::string module_name(key_module->GetName());
73  if (module_name.find("Jpsi") != std::string::npos) {
74  processModule(f, run_dir, key_module, "Jpsi");
75  } else if (module_name.find("Zmumu") != std::string::npos) {
76  processModule(f, run_dir, key_module, "Zmumu");
77  }
78  }
79  } else {
80  delete obj_perf;
81  }
82  }
83  }
84  }
85  // if without top run_directory
86  else if (tdir_run_name.find("DiMuMon") != std::string::npos) {
87  TObject* obj_perf = key_run->ReadObj();
88  TDirectory* tdir_perf = dynamic_cast<TDirectory*>(obj_perf);
89  if (tdir_perf != 0) {
90  run_dir = '/';
91  TIter next_module(tdir_perf->GetListOfKeys());
92  TKey* key_module(0);
93  while ((key_module = dynamic_cast<TKey*>(next_module())) != 0) {
94  std::string module_name(key_module->GetName());
95  if (module_name.find("Jpsi") != std::string::npos) {
96  processModule(f, run_dir, key_module, "Jpsi");
97  } else if (module_name.find("Zmumu") != std::string::npos) {
98  processModule(f, run_dir, key_module, "Zmumu");
99  }
100  }
101  } else {
102  delete obj_perf;
103  }
104  }
105  } else {
106  delete obj_run;
107  }
108  }
109 
110  f->Close();
111  delete f;
112  // if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "DiMuPostProcessing done!" <<endmsg;
113  }
114 
115  void MonitoringFile::processModule(TFile* f, const std::string& run_dir, TKey* key_module,
116  const std::string& moduleName) {
117  TObject* obj_mod = key_module->ReadObj();
118  TDirectory* tdir_mod = dynamic_cast<TDirectory*>(obj_mod);
119 
120  if (tdir_mod != 0) {
121  TIter next_trigger(tdir_mod->GetListOfKeys());
122  TKey* key_trigger(0);
123  while ((key_trigger = dynamic_cast<TKey*>(next_trigger())) != 0) {
124  std::string triggerName = key_trigger->GetName();
125  fitMergedFile_DiMuMonAll(f, run_dir, moduleName, triggerName);
126  }
127  } else {
128  delete obj_mod;
129  }
130  }
131 
132  void
133  MonitoringFile::fitMergedFile_DiMuMonAll(TFile* f, const std::string& run_dir, const std::string& resonName,
134  const std::string& triggerName) {
135  //std::cout<<"fitMergedFile_DiMuMon has been called"<<endl;
136  std::string path;
137  path = run_dir + "DiMuMon/" + resonName + "/" + triggerName;
138 
139  if (f->cd(path.c_str()) == 0) {
140  return;
141  }
142  std::vector<std::string> regions;
143  regions.push_back("All");
144  regions.push_back("BB");
145  regions.push_back("EAEA");
146  regions.push_back("ECEC");
147  std::vector<std::string> vars;
148  vars.push_back("eta");
149  vars.push_back("etaAll");
150  vars.push_back("etaPos");
151  vars.push_back("etaNeg");
152  vars.push_back("phi");
153  vars.push_back("phiAll");
154  vars.push_back("phiPos");
155  vars.push_back("phiNeg");
156  vars.push_back("pt");
157  vars.push_back("ptAll");
158  vars.push_back("ptPos");
159  vars.push_back("ptNeg");
160  vars.push_back("etaDiff");
161  vars.push_back("etaSumm");
162  vars.push_back("phiDiff");
163  vars.push_back("phiSumm");
164  vars.push_back("crtDiff");
165 
166  std::map< std::string, TH1F* > h_invmass;
167  std::map< std::string, std::map< std::string, TH2F*> > h_2DinvmassVSx;
168  std::map< std::string, std::map< std::string, TH1F*> > h_invmassVSx;
169  std::map< std::string, std::map< std::string, TH1F*> > h_widthVSx;
170  TH1F* h_chi2;
171 
172  //loop over all possible 2D histos
173  //check if 2D histo has been filled
174  //if found the 2D histo, then see whether the mean or width or both 1D histos were also made.-->Decide what to refit
175  // `
176  if (CheckHistogram(f, (path + "_detail/chi2").c_str())) {
177  h_chi2 = (TH1F*) (f->Get((path + "_detail/chi2").c_str())->Clone());
178  std::vector<std::string> ::iterator ivar = vars.begin();
179  std::vector<std::string> ::iterator ireg = regions.begin();
180  for (ireg = regions.begin(); ireg != regions.end(); ++ireg) {
181  for (ivar = vars.begin(); ivar != vars.end(); ++ivar) {
182  std::string hname2D = resonName + "_2DinvmassVS" + *ivar + "_" + *ireg;
183  if (CheckHistogram(f, (path + "/" + hname2D).c_str())) {
184  h_2DinvmassVSx[*ireg][*ivar] = (TH2F*) (f->Get((path + "/" + hname2D).c_str())->Clone());
185  std::string hnameMean = resonName + "_invmassVS" + *ivar + "_" + *ireg;
186  std::string hnameWidth = resonName + "_widthVS" + *ivar + "_" + *ireg;
187  std::vector<TH1F*> hfitted;
188  if (CheckHistogram(f, (path + "/" + hnameMean).c_str())) {
189  h_invmassVSx[*ireg][*ivar] = (TH1F*) (f->Get((path + "/" + hnameMean).c_str())->Clone());
190  hfitted.push_back(h_invmassVSx[*ireg][*ivar]);
191  if (CheckHistogram(f, (path + "_detail/" + hnameWidth).c_str())) {
192  h_widthVSx[*ireg][*ivar] = (TH1F*) (f->Get((path + "_detail/" + hnameWidth).c_str())->Clone());
193  hfitted.push_back(h_widthVSx[*ireg][*ivar]);
194  fitHistos(h_2DinvmassVSx[*ireg][*ivar], hfitted, 0, triggerName, resonName, h_chi2);// 0 means to fill
195  // both mean and
196  // width results
197  // from the fit
198  f->cd((path + "/").c_str());
199  h_invmassVSx[*ireg][*ivar]->Write("", TObject::kOverwrite);
200  f->cd((path + "_detail/").c_str());
201  h_widthVSx[*ireg][*ivar]->Write("", TObject::kOverwrite);
202  } else {
203  fitHistos(h_2DinvmassVSx[*ireg][*ivar], hfitted, 1, triggerName, resonName, h_chi2);// 1 means to fill
204  // only mean results
205  // from the fit
206  f->cd((path + "/").c_str());
207  h_invmassVSx[*ireg][*ivar]->Write("", TObject::kOverwrite);
208  }
209  } else {
210  if (CheckHistogram(f, (path + "_detail/" + hnameWidth).c_str())) {
211  h_widthVSx[*ireg][*ivar] = (TH1F*) (f->Get((path + "_detail/" + hnameWidth).c_str())->Clone());
212  hfitted.push_back(h_widthVSx[*ireg][*ivar]);
213  fitHistos(h_2DinvmassVSx[*ireg][*ivar], hfitted, 2, triggerName, resonName, h_chi2);// 2 means to fill
214  // only width
215  // results from the
216  // fit
217  f->cd((path + "_detail/").c_str());
218  h_widthVSx[*ireg][*ivar]->Write((path + "_detail/" + hnameWidth).c_str(), TObject::kOverwrite);
219  }
220  }
221  }
222  }
223  }
224  h_chi2->Write("", TObject::kOverwrite);
225  }
226  f->Write();
227  }
228 
229  void MonitoringFile::fitHistos(TH2F* hin, const std::vector<TH1F*>& hout, int mode, const std::string& triggerName,
230  const std::string& resonName, TH1F* h_chi2) {
231  bool saveHistos = false;
232 
233  // a canvas may be needed when implmenting this into the post-processing file
234  //std::cout<<"The fitHistos method is called"<<endl;
235  std::string hname = hin->GetName();
236  TCanvas* ctemp;
237  char num2str[50];
238  int nbins = hin->GetNbinsX();
239  for (int i = 0; i < nbins; i++) {
240  snprintf(num2str, 50, "%s_%i", (hname).c_str(), i);
241  TH1F* htemp = (TH1F*) (hin->ProjectionY(num2str, i + 1, i + 1));
242  //htemp->SetTitle(projName);
243  htemp->Sumw2();
244  if (htemp->GetEntries() > 50) {
245  double mean = 999., meanErr = 999., sigma = 999., sigmaErr = 999., chi2 = 0.;
246  if (resonName == "Jpsi" || resonName == "Upsi") {
247  mean = htemp->GetMean();
248  sigma = htemp->GetRMS();
249  TF1* fn = new TF1("fn", "gaus", mean - 2 * sigma, mean + 2 * sigma);
250  fn->SetParameters(htemp->GetEntries() / 4, mean, sigma);
251  htemp->Fit("fn", "RQMN");
252  mean = fn->GetParameter(1);
253  sigma = fn->GetParameter(2);
254  fn->SetRange(mean - 1.2 * sigma, mean + 1.2 * sigma);
255 
256  if (saveHistos == false) {
257  htemp->Fit("fn", "RQMN");
258  } else {
259  ctemp = new TCanvas("ctemp", "ctemp", 500, 500);
260  TString psName = num2str + triggerName + ".ps";
261  htemp->Fit("fn", "RQM");
262  ctemp->Print(psName);
263  delete ctemp;
264  }
265  double frange = 2.4 * sigma;
266  double hrange = htemp->GetXaxis()->GetXmax() - htemp->GetXaxis()->GetXmin();
267  double ndf = frange / hrange * (htemp->GetNbinsX()) - 3;//subtract number of fit parameters
268  //fill results ;
269  mean = fn->GetParameter(1);
270  meanErr = fn->GetParError(1);
271  sigma = fn->GetParameter(2);
272  sigmaErr = fn->GetParError(2);
273  chi2 = (fn->GetChisquare()) / ndf;
274  delete fn;
275  } else {
276  //fit Z peak with a convolution of BreitWigner and Crystal Ball fns, fit by Louise, implementation by Jike
277  // taken from IDPerfMon
278  RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
279  RooAbsPdf::verboseEval(-100); //sami in order to make roofit quiet
280  RooRealVar m("mass", "dimuon invariable mass", 91.2, 71., 111., "GeV");
281  RooDataHist* data = 0;
282  data = new RooDataHist("data", "data", m, htemp);
283  RooRealVar bwm0("bw_#mu", "bw_#mu", 91.2, 85.2, 97.2);
284  RooRealVar bwsg("bw_#sigma", "bw_#sigma", 2.4952);
285  RooBreitWigner bw("bw", "bw", m, bwm0, bwsg);
286 
287  RooRealVar cbm0("cb_#mu", "cb_#mu", 0);
288  RooRealVar cbsg("cb_#sigma", "cb_#sigma", 3., 1., 10.);
289  RooRealVar cbal("cb_#alpha", "cb_#alpha", 2.0);
290  RooRealVar cbn("cb_n", "cb_n", 1., 0.05, 3.);
291  RooCBShape cb("cb", "cb", m, cbm0, cbsg, cbal, cbn);
292 
293  m.setBins(5000);
294  RooFFTConvPdf bxc("bxc", "BW (X) CB", m, bw, cb);
295  auto resourceToDelete = bxc.fitTo(*data, RooFit::PrintLevel(-1), RooFit::PrintEvalErrors(-1), RooFit::Warnings(kFALSE));
296  RooPlot* frame;
297  if (saveHistos == true) {
298  frame = m.frame();
299  data->plotOn(frame, RooFit::MarkerSize(0.9));
300  bxc.paramOn(frame, RooFit::Format("NELU", RooFit::AutoPrecision(2)), RooFit::Layout(0.1, 0.4, 0.9));
301  bxc.plotOn(frame, RooFit::LineColor(kBlue));
302  delete frame;
303  }
304  delete resourceToDelete;
305  mean = bwm0.getVal();
306  meanErr = bwm0.getError();
307  sigma = cbsg.getVal();
308  sigmaErr = cbsg.getError();
309  delete data;
310  }
311  //fill results
312  if (mode == 0) {//plot both invmass and width vs variable
313  hout.at(0)->SetBinContent(i + 1, mean);
314  hout.at(0)->SetBinError(i + 1, meanErr);
315  hout.at(1)->SetBinContent(i + 1, sigma);
316  hout.at(1)->SetBinError(i + 1, sigmaErr);
317  } else if (mode == 1) {//plot only invmass vs variable
318  hout.at(0)->SetBinContent(i + 1, mean);
319  hout.at(0)->SetBinError(i + 1, meanErr);
320  } else if (mode == 2) {//plot only width vs variable
321  hout.at(0)->SetBinContent(i + 1, sigma);
322  hout.at(0)->SetBinError(i + 1, sigmaErr);
323  }
324  h_chi2->Fill(chi2);
325  }// more than 50 events
326 
327  delete htemp;
328  }
329  }
330 }
xAOD::iterator
JetConstituentVector::iterator iterator
Definition: JetConstituentVector.cxx:68
covarianceTool.ndf
ndf
Definition: covarianceTool.py:678
dqutils::MonitoringFile::processModule
static void processModule(TFile *f, const std::string &run_dir, TKey *key_module, const std::string &moduleName)
Definition: MonitoringFile_DiMuPostProcess.cxx:115
data
char data[hepevt_bytes_allocation_ATLAS]
Definition: HepEvt.cxx:11
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
python.SystemOfUnits.m
int m
Definition: SystemOfUnits.py:91
athena.path
path
python interpreter configuration --------------------------------------—
Definition: athena.py:128
dqutils::MonitoringFile::fitHistos
static void fitHistos(TH2F *hin, const std::vector< TH1F * > &hout, int mode, const std::string &triggerName, const std::string &resonName, TH1F *m_chi2)
Definition: MonitoringFile_DiMuPostProcess.cxx:229
python.Constants.FATAL
int FATAL
Definition: Control/AthenaCommon/python/Constants.py:19
dqutils::MonitoringFile::CheckHistogram
static bool CheckHistogram(TFile *f, const char *HistoName)
dqt_zlumi_pandas.hname
string hname
Definition: dqt_zlumi_pandas.py:279
StateLessPT_NewConfig.Format
Format
Definition: StateLessPT_NewConfig.py:146
dqutils::MonitoringFile::fitMergedFile_DiMuMonManager
static void fitMergedFile_DiMuMonManager(const std::string &inFileName, bool isIncremental=false)
Definition: MonitoringFile_DiMuPostProcess.cxx:35
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
instance
std::map< std::string, double > instance
Definition: Run_To_Get_Tags.h:8
dqutils::MonitoringFile::fitMergedFile_DiMuMonAll
static void fitMergedFile_DiMuMonAll(TFile *f, const std::string &run_dir, const std::string &resonName, const std::string &triggerName)
Definition: MonitoringFile_DiMuPostProcess.cxx:133
python.getCurrentFolderTag.fn
fn
Definition: getCurrentFolderTag.py:65
lumiFormat.i
int i
Definition: lumiFormat.py:85
plotBeamSpotCompare.ivar
int ivar
Definition: plotBeamSpotCompare.py:383
chi2
double chi2(TH1 *h0, TH1 *h1)
Definition: comparitor.cxx:523
Preparation.mode
mode
Definition: Preparation.py:94
hist_file_dump.f
f
Definition: hist_file_dump.py:135
dqutils
Definition: CoolMdt.h:76
SCT_CalibAlgs::nbins
@ nbins
Definition: SCT_CalibNumbers.h:10
MonitoringFile.h
StateLessPT_NewConfig.Layout
Layout
Definition: StateLessPT_NewConfig.py:159
python.TrigEgammaMonitorHelper.TH1F
def TH1F(name, title, nxbins, bins_par2, bins_par3=None, path='', **kwargs)
Definition: TrigEgammaMonitorHelper.py:24
DQHistogramMerge.isIncremental
isIncremental
Definition: DQHistogramMerge.py:37
Generated on Sun Dec 22 2024 21:14:41 for ATLAS Offline Software by   doxygen 1.8.18