AliPhysics  1168478 (1168478)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MuonTrackingEfficiency.C
Go to the documentation of this file.
1 
29 #include <TROOT.h>
30 #include <TH1.h>
31 #include <TH2.h>
32 #include <THnSparse.h>
33 #include <TAxis.h>
34 #include <TString.h>
35 #include <TObjString.h>
36 #include <Riostream.h>
37 #include <TFile.h>
38 #include <TList.h>
39 #include <TCanvas.h>
40 #include <TGraphAsymmErrors.h>
41 #include <TMath.h>
42 #include <TArrayD.h>
43 #include <TStyle.h>
44 #include <THashList.h>
45 #include <TParameter.h>
46 
47 //const Char_t *effErrMode = "cp"; // Clopper-Pearson
48 const Char_t *effErrMode = "b(1,1)mode"; // Bayesian with uniform prior
49 
50 Double_t centMin = -999.;
54 Int_t charge = 0; // 0 selects + and -, -1 and +1 selects - or + muons respectively
55 
57 
58 THashList *runWeights = 0x0;
60 
61 
62 void PlotMuonEfficiencyVsX(TString var, TString fileNameData, TString fileNameSave, Bool_t saveEdges, Bool_t print, Bool_t draw);
63 void PlotIntegratedMuonEfficiencyVsX(TString var, TString runList, TString fileNameWeights,
64  TString fileNameData, TString fileNameSave, Bool_t print, Bool_t draw);
65 
66 void PlotMuonEfficiencyVsXY(TString xVar, TString yVar, TString fileNameData, TString fileNameSave, Bool_t draw, Bool_t rap = kFALSE);
67 
68 void PlotMuonEfficiency(TString fileNameData, TString fileNameSave, Bool_t saveEdges, Bool_t print, Bool_t draw);
69 void PlotMuonEfficiencyVsRun(TString runList, TString fileNameData, TString fileNameSave, Bool_t print, Bool_t draw);
70 void PlotIntegratedMuonEfficiency(TString fileNameWeights, TString fileNameSave, Bool_t print, Bool_t draw);
71 
72 void PlotMuonEfficiencyPerDE(TString fileNameData, TString fileNameSave, Bool_t saveEdges);
73 void PlotMuonEfficiencyPerDEVsRun(TString runList, TString fileNameData, TString fileNameSave);
74 void PlotIntegratedMuonEfficiencyPerDE(TString fileNameWeights, TString fileNameSave);
75 
76 Bool_t GetChamberEfficiency(THnSparse &TT, THnSparse &TD, TArrayD &chEff, TArrayD chEffErr[2], Bool_t printError = kFALSE);
77 void GetDEEfficiency(THnSparse &TT, THnSparse &TD, TGraphAsymmErrors &effVsDE);
78 void ComputeStationEfficiency(TArrayD &chEff, TArrayD chEffErr[2], Int_t iSt, Double_t &stEff, Double_t stEffErr[2]);
79 void GetStationEfficiency(TArrayD &chEff, TArrayD chEffErr[2], Int_t iSt, TGraphAsymmErrors *effVsX[3], Int_t ip, Double_t xp);
80 void ComputeStation45Efficiency(TArrayD &chEff, TArrayD chEffErr[2], Double_t &st45Eff, Double_t st45EffErr[2]);
81 void GetStation45Efficiency(TArrayD &chEff, TArrayD chEffErr[2], TGraphAsymmErrors *effVsX[3], Int_t ip, Double_t xp);
82 void ComputeTrackingEfficiency(Double_t stEff[6], Double_t stEffErr[6][2], Double_t &spectroEff, Double_t spectroEffErr[2]);
83 void GetTrackingEfficiency(TArrayD &chEff, TArrayD chEffErr[2], TGraphAsymmErrors *effVsSt[3],
84  TGraphAsymmErrors *effVsX[3], Int_t ip, Double_t xp, Bool_t print = kFALSE);
85 void IntegrateMuonEfficiency(TGraphAsymmErrors &effVsRunLow, TGraphAsymmErrors &effVsRunUp,
86  TGraphAsymmErrors &effVsX, Int_t ip, Double_t xp);
87 
89 void SetCentPtCh(THnSparse& SparseData);
90 TGraphAsymmErrors* CreateGraph(const char* name, const char* title, int value=-1);
91 void BeautifyGraph(TGraphAsymmErrors &g, const char* xAxisName, const char* yAxisName);
92 void BeautifyGraphs(TObjArray& array, const char* xAxisName, const char* yAxisName);
93 void SetRunLabel(TGraphAsymmErrors &g, Int_t irun, const TList& runs);
94 void SetRunLabel(TObjArray& array, Int_t irun, const TList& runs);
95 
96 
97 //---------------------------------------------------------------------------
98 void MuonTrackingEfficiency(TString runList = "runList.txt",
99  TString fileNameWeights = "",
100  TString objNameExtension = "",
101  TString fileNameData ="AnalysisResults.root",
102  TString fileNameSave = "efficiency_new.root")
103 {
105 
106  gObjNameExtension = objNameExtension;
107 
108  PlotMuonEfficiencyVsX("centrality", fileNameData, fileNameSave, kFALSE, kFALSE, kTRUE);
109  PlotMuonEfficiencyVsX("pt", fileNameData, fileNameSave, kFALSE, kFALSE, kTRUE);
110  PlotMuonEfficiencyVsX("y", fileNameData, fileNameSave, kFALSE, kFALSE, kTRUE);
111  PlotMuonEfficiencyVsX("phi", fileNameData, fileNameSave, kFALSE, kFALSE, kTRUE);
112  PlotMuonEfficiencyVsX("charge", fileNameData, fileNameSave, kFALSE, kFALSE, kTRUE);
113 
114  PlotIntegratedMuonEfficiencyVsX("centrality", runList, fileNameWeights, fileNameData, fileNameSave, kFALSE, kTRUE);
115  PlotIntegratedMuonEfficiencyVsX("pt", runList, fileNameWeights, fileNameData, fileNameSave, kFALSE, kTRUE);
116  PlotIntegratedMuonEfficiencyVsX("y", runList, fileNameWeights, fileNameData, fileNameSave, kFALSE, kTRUE);
117  PlotIntegratedMuonEfficiencyVsX("phi", runList, fileNameWeights, fileNameData, fileNameSave, kFALSE, kTRUE);
118  PlotIntegratedMuonEfficiencyVsX("charge", runList, fileNameWeights, fileNameData, fileNameSave, kFALSE, kTRUE);
119 
120  PlotMuonEfficiencyVsXY("pt", "centrality", fileNameData, fileNameSave, kTRUE);
121  PlotMuonEfficiencyVsXY("y", "centrality", fileNameData, fileNameSave, kTRUE);
122  PlotMuonEfficiencyVsXY("pt", "y", fileNameData, fileNameSave, kTRUE);
123  PlotMuonEfficiencyVsXY("phi", "y", fileNameData, fileNameSave, kTRUE, kTRUE);
124 
125  PlotMuonEfficiency(fileNameData, fileNameSave, kFALSE, kTRUE, kTRUE);
126  PlotMuonEfficiencyVsRun(runList, fileNameData, fileNameSave, kFALSE, kTRUE);
127  PlotIntegratedMuonEfficiency(fileNameWeights, fileNameSave, kTRUE, kTRUE);
128 
129  PlotMuonEfficiencyPerDE(fileNameData, fileNameSave, kFALSE);
130  PlotMuonEfficiencyPerDEVsRun(runList, fileNameData, fileNameSave);
131  PlotIntegratedMuonEfficiencyPerDE(fileNameWeights, fileNameSave);
132 
133 }
134 
135 
136 //---------------------------------------------------------------------------
137 void PlotMuonEfficiencyVsX(TString var, TString fileNameData, TString fileNameSave, Bool_t saveEdges, Bool_t print, Bool_t draw)
138 {
140 
141  printf("plotting efficiency versus %s...\n", var.Data());
142 
143  Int_t xDim = -1;
144  if (var == "centrality") xDim = 1;
145  else if (var == "pt") xDim = 2;
146  else if (var == "y") xDim = 3;
147  else if (var == "phi") xDim = 4;
148  else if (var == "charge") xDim = 5;
149  else {
150  printf("incorrect variable. Choices are centrality, pt, y, phi and charge.\n");
151  return;
152  }
153 
154  // get input hists
155  TFile *file = new TFile(fileNameData.Data(), "read");
156  if (!file || !file->IsOpen()) {
157  printf("cannot open file %s \n",fileNameData.Data());
158  return;
159  }
160  TList *listTT = static_cast<TList*>(file->FindObjectAny(Form("TotalTracksPerChamber%s", gObjNameExtension.Data())));
161  TList *listTD = static_cast<TList*>(file->FindObjectAny(Form("TracksDetectedPerChamber%s", gObjNameExtension.Data())));
162  THnSparse *TT = static_cast<THnSparse*>(listTT->At(10));
163  THnSparse *TD = static_cast<THnSparse*>(listTD->At(10));
164 
165  // output graph
166  TGraphAsymmErrors *effVsX[3] = {0x0, 0x0, 0x0};
167  TString nameAdd[3] = {"", "Low", "Up"};
168  TString titleAdd[3] = {"", " - lower limit", " - upper limit"};
169  for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i) {
170  effVsX[i] = CreateGraph(Form("trackingEffVs%s%s",var.Data(),nameAdd[i].Data()),
171  Form("Measured tracking efficiency versus %s%s",var.Data(),titleAdd[i].Data()));
172  }
173 
174  // set the centrality and pT range for integration
175  SetCentPtCh(*TT);
176  SetCentPtCh(*TD);
177 
178  // loop over X bins
179  TArrayD chEff(11);
180  TArrayD chEffErr[2];
181  chEffErr[0].Set(11);
182  chEffErr[1].Set(11);
183  for (Int_t ix = 1; ix <= TT->GetAxis(xDim)->GetNbins(); ix++) {
184 
185  if (print) cout << var.Data() << " " << TT->GetAxis(xDim)->GetBinLowEdge(ix) << "-" << TT->GetAxis(xDim)->GetBinUpEdge(ix) << ":" << endl;
186 
187  // set the var range to the current bin
188  TT->GetAxis(xDim)->SetRange(ix, ix);
189  TD->GetAxis(xDim)->SetRange(ix, ix);
190 
191  // compute chamber efficency and errors
192  if (GetChamberEfficiency(*TT, *TD, chEff, chEffErr, print)) {
193 
194  // compute the overall tracking efficiency
195  TGraphAsymmErrors *dummy[3] = {0x0, 0x0, 0x0};
196  GetTrackingEfficiency(chEff, chEffErr, dummy, effVsX, ix-1, TT->GetAxis(xDim)->GetBinCenter(ix), print);
197 
198  } else {
199 
200  // fill graph with 1 +0 -1
201  effVsX[0]->SetPoint(ix-1,TT->GetAxis(xDim)->GetBinCenter(ix),1.);
202  effVsX[0]->SetPointError(ix-1,0.,0.,1.,0.);
203 
204  if (saveEdges) {
205 
206  // lower = 0 ± 0
207  effVsX[1]->SetPoint(ix-1,TT->GetAxis(xDim)->GetBinCenter(ix),0.);
208  effVsX[1]->SetPointError(ix-1,0.,0.,0.,0.);
209 
210  // upper = 1 ± 0
211  effVsX[2]->SetPoint(ix-1,TT->GetAxis(xDim)->GetBinCenter(ix),1.);
212  effVsX[2]->SetPointError(ix-1,0.,0.,0.,0.);
213 
214  }
215 
216  }
217 
218  }
219 
220  // close input file
221  file->Close();
222 
223  // display
224  for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i) BeautifyGraph(*effVsX[i], var.Data(), "efficiency");
225 
226  if (draw) {
227  new TCanvas(Form("cTrackingEffVs%s",var.Data()), Form("Measured tracking efficiency versus %s",var.Data()),1000,400);
228  effVsX[0]->DrawClone("ap");
229  }
230 
231  // save output
232  file = new TFile(fileNameSave.Data(),"update");
233  for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i) effVsX[i]->Write(0x0, TObject::kOverwrite);
234  file->Close();
235 
236  // clean memory
237  for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i) delete effVsX[i];
238 
239 }
240 
241 
242 //---------------------------------------------------------------------------
243 void PlotIntegratedMuonEfficiencyVsX(TString var, TString runList, TString fileNameWeights,
244  TString fileNameData, TString fileNameSave, Bool_t print, Bool_t draw)
245 {
247 
248  printf("plotting integrated efficiency versus %s...\n", var.Data());
249 
250  // load run weights
251  if (!fileNameWeights.IsNull()) LoadRunWeights(fileNameWeights);
252  if (!runWeights) {
253  printf("Cannot compute integrated efficiency without run-by-run weights\n");
254  return;
255  }
256 
257  // open run list
258  ifstream inFile(runList.Data());
259  if (!inFile.is_open()) {
260  printf("cannot open file %s\n",runList.Data());
261  return;
262  }
263 
264  // output graph
265  TGraphAsymmErrors *intEffVsX = CreateGraph(Form("integratedTrackingEffVs%s",var.Data()),
266  Form("Integrated tracking efficiency versus %s",var.Data()));
267 
268  Int_t n = -1;
269  TArrayD x;
270  TArrayD rec[2];
271  Double_t gen = 0.;
272  TArrayD effErr[2];
273 
274  // loop over runs
275  while (!inFile.eof()) {
276 
277  // get current run number
278  TString currRun;
279  currRun.ReadLine(inFile,kTRUE);
280  if(currRun.IsNull() || !currRun.IsDec()) continue;
281  Int_t run = currRun.Atoi();
282 
283  printf("run %d: ", run);
284 
285  // compute efficiency vs var
286  TString dataFile = Form("runs/%d/%s", run, fileNameData.Data());
287  TString outFile = Form("runs/%d/%s", run, fileNameSave.Data());
288  PlotMuonEfficiencyVsX(var, dataFile, outFile, kTRUE, print, kFALSE);
289 
290  // get run weight
291  TParameter<Double_t> *weight = static_cast<TParameter<Double_t>*>(runWeights->FindObject(currRun.Data()));
292  if (!weight) {
293  printf("weight not found for run %s\n", currRun.Data());
294  continue;
295  }
296  Double_t w = weight->GetVal();
297  Double_t w2 = w*w;
298 
299  // get result
300  TFile *file = new TFile(outFile.Data(), "read");
301  if (!file || !file->IsOpen()) {
302  printf("cannot open file\n");
303  continue;
304  }
305  TGraphAsymmErrors *effVsX[2];
306  effVsX[0] = static_cast<TGraphAsymmErrors*>(file->FindObjectAny(Form("trackingEffVs%sLow",var.Data())));
307  effVsX[1] = static_cast<TGraphAsymmErrors*>(file->FindObjectAny(Form("trackingEffVs%sUp",var.Data())));
308  if (!effVsX[0] || !effVsX[1]) {
309  printf("trackingEffVs%sLow(Up) object not found\n", var.Data());
310  continue;
311  }
312 
313  // prepare the arrays if not already done
314  if (n < 0) {
315  n = effVsX[0]->GetN();
316  x.Set(n, effVsX[0]->GetX());
317  for (Int_t i = 0; i < 2; ++i) {
318  rec[i].Set(n);
319  effErr[i].Set(n);
320  }
321  } else if (n != effVsX[0]->GetN()) {
322  printf("number of points in graph trackingEffVs%sLow(Up) for run %d is different than from previous runs\n", var.Data(), run);
323  continue;
324  }
325 
326  // integrate for all bins
327  gen += w;
328  for (Int_t ix = 0; ix < n; ++ix) {
329  Double_t ieffErr[2] = {effVsX[0]->GetErrorYlow(ix), effVsX[1]->GetErrorYhigh(ix)};
330  for (Int_t i = 0; i < 2; ++i) {
331  rec[i][ix] += w*effVsX[i]->GetY()[ix];
332  effErr[i][ix] += w2*ieffErr[i]*ieffErr[i];
333  }
334  }
335 
336  // close input file
337  file->Close();
338 
339  }
340 
341  inFile.close();
342 
343  // fill output graph
344  if (gen > 0.) {
345 
346  for (Int_t ix = 0; ix < n; ++ix) {
347 
348  intEffVsX->SetPoint(ix, x[ix], rec[1][ix]/gen);
349  intEffVsX->SetPointError(ix, 0., 0., (rec[1][ix]-rec[0][ix]+TMath::Sqrt(effErr[0][ix]))/gen, TMath::Sqrt(effErr[1][ix])/gen);
350 
351  }
352 
353  } else {
354 
355  for (Int_t ix = 0; ix < n; ++ix) {
356 
357  printf("impossible to integrate, all weights = 0 or unknown ?!?\n");
358 
359  intEffVsX->SetPoint(ix, x[ix], -1.);
360  intEffVsX->SetPointError(ix, 0., 0., 0., 0.);
361 
362  }
363 
364  }
365 
366  // display
367  BeautifyGraph(*intEffVsX, var.Data(), "efficiency");
368 
369  if (draw) {
370  new TCanvas(Form("cIntegratedTrackingEffVs%s",var.Data()), Form("Integrated tracking efficiency versus %s",var.Data()),1000,400);
371  intEffVsX->DrawClone("ap");
372  }
373 
374  // save output
375  TFile *file = new TFile(fileNameSave.Data(),"update");
376  intEffVsX->Write(0x0, TObject::kOverwrite);
377  file->Close();
378 
379  // clean memory
380  delete intEffVsX;
381 
382 }
383 
384 
385 //---------------------------------------------------------------------------
386 void PlotMuonEfficiencyVsXY(TString xVar, TString yVar, TString fileNameData, TString fileNameSave, Bool_t draw, Bool_t rap)
387 {
389 
390  printf("plotting efficiency versus %s/%s...\n", xVar.Data(), yVar.Data());
391 
392  Int_t xDim = -1;
393  if (xVar == "centrality") xDim = 1;
394  else if (xVar == "pt") xDim = 2;
395  else if (xVar == "y") xDim = 3;
396  else if (xVar == "phi") xDim = 4;
397  else if (xVar == "charge") xDim = 5;
398  else {
399  printf("incorrect variable. Choices are centrality, pt, y, phi and charge.\n");
400  return;
401  }
402  Int_t yDim = -1;
403  if (yVar == "centrality") yDim = 1;
404  else if (yVar == "pt") yDim = 2;
405  else if (yVar == "y") yDim = 3;
406  else if (yVar == "phi") yDim = 4;
407  else if (yVar == "charge") yDim = 5;
408  else {
409  printf("incorrect variable. Choices are centrality, pt, y, phi and charge.\n");
410  return;
411  }
412 
413  // get input hists
414  TFile *file = new TFile(fileNameData.Data(), "read");
415  if (!file || !file->IsOpen()) {
416  printf("cannot open file %s\n",fileNameData.Data());
417  return;
418  }
419  TList *listTT = static_cast<TList*>(file->FindObjectAny(Form("TotalTracksPerChamber%s", gObjNameExtension.Data())));
420  TList *listTD = static_cast<TList*>(file->FindObjectAny(Form("TracksDetectedPerChamber%s", gObjNameExtension.Data())));
421  THnSparse *TT = static_cast<THnSparse*>(listTT->At(10));
422  THnSparse *TD = static_cast<THnSparse*>(listTD->At(10));
423 
424  // output map
425  Int_t nxBins = TT->GetAxis(xDim)->GetNbins();
426  Int_t nyBins = TT->GetAxis(yDim)->GetNbins();
427  TH2F *effVsXY = new TH2F(Form("trackingEffVs%s-%s",xVar.Data(),yVar.Data()),
428  Form("Measured tracking efficiency versus %s and %s",xVar.Data(),yVar.Data()),
429  nxBins, TT->GetAxis(xDim)->GetBinLowEdge(1), TT->GetAxis(xDim)->GetBinUpEdge(nxBins),
430  nyBins, TT->GetAxis(yDim)->GetBinLowEdge(1), TT->GetAxis(yDim)->GetBinUpEdge(nyBins));
431  effVsXY->SetDirectory(0);
432 
433  // set the centrality and pT range for integration
434  SetCentPtCh(*TT);
435  SetCentPtCh(*TD);
436 
437  // loop over X/Y bins
438  TArrayD chEff(11);
439  TArrayD chEffErr[2];
440  chEffErr[0].Set(11);
441  chEffErr[1].Set(11);
442  for (Int_t ix = 1; ix <= nxBins; ++ix) {
443 
444  // set x range
445  TT->GetAxis(xDim)->SetRange(ix, ix);
446  TD->GetAxis(xDim)->SetRange(ix, ix);
447 
448  for (Int_t iy = 1; iy <= nyBins; ++iy) {
449 
450  // set y range
451  TT->GetAxis(yDim)->SetRange(iy, iy);
452  TD->GetAxis(yDim)->SetRange(iy, iy);
453 
454  // compute chamber efficency and errors
455  if (GetChamberEfficiency(*TT, *TD, chEff, chEffErr, kFALSE)) {
456 
457  // compute the overall tracking efficiency
458  TGraphAsymmErrors *dummy[3] = {0x0, 0x0, 0x0};
459  GetTrackingEfficiency(chEff, chEffErr, dummy, dummy, 0, 0.);
460 
461  // fill histo
462  effVsXY->Fill(TT->GetAxis(xDim)->GetBinCenter(ix),TT->GetAxis(yDim)->GetBinCenter(iy),chEff[0]);
463  effVsXY->SetBinError(ix,iy,TMath::Max(chEffErr[0][0], chEffErr[1][0]));
464 
465  } else {
466 
467  // fill histo with 0 ± 1
468  effVsXY->Fill(TT->GetAxis(xDim)->GetBinCenter(ix),TT->GetAxis(yDim)->GetBinCenter(iy),0.);
469  effVsXY->SetBinError(ix,iy,1.);
470 
471  }
472 
473  }
474 
475  }
476 
477  // close input file
478  file->Close();
479 
480  // display
481  effVsXY->GetXaxis()->SetTitle(xVar.Data());
482  effVsXY->GetXaxis()->CenterTitle(kTRUE);
483  effVsXY->GetXaxis()->SetLabelFont(22);
484  effVsXY->GetXaxis()->SetTitleFont(22);
485  effVsXY->GetYaxis()->SetTitle(yVar.Data());
486  effVsXY->GetYaxis()->CenterTitle(kTRUE);
487  effVsXY->GetYaxis()->SetLabelFont(22);
488  effVsXY->GetYaxis()->SetTitleFont(22);
489  effVsXY->GetZaxis()->SetTitle("efficiency");
490  effVsXY->GetZaxis()->SetLabelFont(22);
491  effVsXY->GetZaxis()->SetTitleFont(22);
492 
493  if (draw) {
494  new TCanvas(Form("cTrackingEffVs%s-%s",xVar.Data(),yVar.Data()), Form("Measured tracking efficiency versus %s and %s",xVar.Data(),yVar.Data()),700,600);
495  effVsXY->DrawClone("surf1");
496  }
497 
498  // save output
499  file = new TFile(fileNameSave.Data(),"update");
500  effVsXY->Write(0x0, TObject::kOverwrite);
501 
502  // add an histo with variable size rapidity bins
503  if (yDim == 3 && rap) {
504 
505  TH2F* effVsXYrap = new TH2F();
506  TString rapName = Form("trackingEffVs%s-%sRapBins", xVar.Data(), yVar.Data());
507  TString rapTitle = Form("Measured tracking efficiency versus %s and %s", xVar.Data(), yVar.Data());
508  effVsXYrap->SetTitle(rapTitle.Data());
509  effVsXYrap->SetName(rapName.Data());
510 
511  Double_t xBinEdge[nxBins+1];
512  for (Int_t xbin = 0; xbin <= nxBins; ++xbin)
513  xBinEdge[xbin] = effVsXY->GetXaxis()->GetBinLowEdge(xbin+1);
514 
515  Double_t yBinEdge[nyBins+1];
516  for (Int_t ybin = 0; ybin <= nyBins; ++ybin)
517  yBinEdge[ybin] = 2*TMath::ATan(TMath::Exp((effVsXY->GetYaxis()->GetBinLowEdge(ybin+1))));
518 
519  effVsXYrap->SetBins(nxBins, xBinEdge, nyBins, yBinEdge);
520 
521  for (Int_t xbin = 1; xbin <= nxBins; ++xbin)
522  for (Int_t ybin = 1; ybin <= nyBins; ++ybin)
523  effVsXYrap->SetBinContent(xbin, ybin, effVsXY->GetBinContent(xbin,ybin));
524 
525  effVsXYrap->Write(0x0, TObject::kOverwrite);
526 
527  delete effVsXYrap;
528 
529  }
530 
531  file->Close();
532 
533  // clean memory
534  delete effVsXY;
535 
536 }
537 
538 
539 //---------------------------------------------------------------------------
540 void PlotMuonEfficiency(TString fileNameData, TString fileNameSave, Bool_t saveEdges, Bool_t print, Bool_t draw)
541 {
543 
544  printf("plotting efficiency...\n");
545 
546  // get input hists
547  TFile *file = new TFile(fileNameData.Data(), "read");
548  if (!file || !file->IsOpen()) {
549  printf("cannot open file %s \n",fileNameData.Data());
550  return;
551  }
552  TList *listTT = static_cast<TList*>(file->FindObjectAny(Form("TotalTracksPerChamber%s", gObjNameExtension.Data())));
553  TList *listTD = static_cast<TList*>(file->FindObjectAny(Form("TracksDetectedPerChamber%s", gObjNameExtension.Data())));
554  THnSparse *TT = static_cast<THnSparse*>(listTT->At(10));
555  THnSparse *TD = static_cast<THnSparse*>(listTD->At(10));
556 
557  // output graphs
558  TGraphAsymmErrors *effVsCh[3] = {0x0, 0x0, 0x0};
559  TGraphAsymmErrors *effVsSt[3] = {0x0, 0x0, 0x0};
560  TString nameAdd[3] = {"", "Low", "Up"};
561  TString titleAdd[3] = {"", " - lower limit", " - upper limit"};
562  for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i) {
563  effVsCh[i] = CreateGraph(Form("chamberEff%s",nameAdd[i].Data()),
564  Form("Measured efficiency per chamber (0 = spectro)%s",titleAdd[i].Data()));
565  effVsSt[i] = CreateGraph(Form("stationEff%s",nameAdd[i].Data()),
566  Form("Measured efficiency per station (6 = st4&5)%s",titleAdd[i].Data()));
567  }
568 
569  // set the centrality and pT ranges for integration
570  SetCentPtCh(*TT);
571  SetCentPtCh(*TD);
572 
573  TArrayD chEff(11);
574  TArrayD chEffErr[2];
575  chEffErr[0].Set(11);
576  chEffErr[1].Set(11);
577 
578  // compute chamber efficency and errors
579  if (GetChamberEfficiency(*TT, *TD, chEff, chEffErr, print)) {
580 
581  // compute the overall tracking efficiency
582  GetTrackingEfficiency(chEff, chEffErr, effVsSt, effVsCh, 0, 0., print);
583 
584  } else {
585 
586  // set tracking efficiency to 1 +0 -1
587  effVsCh[0]->SetPoint(0,0.,1.);
588  effVsCh[0]->SetPointError(0,0.,0.,1.,0.);
589 
590  for (Int_t iSt = 0; iSt < 6; ++iSt) {
591  effVsSt[0]->SetPoint(iSt,iSt+1,1.);
592  effVsSt[0]->SetPointError(iSt,0.,0.,1.,0.);
593  }
594 
595  if (saveEdges) {
596 
597  // lower = 0 ± 0
598  effVsCh[1]->SetPoint(0,0.,0.);
599  effVsCh[1]->SetPointError(0,0.,0.,0.,0.);
600 
601  // upper = 1 ± 0
602  effVsCh[2]->SetPoint(0,0.,1.);
603  effVsCh[2]->SetPointError(0,0.,0.,0.,0.);
604 
605  for (Int_t iSt = 0; iSt < 6; ++iSt) {
606 
607  effVsSt[1]->SetPoint(iSt,iSt+1,0.);
608  effVsSt[1]->SetPointError(iSt,0.,0.,0.,0.);
609 
610  effVsSt[2]->SetPoint(iSt,iSt+1,1.);
611  effVsSt[2]->SetPointError(iSt,0.,0.,0.,0.);
612 
613  }
614 
615  }
616 
617  }
618 
619  // fill graph vs chamber
620  for (Int_t iCh = 1; iCh < 11; iCh++) {
621  for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i) {
622  effVsCh[i]->SetPoint(iCh,iCh,chEff[iCh]);
623  effVsCh[i]->SetPointError(iCh,0.,0.,chEffErr[0][iCh],chEffErr[1][iCh]);
624  }
625  }
626 
627  // close input file
628  file->Close();
629 
630  // display
631  for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i) {
632  effVsCh[i]->GetXaxis()->Set(12, -0.5, 10.5);
633  effVsCh[i]->GetXaxis()->SetNdivisions(11);
634  BeautifyGraph(*effVsCh[i], "chamber", "efficiency");
635  effVsSt[i]->GetXaxis()->Set(7, 0.5, 6.5);
636  effVsSt[i]->GetXaxis()->SetNdivisions(6);
637  BeautifyGraph(*effVsSt[i], "station", "efficiency");
638  }
639 
640  if (draw) {
641  TCanvas *c = new TCanvas("cEfficiency", "Measured tracking efficiency" , 1000, 400);
642  c->Divide(2,1);
643  gROOT->SetSelectedPad(c->cd(1));
644  effVsCh[0]->DrawClone("ap");
645  gROOT->SetSelectedPad(c->cd(2));
646  effVsSt[0]->DrawClone("ap");
647  }
648 
649  // save output
650  file = new TFile(fileNameSave.Data(),"update");
651  for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i) effVsCh[i]->Write(0x0, TObject::kOverwrite);
652  for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i) effVsSt[i]->Write(0x0, TObject::kOverwrite);
653  file->Close();
654 
655  // clean memory
656  for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i) {
657  delete effVsCh[i];
658  delete effVsSt[i];
659  }
660 
661 }
662 
663 
664 //---------------------------------------------------------------------------
665 void PlotMuonEfficiencyVsRun(TString runList, TString fileNameData, TString fileNameSave, Bool_t print, Bool_t draw)
666 {
668 
669  printf("plotting efficiency versus run...\n");
670 
671  // open run list
672  ifstream inFile(runList.Data());
673  if (!inFile.is_open()) {
674  printf("cannot open file %s\n",runList.Data());
675  return;
676  }
677 
678  // output graphs
679  TObjArray chamberVsRunGraphs; // 10 graphs: 1 for each chamber
680  chamberVsRunGraphs.SetOwner(kTRUE);
681  for ( Int_t iCh = 1; iCh < 11; ++iCh)
682  chamberVsRunGraphs.Add(CreateGraph("effCh%dVsRun", "Measured efficiency for chamber %d versus run",iCh));
683 
684  TObjArray stationVsRunGraphs[3]; // 6 graphs: 1 for each station, and 1 for station 4&5
685  TGraphAsymmErrors *trkVsRun[3];
686  TString nameAdd[3] = {"", "Low", "Up"};
687  TString titleAdd[3] = {"", " - lower limit", " - upper limit"};
688  for (Int_t i = 0; i < 3; ++i) {
689 
690  stationVsRunGraphs[i].SetOwner(kTRUE);
691  for ( Int_t iSt = 1; iSt < 6; ++iSt)
692  stationVsRunGraphs[i].Add(CreateGraph(Form("effSt%%dVsRun%s",nameAdd[i].Data()),
693  Form("Measured efficiency for station %%d versus run%s",titleAdd[i].Data()),iSt));
694  stationVsRunGraphs[i].Add(CreateGraph(Form("effSt4&5VsRun%s",nameAdd[i].Data()),
695  Form("Measured efficiency for station 4&5 versus run%s",titleAdd[i].Data())));
696 
697  trkVsRun[i] = CreateGraph(Form("trackingEffVsRun%s",nameAdd[i].Data()),
698  Form("Measured tracking efficiency versus run%s", titleAdd[i].Data()));
699 
700  }
701 
702  Int_t irun = -1;
703  TList runs;
704  runs.SetOwner();
705 
706  // loop over runs
707  while (!inFile.eof()) {
708 
709  // get current run number
710  TString currRun;
711  currRun.ReadLine(inFile,kTRUE);
712  if(currRun.IsNull()) continue;
713  runs.AddLast(new TObjString(currRun));
714  irun++;
715 
716  Int_t run = currRun.Atoi();
717 
718  printf("run %d: ", run);
719 
720  // compute efficiencies for this run
721  TString dataFile = Form("runs/%d/%s", run, fileNameData.Data());
722  TString outFile = Form("runs/%d/%s", run, fileNameSave.Data());
723  PlotMuonEfficiency(dataFile, outFile, kTRUE, print, kFALSE);
724 
725  TFile *file = new TFile(outFile.Data(), "read");
726  if (!file || !file->IsOpen()) {
727  printf("cannot open file\n");
728  continue;
729  }
730 
731  // loop over central value and edges
732  for (Int_t i = 0; i < 3; ++i) {
733 
734  // get results
735  TGraphAsymmErrors *effVsCh;
736  TGraphAsymmErrors *effVsSt;
737  effVsCh = static_cast<TGraphAsymmErrors*>(file->FindObjectAny(Form("chamberEff%s",nameAdd[i].Data())));
738  effVsSt = static_cast<TGraphAsymmErrors*>(file->FindObjectAny(Form("stationEff%s",nameAdd[i].Data())));
739  if (!effVsCh || !effVsSt) {
740  printf("efficiency graph not found\n");
741  continue;
742  }
743 
744  // fill chamber efficiency
745  if (i == 0) for ( Int_t iCh = 0; iCh < 10; ++iCh) {
746  TGraphAsymmErrors *g = static_cast<TGraphAsymmErrors*>(chamberVsRunGraphs.UncheckedAt(iCh));
747  g->SetPoint(irun,irun,effVsCh->GetY()[iCh+1]);
748  g->SetPointError(irun,0.,0.,effVsCh->GetErrorYlow(iCh+1),effVsCh->GetErrorYhigh(iCh+1));
749  }
750 
751  // fill station efficiency
752  for ( Int_t iSt = 0; iSt < 6; ++iSt) {
753  TGraphAsymmErrors *g = static_cast<TGraphAsymmErrors*>(stationVsRunGraphs[i].UncheckedAt(iSt));
754  g->SetPoint(irun,irun,effVsSt->GetY()[iSt]);
755  g->SetPointError(irun,0.,0.,effVsSt->GetErrorYlow(iSt),effVsSt->GetErrorYhigh(iSt));
756  }
757 
758  // fill spectrometer efficiency
759  trkVsRun[i]->SetPoint(irun,irun,effVsCh->GetY()[0]);
760  trkVsRun[i]->SetPointError(irun,0.,0.,effVsCh->GetErrorYlow(0),effVsCh->GetErrorYhigh(0));
761 
762  }
763 
764  file->Close();
765 
766  }
767 
768  inFile.close();
769 
770  // display
771  BeautifyGraphs(chamberVsRunGraphs,"run number","efficiency");
772  SetRunLabel(chamberVsRunGraphs,irun,runs);
773  for (Int_t i = 0; i < 3; ++i) {
774  BeautifyGraph(*trkVsRun[i],"run number","efficiency");
775  SetRunLabel(*trkVsRun[i],irun,runs);
776  BeautifyGraphs(stationVsRunGraphs[i],"run number","efficiency");
777  SetRunLabel(stationVsRunGraphs[i],irun,runs);
778  }
779 
780  if (draw) {
781  new TCanvas("cTrackingEffVsRun", "Tracking efficiency versus run",1000,400);
782  trkVsRun[0]->DrawClone("ap");
783  }
784 
785  // save output
786  TFile* file = new TFile(fileNameSave.Data(),"update");
787  chamberVsRunGraphs.Write("ChamberEffVsRun", TObject::kOverwrite | TObject::kSingleKey);
788  for (Int_t i = 0; i < 3; ++i)
789  stationVsRunGraphs[i].Write(Form("StationEffVsRun%s",nameAdd[i].Data()), TObject::kOverwrite | TObject::kSingleKey);
790  for (Int_t i = 0; i < 3; ++i) trkVsRun[i]->Write(0x0, TObject::kOverwrite);
791  file->Close();
792 
793  // clean memory
794  for (Int_t i = 0; i < 3; ++i) delete trkVsRun[i];
795 
796 }
797 
798 
799 //---------------------------------------------------------------------------
800 void PlotIntegratedMuonEfficiency(TString fileNameWeights, TString fileNameSave, Bool_t print, Bool_t draw)
801 {
803 
804  printf("plotting integrated efficiency...\n");
805 
806  // load run weights
807  if (!fileNameWeights.IsNull()) LoadRunWeights(fileNameWeights);
808  if (!runWeights) {
809  printf("Cannot compute integrated efficiency without run-by-run weights\n");
810  return;
811  }
812 
813  // get input hists
814  TFile *file = new TFile(fileNameSave.Data(), "update");
815  if (!file || !file->IsOpen()) {
816  printf("cannot open file\n");
817  return;
818  }
819  TObjArray *chamberVsRunGraphs = static_cast<TObjArray*>(file->FindObjectAny("ChamberEffVsRun"));
820  TObjArray *stationVsRunGraphs[2];
821  stationVsRunGraphs[0] = static_cast<TObjArray*>(file->FindObjectAny("StationEffVsRunLow"));
822  stationVsRunGraphs[1] = static_cast<TObjArray*>(file->FindObjectAny("StationEffVsRunUp"));
823  TGraphAsymmErrors *trkVsRun[2];
824  trkVsRun[0] = static_cast<TGraphAsymmErrors*>(file->FindObjectAny("trackingEffVsRunLow"));
825  trkVsRun[1] = static_cast<TGraphAsymmErrors*>(file->FindObjectAny("trackingEffVsRunUp"));
826  if (!chamberVsRunGraphs || !stationVsRunGraphs[0] || !stationVsRunGraphs[1] || !trkVsRun[0] || !trkVsRun[1]) {
827  printf("object not found --> you must first plot the efficiency versus run\n");
828  return;
829  }
830 
831  // output graphs
832  TGraphAsymmErrors *effVsCh = CreateGraph("integratedChamberEff", "Integrated efficiency per chamber (0 = spectro)");
833  TGraphAsymmErrors *effVsSt = CreateGraph("integratedStationEff", "Integrated efficiency per station (6 = st4&5)");
834 
835  // integrate spectrometer efficiency
836  IntegrateMuonEfficiency(*trkVsRun[0], *trkVsRun[1], *effVsCh, 0, 0.);
837 
838  // integrate chamber efficiency
839  for ( Int_t iCh = 0; iCh < 10; ++iCh) {
840  TGraphAsymmErrors *g = static_cast<TGraphAsymmErrors*>(chamberVsRunGraphs->UncheckedAt(iCh));
841  IntegrateMuonEfficiency(*g, *g, *effVsCh, iCh+1, iCh+1);
842  }
843 
844  // integrate station efficiency
845  for ( Int_t iSt = 0; iSt < 6; ++iSt) {
846  TGraphAsymmErrors *gLow = static_cast<TGraphAsymmErrors*>(stationVsRunGraphs[0]->UncheckedAt(iSt));
847  TGraphAsymmErrors *gUp = static_cast<TGraphAsymmErrors*>(stationVsRunGraphs[1]->UncheckedAt(iSt));
848  IntegrateMuonEfficiency(*gLow, *gUp, *effVsSt, iSt, iSt+1);
849  }
850 
851  // print results
852  if (print) {
853  for (Int_t iCh = 1; iCh < 11; ++iCh) {
854  cout << "Efficiency chamber " << iCh << " : ";
855  cout << effVsCh->GetY()[iCh] << " + " << effVsCh->GetErrorYhigh(iCh) << " - " << effVsCh->GetErrorYlow(iCh) << endl;
856  }
857  for (Int_t iSt = 0; iSt < 6; ++iSt) {
858  if (iSt < 5) cout << "Station " << iSt+1 << " = ";
859  else cout << "Station 45 = ";
860  cout << effVsSt->GetY()[iSt] << " + " << effVsSt->GetErrorYhigh(iSt) << " - " << effVsSt->GetErrorYlow(iSt) << endl;
861  }
862  cout << "Total tracking efficiency : ";
863  cout << effVsCh->GetY()[0] << " + " << effVsCh->GetErrorYhigh(0) << " - " << effVsCh->GetErrorYlow(0) << endl << endl;
864  }
865 
866  // display
867  effVsCh->GetXaxis()->Set(12, -0.5, 10.5);
868  effVsCh->GetXaxis()->SetNdivisions(11);
869  BeautifyGraph(*effVsCh, "chamber", "efficiency");
870  effVsSt->GetXaxis()->Set(7, 0.5, 6.5);
871  effVsSt->GetXaxis()->SetNdivisions(6);
872  BeautifyGraph(*effVsSt, "station", "efficiency");
873 
874  if (draw) {
875  TCanvas *c = new TCanvas("cIntegratedEfficiency", "Integrated tracking efficiency" , 1000, 400);
876  c->Divide(2,1);
877  gROOT->SetSelectedPad(c->cd(1));
878  effVsCh->DrawClone("ap");
879  gROOT->SetSelectedPad(c->cd(2));
880  effVsSt->DrawClone("ap");
881  }
882 
883  // save output
884  effVsCh->Write(0x0, TObject::kOverwrite);
885  effVsSt->Write(0x0, TObject::kOverwrite);
886  file->Close();
887 
888  // clean memory
889  delete effVsCh;
890  delete effVsSt;
891 
892 }
893 
894 
895 //---------------------------------------------------------------------------
896 void PlotMuonEfficiencyPerDE(TString fileNameData, TString fileNameSave, Bool_t saveEdges)
897 {
899 
900  printf("plotting efficiency per DE...\n");
901 
902  // get input hists
903  TFile *file = new TFile(fileNameData.Data(), "read");
904  if (!file || !file->IsOpen()) {
905  printf("cannot open file %s \n",fileNameData.Data());
906  return;
907  }
908  TList *listTT = static_cast<TList*>(file->FindObjectAny(Form("TotalTracksPerChamber%s", gObjNameExtension.Data())));
909  TList *listTD = static_cast<TList*>(file->FindObjectAny(Form("TracksDetectedPerChamber%s", gObjNameExtension.Data())));
910 
911  // output graph arrays
912  TObjArray chamberVsDEGraphs; // 10 graphs: 1 for each chamber
913  chamberVsDEGraphs.SetOwner(kTRUE);
914  TGraphAsymmErrors *gCh[10]; // shortcut
915 
916  TObjArray stationVsDEGraphs[3]; // 6 graphs: 1 for each station, and 1 for station 4&5
917  TString nameAdd[3] = {"", "Low", "Up"};
918  TString titleAdd[3] = {"", " - lower limit", " - upper limit"};
919  for (Int_t i = 0; i < 3; ++i) stationVsDEGraphs[i].SetOwner(kTRUE);
920  TGraphAsymmErrors *gSt[6][3]; // shortcut
921  for (Int_t iSt = 0; iSt < 6; ++iSt) for (Int_t i = 0; i < 3; ++i) gSt[iSt][i] = 0x0;
922 
923  // loop over chambers
924  for (Int_t iCh = 0; iCh < 10; ++iCh) {
925 
926  // output graphs
927  chamberVsDEGraphs.Add(CreateGraph("effCh%dVsDE","Measured efficiency for chamber %d per DE",iCh+1));
928  gCh[iCh] = static_cast<TGraphAsymmErrors*>(chamberVsDEGraphs.UncheckedAt(iCh));
929 
930  // get input hists
931  THnSparse *TT = static_cast<THnSparse*>(listTT->At(iCh));
932  THnSparse *TD = static_cast<THnSparse*>(listTD->At(iCh));
933 
934  // set the centrality and pT range for integration
935  SetCentPtCh(*TT);
936  SetCentPtCh(*TD);
937 
938  // compute DE efficency and errors
939  GetDEEfficiency(*TT, *TD, *gCh[iCh]);
940 
941  }
942 
943  // close input file
944  file->Close();
945 
946  TArrayD chEff(11);
947  TArrayD chEffErr[2];
948  chEffErr[0].Set(11);
949  chEffErr[1].Set(11);
950 
951  // loop over the first 3 stations
952  for (Int_t iSt = 0; iSt < 3; ++iSt) {
953 
954  // output graphs
955  for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i) {
956  stationVsDEGraphs[i].Add(CreateGraph(Form("effSt%%dVsDE%s",nameAdd[i].Data()),
957  Form("Measured efficiency for station %%d per DE%s",titleAdd[i].Data()),iSt+1));
958  gSt[iSt][i] = static_cast<TGraphAsymmErrors*>(stationVsDEGraphs[i].UncheckedAt(iSt));
959  }
960 
961  // loop over DE
962  Int_t nDE = gCh[2*iSt]->GetN();
963  for (Int_t iDE = 0; iDE < nDE; ++iDE) {
964 
965  // copy efficiency
966  for (Int_t iCh = 0; iCh < 2; ++iCh) {
967  chEff[2*iSt+iCh+1] = gCh[2*iSt+iCh]->GetY()[iDE];
968  chEffErr[0][2*iSt+iCh+1] = gCh[2*iSt+iCh]->GetErrorYlow(iDE);
969  chEffErr[1][2*iSt+iCh+1] = gCh[2*iSt+iCh]->GetErrorYhigh(iDE);
970  }
971 
972  // compute station efficiency
973  GetStationEfficiency(chEff, chEffErr, iSt, gSt[iSt], iDE, iDE);
974 
975  }
976 
977  }
978 
979  // output graphs for last 2 stations
980  for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i) {
981  for (Int_t iSt = 3; iSt < 5; ++iSt) {
982  stationVsDEGraphs[i].Add(CreateGraph(Form("effSt%%dVsDE%s",nameAdd[i].Data()),
983  Form("Measured efficiency for station %%d per DE%s",titleAdd[i].Data()),iSt+1));
984  gSt[iSt][i] = static_cast<TGraphAsymmErrors*>(stationVsDEGraphs[i].UncheckedAt(iSt));
985  }
986  stationVsDEGraphs[i].Add(CreateGraph(Form("effSt4&5VsDE%s",nameAdd[i].Data()),
987  Form("Measured efficiency for station 4&5 per DE%s",titleAdd[i].Data())));
988  gSt[5][i] = static_cast<TGraphAsymmErrors*>(stationVsDEGraphs[i].UncheckedAt(5));
989  }
990 
991  // loop over DE
992  Int_t nDE = gCh[6]->GetN();
993  for (Int_t iDE = 0; iDE < nDE; ++iDE) {
994 
995  // copy efficiency
996  for (Int_t iCh = 6; iCh < 10; ++iCh) {
997  chEff[iCh+1] = gCh[iCh]->GetY()[iDE];
998  chEffErr[0][iCh+1] = gCh[iCh]->GetErrorYlow(iDE);
999  chEffErr[1][iCh+1] = gCh[iCh]->GetErrorYhigh(iDE);
1000  }
1001 
1002  // compute station 4&5 efficiency individually
1003  for (Int_t iSt = 3; iSt < 5; ++iSt) GetStationEfficiency(chEff, chEffErr, iSt, gSt[iSt], iDE, iDE);
1004 
1005  // compute station 4&5 efficiency together
1006  GetStation45Efficiency(chEff, chEffErr, gSt[5], iDE, iDE);
1007 
1008  }
1009 
1010  // display
1011  for (Int_t iCh = 0; iCh < 10; ++iCh) {
1012  Int_t nDE = gCh[iCh]->GetN();
1013  gCh[iCh]->GetXaxis()->Set(nDE+1, -0.5, nDE-0.5);
1014  gCh[iCh]->GetXaxis()->SetNdivisions(nDE);
1015  }
1016  BeautifyGraphs(chamberVsDEGraphs,"Detection Element","efficiency");
1017  for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i) {
1018  for (Int_t iSt = 0; iSt < 6; ++iSt) {
1019  Int_t nDE = gSt[iSt][i]->GetN();
1020  gSt[iSt][i]->GetXaxis()->Set(nDE+1, -0.5, nDE-0.5);
1021  gSt[iSt][i]->GetXaxis()->SetNdivisions(nDE);
1022  }
1023  BeautifyGraphs(stationVsDEGraphs[i],"Detection Element","efficiency");
1024 
1025  }
1026 
1027  // save Output
1028  file = new TFile(fileNameSave.Data(),"update");
1029  chamberVsDEGraphs.Write("ChamberEffVsDE", TObject::kOverwrite | TObject::kSingleKey);
1030  for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i)
1031  stationVsDEGraphs[i].Write(Form("StationEffVsDE%s",nameAdd[i].Data()), TObject::kOverwrite | TObject::kSingleKey);
1032  file->Close();
1033 
1034 }
1035 
1036 
1037 //---------------------------------------------------------------------------
1038 void PlotMuonEfficiencyPerDEVsRun(TString runList, TString fileNameData, TString fileNameSave)
1039 {
1041 
1042  printf("plotting efficiency per DE versus run...\n");
1043 
1044  // open run list
1045  ifstream inFile(runList.Data());
1046  if (!inFile.is_open()) {
1047  printf("cannot open file %s\n",runList.Data());
1048  return;
1049  }
1050 
1051  // output graphs
1052  TObjArray deVsRunGraphs; // 1 graph per DE
1053  deVsRunGraphs.SetOwner(kTRUE);
1054  TObjArray stDEVsRunGraphs[3]; // 1 graph per pair (quartet) of DE in individual stations (stations 4&5)
1055  TString nameAdd[3] = {"", "Low", "Up"};
1056  TString titleAdd[3] = {"", " - lower limit", " - upper limit"};
1057  for (Int_t i = 0; i < 3; ++i) stDEVsRunGraphs[i].SetOwner(kTRUE);
1058  Bool_t createGraph = kTRUE;
1059 
1060  Int_t irun = -1;
1061  TList runs;
1062  runs.SetOwner();
1063 
1064  // loop over runs
1065  while (!inFile.eof()) {
1066 
1067  // get current run number
1068  TString currRun;
1069  currRun.ReadLine(inFile,kTRUE);
1070  if(currRun.IsNull()) continue;
1071  runs.AddLast(new TObjString(currRun));
1072  irun++;
1073 
1074  Int_t run = currRun.Atoi();
1075 
1076  printf("run %d: ", run);
1077 
1078  // compute efficiencies for this run
1079  TString dataFile = Form("runs/%d/%s", run, fileNameData.Data());
1080  TString outFile = Form("runs/%d/%s", run, fileNameSave.Data());
1081  PlotMuonEfficiencyPerDE(dataFile, outFile, kTRUE);
1082 
1083  TFile *file = new TFile(outFile.Data(), "read");
1084  if (!file || !file->IsOpen()) {
1085  printf("cannot open file\n");
1086  continue;
1087  }
1088 
1089  // get results
1090  TObjArray *chamberVsDEGraphs = static_cast<TObjArray*>(file->FindObjectAny("ChamberEffVsDE"));
1091  if (!chamberVsDEGraphs) {
1092  printf("efficiency graph not found\n");
1093  continue;
1094  }
1095 
1096  Int_t currentDE = 0;
1097 
1098  // loop over chambers
1099  for ( Int_t iCh = 0; iCh < 10; ++iCh) {
1100 
1101  // input graph
1102  TGraphAsymmErrors *g = static_cast<TGraphAsymmErrors*>(chamberVsDEGraphs->UncheckedAt(iCh));
1103  Int_t nDE = g->GetN();
1104 
1105  // loop over DE
1106  for (Int_t iDE = 0; iDE < nDE; ++iDE) {
1107 
1108  // output graph
1109  if (createGraph) deVsRunGraphs.Add(CreateGraph("effDE%dVsRun","Measured efficiency for DE %d versus run",100*(iCh+1)+iDE));
1110  TGraphAsymmErrors *gDE= static_cast<TGraphAsymmErrors*>(deVsRunGraphs.UncheckedAt(currentDE++));
1111 
1112  // fill DE efficiency
1113  gDE->SetPoint(irun,irun,g->GetY()[iDE]);
1114  gDE->SetPointError(irun,0.,0.,g->GetErrorYlow(iDE),g->GetErrorYhigh(iDE));
1115 
1116  }
1117 
1118  }
1119 
1120  // loop over central value and edges
1121  for (Int_t i = 0; i < 3; ++i) {
1122 
1123  // get results
1124  TObjArray *stationVsDEGraphs = static_cast<TObjArray*>(file->FindObjectAny(Form("StationEffVsDE%s",nameAdd[i].Data())));
1125  if (!stationVsDEGraphs) {
1126  printf("efficiency graph not found\n");
1127  continue;
1128  }
1129 
1130  Int_t currentStDE = 0;
1131 
1132  // loop over stations
1133  for ( Int_t iSt = 0; iSt < 6; ++iSt) {
1134 
1135  // input graph
1136  TGraphAsymmErrors *g = static_cast<TGraphAsymmErrors*>(stationVsDEGraphs->UncheckedAt(iSt));
1137  Int_t nDE = g->GetN();
1138 
1139  // loop over DE
1140  for (Int_t iDE = 0; iDE < nDE; ++iDE) {
1141 
1142  // output graph
1143  if (createGraph) {
1144  TString sSt = (iSt<5) ? Form("%d",iSt+1) : "4&5";
1145  stDEVsRunGraphs[i].Add(CreateGraph(Form("effSt%sDE%%dVsRun%s",sSt.Data(),nameAdd[i].Data()),
1146  Form("Measured efficiency for DE %%d in station %s versus run%s",sSt.Data(),titleAdd[i].Data()),iDE));
1147  }
1148  TGraphAsymmErrors *gDE= static_cast<TGraphAsymmErrors*>(stDEVsRunGraphs[i].UncheckedAt(currentStDE++));
1149 
1150  // fill DE efficiency
1151  gDE->SetPoint(irun,irun,g->GetY()[iDE]);
1152  gDE->SetPointError(irun,0.,0.,g->GetErrorYlow(iDE),g->GetErrorYhigh(iDE));
1153 
1154  }
1155 
1156  }
1157 
1158  }
1159 
1160  file->Close();
1161 
1162  createGraph = kFALSE;
1163 
1164  }
1165 
1166  inFile.close();
1167 
1168  // display
1169  BeautifyGraphs(deVsRunGraphs,"run number","efficiency");
1170  SetRunLabel(deVsRunGraphs,irun,runs);
1171  for (Int_t i = 0; i < 3; ++i) {
1172  BeautifyGraphs(stDEVsRunGraphs[i],"run number","efficiency");
1173  SetRunLabel(stDEVsRunGraphs[i],irun,runs);
1174  }
1175 
1176  // save output
1177  TFile* file = new TFile(fileNameSave.Data(),"update");
1178  deVsRunGraphs.Write("DEEffVsRun", TObject::kOverwrite | TObject::kSingleKey);
1179  for (Int_t i = 0; i < 3; ++i)
1180  stDEVsRunGraphs[i].Write(Form("DEEffPerStationVsRun%s",nameAdd[i].Data()), TObject::kOverwrite | TObject::kSingleKey);
1181  file->Close();
1182 
1183 }
1184 
1185 
1186 //---------------------------------------------------------------------------
1187 void PlotIntegratedMuonEfficiencyPerDE(TString fileNameWeights, TString fileNameSave)
1188 {
1190 
1191  printf("plotting integrated efficiency per DE...\n");
1192 
1193  // load run weights
1194  if (!fileNameWeights.IsNull()) LoadRunWeights(fileNameWeights);
1195  if (!runWeights) {
1196  printf("Cannot compute integrated efficiency without run-by-run weights\n");
1197  return;
1198  }
1199 
1200  // get input hists
1201  TFile *file = new TFile(fileNameSave.Data(), "update");
1202  if (!file || !file->IsOpen()) {
1203  printf("cannot open file\n");
1204  return;
1205  }
1206  TObjArray *deVsRunGraphs = static_cast<TObjArray*>(file->FindObjectAny("DEEffVsRun"));
1207  TObjArray *stDEVsRunGraphs[2];
1208  stDEVsRunGraphs[0] = static_cast<TObjArray*>(file->FindObjectAny("DEEffPerStationVsRunLow"));
1209  stDEVsRunGraphs[1] = static_cast<TObjArray*>(file->FindObjectAny("DEEffPerStationVsRunUp"));
1210  if (!deVsRunGraphs || !stDEVsRunGraphs[0] || !stDEVsRunGraphs[1]) {
1211  printf("object not found --> you must first plot the efficiency versus run\n");
1212  return;
1213  }
1214 
1215  // output graph arrays
1216  TObjArray chamberVsDEGraphs; // 10 graphs: 1 for each chamber
1217  chamberVsDEGraphs.SetOwner(kTRUE);
1218  for ( Int_t iCh = 0; iCh < 10; ++iCh)
1219  chamberVsDEGraphs.Add(CreateGraph("integratedEffCh%dVsDE","Integrated efficiency for chamber %d per DE",iCh+1));
1220 
1221  TObjArray stationVsDEGraphs; // 6 graphs: 1 for each station, and 1 for station 4&5
1222  stationVsDEGraphs.SetOwner(kTRUE);
1223  for ( Int_t iSt = 0; iSt < 5; ++iSt)
1224  stationVsDEGraphs.Add(CreateGraph("integratedEffSt%dVsDE","Integrated efficiency for station %d per DE",iSt+1));
1225  stationVsDEGraphs.Add(CreateGraph("integratedEffSt4&5VsDE","Integrated efficiency for station 4&5 per DE"));
1226 
1227  // Loop over DE
1228  TIter nextDE(deVsRunGraphs);
1229  TGraphAsymmErrors *gDE = 0x0;
1230  while ((gDE = static_cast<TGraphAsymmErrors*>(nextDE()))) {
1231 
1232  // get chamber and DE indices
1233  Int_t deId;
1234  sscanf(gDE->GetName(), "effDE%dVsRun", &deId);
1235  Int_t iCh = deId/100-1;
1236  Int_t iDE = deId%100;
1237 
1238  // integrate DE efficiency
1239  TGraphAsymmErrors *g = static_cast<TGraphAsymmErrors*>(chamberVsDEGraphs.UncheckedAt(iCh));
1240  IntegrateMuonEfficiency(*gDE, *gDE, *g, iDE, iDE);
1241 
1242  }
1243 
1244  // Loop over DE per station
1245  Int_t ng = stDEVsRunGraphs[0]->GetEntries();
1246  for (Int_t ig = 0; ig < ng; ++ig) {
1247 
1248  TGraphAsymmErrors *gDELow = static_cast<TGraphAsymmErrors*>(stDEVsRunGraphs[0]->UncheckedAt(ig));
1249  TGraphAsymmErrors *gDEUp = static_cast<TGraphAsymmErrors*>(stDEVsRunGraphs[1]->UncheckedAt(ig));
1250 
1251  // get station and DE indices
1252  Int_t iSt, iDE;
1253  if (strstr(gDELow->GetName(), "4&5")) {
1254  iSt = 5;
1255  sscanf(gDELow->GetName(), "effSt4&5DE%dVsRunLow", &iDE);
1256  } else {
1257  sscanf(gDELow->GetName(), "effSt%dDE%dVsRunLow", &iSt, &iDE);
1258  iSt--;
1259  }
1260 
1261  // Integrate DE efficiency per station
1262  TGraphAsymmErrors *g = static_cast<TGraphAsymmErrors*>(stationVsDEGraphs.UncheckedAt(iSt));
1263  IntegrateMuonEfficiency(*gDELow, *gDEUp, *g, iDE, iDE);
1264 
1265  }
1266 
1267  // display
1268  for ( Int_t iCh = 0; iCh < 10; ++iCh) {
1269  TGraphAsymmErrors *g = static_cast<TGraphAsymmErrors*>(chamberVsDEGraphs.UncheckedAt(iCh));
1270  Int_t nDE = g->GetN();
1271  g->GetXaxis()->Set(nDE+1, -0.5, nDE-0.5);
1272  g->GetXaxis()->SetNdivisions(nDE);
1273  }
1274  BeautifyGraphs(chamberVsDEGraphs,"Detection Element","efficiency");
1275  for ( Int_t iSt = 0; iSt < 6; ++iSt) {
1276  TGraphAsymmErrors *g = static_cast<TGraphAsymmErrors*>(stationVsDEGraphs.UncheckedAt(iSt));
1277  Int_t nDE = g->GetN();
1278  g->GetXaxis()->Set(nDE+1, -0.5, nDE-0.5);
1279  g->GetXaxis()->SetNdivisions(nDE);
1280  }
1281  BeautifyGraphs(stationVsDEGraphs,"Detection Element","efficiency");
1282 
1283  // save Output
1284  chamberVsDEGraphs.Write("IntegratedChamberEffVsDE", TObject::kOverwrite | TObject::kSingleKey);
1285  stationVsDEGraphs.Write("IntegratedStationEffVsDE", TObject::kOverwrite | TObject::kSingleKey);
1286  file->Close();
1287 
1288 }
1289 
1290 
1291 //---------------------------------------------------------------------------
1292 Bool_t GetChamberEfficiency(THnSparse &TT, THnSparse &TD, TArrayD &chEff, TArrayD chEffErr[2], Bool_t printError)
1293 {
1296 
1297  // project track hists over the chamber axis
1298  TH1D *TTdraw = TT.Projection(0,"e");
1299  TH1D *TDdraw = TD.Projection(0,"e");
1300 
1301  // compute chamber efficiency and errors
1302  TGraphAsymmErrors *efficiency = (TTdraw->GetEntries() > 0) ? new TGraphAsymmErrors(TDdraw, TTdraw,Form("%se0",effErrMode)) : 0x0;
1303  Bool_t ok = (efficiency);
1304 
1305  // fill arrays
1306  if (ok) {
1307 
1308  Bool_t missingEff = kFALSE;
1309 
1310  for (Int_t i = 0; i < 10; i++) {
1311 
1312  if (TTdraw->GetBinContent(i+1) > 0) {
1313 
1314  chEff[i+1] = efficiency->GetY()[i];
1315  chEffErr[0][i+1] = efficiency->GetErrorYlow(i);
1316  chEffErr[1][i+1] = efficiency->GetErrorYhigh(i);
1317 
1318  } else {
1319 
1320  chEff[i+1] = -1.;
1321  chEffErr[0][i+1] = 0.;
1322  chEffErr[1][i+1] = 0.;
1323 
1324  missingEff = kTRUE;
1325 
1326  }
1327 
1328  }
1329 
1330  if (missingEff && printError) cout << "efficiency partially unknown" << endl;
1331 
1332  } else {
1333 
1334  for (Int_t i = 0; i < 10; i++) {
1335 
1336  chEff[i+1] = -1.;
1337  chEffErr[0][i+1] = 0.;
1338  chEffErr[1][i+1] = 0.;
1339 
1340  }
1341 
1342  if (printError) cout << "efficiency unknown" << endl << endl;
1343 
1344  }
1345 
1346  // clean memory
1347  delete TTdraw;
1348  delete TDdraw;
1349  delete efficiency;
1350 
1351  return ok;
1352 
1353 }
1354 
1355 
1356 //---------------------------------------------------------------------------
1357 void GetDEEfficiency(THnSparse &TT, THnSparse &TD, TGraphAsymmErrors &effVsDE)
1358 {
1360 
1361  // project track hists over the chamber axis
1362  TH1D *TTdraw = TT.Projection(0,"e");
1363  TH1D *TDdraw = TD.Projection(0,"e");
1364 
1365  // compute DE efficiency and errors
1366  TGraphAsymmErrors *efficiency = (TTdraw->GetEntries() > 0) ? new TGraphAsymmErrors(TDdraw, TTdraw,Form("%se0",effErrMode)) : 0x0;
1367  Int_t nDE = TTdraw->GetNbinsX();
1368 
1369  if (efficiency) {
1370 
1371  for (Int_t iDE = 0; iDE < nDE; ++iDE) {
1372 
1373  if (TTdraw->GetBinContent(iDE+1) > 0) {
1374 
1375  effVsDE.SetPoint(iDE,iDE,efficiency->GetY()[iDE]);
1376  effVsDE.SetPointError(iDE,0,0,efficiency->GetErrorYlow(iDE),efficiency->GetErrorYhigh(iDE));
1377 
1378  } else {
1379 
1380  effVsDE.SetPoint(iDE,iDE,-1);
1381  effVsDE.SetPointError(iDE,0,0,0,0);
1382 
1383  }
1384 
1385  }
1386 
1387  } else {
1388 
1389  for (Int_t iDE = 0; iDE < nDE; ++iDE) {
1390 
1391  effVsDE.SetPoint(iDE,iDE,-1);
1392  effVsDE.SetPointError(iDE,0,0,0,0);
1393 
1394  }
1395 
1396  }
1397 
1398  // clean memory
1399  delete TTdraw;
1400  delete TDdraw;
1401  delete efficiency;
1402 
1403 }
1404 
1405 
1406 //---------------------------------------------------------------------------
1407 void ComputeStationEfficiency(TArrayD &chEff, TArrayD chEffErr[2], Int_t iSt, Double_t &stEff, Double_t stEffErr[2])
1408 {
1410 
1411  stEff = 1.-(1.-chEff[2*iSt+1])*(1.-chEff[2*iSt+2]);
1412 
1413  Double_t d1 = (1. - chEff[2*iSt+2]); d1 *= d1;
1414  Double_t d2 = (1. - chEff[2*iSt+1]); d2 *= d2;
1415 
1416  for (Int_t i = 0; i < 2; ++i) {
1417  Double_t s1 = chEffErr[i][2*iSt+1] * chEffErr[i][2*iSt+1];
1418  Double_t s2 = chEffErr[i][2*iSt+2] * chEffErr[i][2*iSt+2];
1419  stEffErr[i] = TMath::Sqrt(d1*s1 + d2*s2 + s1*s2);
1420  }
1421 
1422  stEffErr[0] = TMath::Min(stEff, stEffErr[0]);
1423  stEffErr[1] = TMath::Min(1.-stEff, stEffErr[1]);
1424 
1425 }
1426 
1427 
1428 //---------------------------------------------------------------------------
1429 void GetStationEfficiency(TArrayD &chEff, TArrayD chEffErr[2], Int_t iSt, TGraphAsymmErrors *effVsX[3], Int_t ip, Double_t xp)
1430 {
1432 
1433  if (chEff[2*iSt+1] >= 0 && chEff[2*iSt+2] >= 0) {
1434 
1435  // compute station efficiency from known chamber efficiency
1436  Double_t stEff, stEffErr[2];
1437  ComputeStationEfficiency(chEff, chEffErr, iSt, stEff, stEffErr);
1438 
1439  // fill graphs if required
1440  for (Int_t i = 0; i < 3; ++i) {
1441  if (effVsX[i]) {
1442  effVsX[i]->SetPoint(ip,xp,stEff);
1443  effVsX[i]->SetPointError(ip,0.,0.,stEffErr[0],stEffErr[1]);
1444  }
1445  }
1446 
1447  } else {
1448 
1449  Double_t edge[2] = {0., 1.};
1450  TArrayD chEffEdge[2];
1451  Double_t stEffEdge[2];
1452  Double_t stEffEdgeErr[2][2];
1453 
1454  for (Int_t i = 0; i < 2; ++i) {
1455 
1456  // set lower(upper) limit of chamber efficiency
1457  chEffEdge[i].Set(11);
1458  for (Int_t iCh = 1; iCh < 3; ++iCh) chEffEdge[i][2*iSt+iCh] = (chEff[2*iSt+iCh] < 0) ? edge[i] : chEff[2*iSt+iCh];
1459 
1460  // compute station efficiency
1461  ComputeStationEfficiency(chEffEdge[i], chEffErr, iSt, stEffEdge[i], stEffEdgeErr[i]);
1462 
1463  // fill graph if required
1464  if (effVsX[i+1]) {
1465  effVsX[i+1]->SetPoint(ip,xp,stEffEdge[i]);
1466  effVsX[i+1]->SetPointError(ip,0.,0.,stEffEdgeErr[i][0],stEffEdgeErr[i][1]);
1467  }
1468 
1469  }
1470 
1471  // combine extreme cases to get station efficiency and fill graph if required
1472  if (effVsX[0]) {
1473  effVsX[0]->SetPoint(ip,xp,stEffEdge[1]);
1474  effVsX[0]->SetPointError(ip,0.,0.,stEffEdge[1]-stEffEdge[0]+stEffEdgeErr[0][0],stEffEdgeErr[1][1]);
1475  }
1476 
1477  }
1478 
1479 }
1480 
1481 
1482 //---------------------------------------------------------------------------
1483 void ComputeStation45Efficiency(TArrayD &chEff, TArrayD chEffErr[2], Double_t &st45Eff, Double_t st45EffErr[2])
1484 {
1486 
1487  st45Eff = chEff[7]*chEff[8]*chEff[9] + chEff[7]*chEff[8]*chEff[10] + chEff[7]*chEff[9]*chEff[10] + chEff[8]*chEff[9]*chEff[10] - 3.*chEff[7]*chEff[8]*chEff[9]*chEff[10];
1488 
1489  Double_t d1 = chEff[8]*chEff[9] + chEff[8]*chEff[10] + chEff[9]*chEff[10] - 3.*chEff[8]*chEff[9]*chEff[10]; d1 *= d1;
1490  Double_t d2 = chEff[7]*chEff[9] + chEff[7]*chEff[10] + chEff[9]*chEff[10] - 3.*chEff[7]*chEff[9]*chEff[10]; d2 *= d2;
1491  Double_t d3 = chEff[7]*chEff[8] + chEff[7]*chEff[10] + chEff[8]*chEff[10] - 3.*chEff[7]*chEff[8]*chEff[10]; d3 *= d3;
1492  Double_t d4 = chEff[7]*chEff[8] + chEff[7]*chEff[9] + chEff[8]*chEff[9] - 3.*chEff[7]*chEff[8]*chEff[9]; d4 *= d4;
1493  Double_t d12 = chEff[9] + chEff[10] - 3.*chEff[9]*chEff[10]; d12 *= d12;
1494  Double_t d13 = chEff[8] + chEff[10] - 3.*chEff[8]*chEff[10]; d13 *= d13;
1495  Double_t d14 = chEff[8] + chEff[9] - 3.*chEff[8]*chEff[9]; d14 *= d14;
1496  Double_t d23 = chEff[7] + chEff[10] - 3.*chEff[7]*chEff[10]; d23 *= d23;
1497  Double_t d24 = chEff[7] + chEff[9] - 3.*chEff[7]*chEff[9]; d24 *= d24;
1498  Double_t d34 = chEff[7] + chEff[8] - 3.*chEff[7]*chEff[8]; d34 *= d34;
1499  Double_t d123 = 1. - 3.*chEff[10]; d123 *= d123;
1500  Double_t d124 = 1. - 3.*chEff[9]; d124 *= d124;
1501  Double_t d134 = 1. - 3.*chEff[8]; d134 *= d134;
1502  Double_t d234 = 1. - 3.*chEff[7]; d234 *= d234;
1503  Double_t d1234 = - 3.; d1234 *= d1234;
1504 
1505  for (Int_t i = 0; i < 2; ++i) {
1506  Double_t s1 = chEffErr[i][7] * chEffErr[i][7];
1507  Double_t s2 = chEffErr[i][8] * chEffErr[i][8];
1508  Double_t s3 = chEffErr[i][9] * chEffErr[i][9];
1509  Double_t s4 = chEffErr[i][10] * chEffErr[i][10];
1510  st45EffErr[i] = TMath::Sqrt(d1*s1 + d2*s2 + d3*s3 + d4*s4 + d12*s1*s2 + d13*s1*s3 + d14*s1*s4 + d23*s2*s3 + d24*s2*s4 + d34*s3*s4 + d123*s1*s2*s3 + d124*s1*s2*s4 + d134*s1*s3*s4 + d234*s2*s3*s4 + d1234*s1*s2*s3*s4);
1511  }
1512 
1513  st45EffErr[0] = TMath::Min(st45Eff, st45EffErr[0]);
1514  st45EffErr[1] = TMath::Min(1.-st45Eff, st45EffErr[1]);
1515 
1516 }
1517 
1518 
1519 //---------------------------------------------------------------------------
1520 void GetStation45Efficiency(TArrayD &chEff, TArrayD chEffErr[2], TGraphAsymmErrors *effVsX[3], Int_t ip, Double_t xp)
1521 {
1523 
1524  if (chEff[7] >= 0 && chEff[8] >= 0 && chEff[9] >= 0 && chEff[10] >= 0) {
1525 
1526  // compute station efficiency from known chamber efficiency
1527  Double_t stEff, stEffErr[2];
1528  ComputeStation45Efficiency(chEff, chEffErr, stEff, stEffErr);
1529 
1530  // fill graphs if required
1531  for (Int_t i = 0; i < 3; ++i) {
1532  if (effVsX[i]) {
1533  effVsX[i]->SetPoint(ip,xp,stEff);
1534  effVsX[i]->SetPointError(ip,0.,0.,stEffErr[0],stEffErr[1]);
1535  }
1536  }
1537 
1538  } else {
1539 
1540  Double_t edge[2] = {0., 1.};
1541  TArrayD chEffEdge[2];
1542  Double_t stEffEdge[2];
1543  Double_t stEffEdgeErr[2][2];
1544 
1545  for (Int_t i = 0; i < 2; ++i) {
1546 
1547  // set lower(upper) limit of chamber efficiency
1548  chEffEdge[i].Set(11);
1549  for (Int_t iCh = 7; iCh < 11; ++iCh) chEffEdge[i][iCh] = (chEff[iCh] < 0) ? edge[i] : chEff[iCh];
1550 
1551  // compute station efficiency
1552  ComputeStation45Efficiency(chEffEdge[i], chEffErr, stEffEdge[i], stEffEdgeErr[i]);
1553 
1554  // fill graph if required
1555  if (effVsX[i+1]) {
1556  effVsX[i+1]->SetPoint(ip,xp,stEffEdge[i]);
1557  effVsX[i+1]->SetPointError(ip,0.,0.,stEffEdgeErr[i][0],stEffEdgeErr[i][1]);
1558  }
1559 
1560  }
1561 
1562  // combine extreme cases to get station efficiency and fill graph if required
1563  if (effVsX[0]) {
1564  effVsX[0]->SetPoint(ip,xp,stEffEdge[1]);
1565  effVsX[0]->SetPointError(ip,0.,0.,stEffEdge[1]-stEffEdge[0]+stEffEdgeErr[0][0],stEffEdgeErr[1][1]);
1566  }
1567 
1568  }
1569 
1570 }
1571 
1572 
1573 //---------------------------------------------------------------------------
1574 void ComputeTrackingEfficiency(Double_t stEff[6], Double_t stEffErr[6][2], Double_t &spectroEff, Double_t spectroEffErr[2])
1575 {
1577 
1578  Double_t de[6][2];
1579  for (Int_t iSt = 0; iSt < 6; iSt++) de[iSt][0] = stEff[iSt]*stEff[iSt];
1580 
1581  if (moreTrackCandidates) {
1582 
1583  spectroEff = stEff[0] * stEff[1] * stEff[2] * stEff[3] * stEff[4];
1584 
1585  for (Int_t i = 0; i < 2; i++) {
1586 
1587  for (Int_t iSt = 0; iSt < 6; iSt++) de[iSt][1] = stEffErr[iSt][i]*stEffErr[iSt][i];
1588 
1589  spectroEffErr[i] = 0.;
1590  for (Int_t j = 1; j < 32; j++) {
1591  Double_t sigmaAdd = 1.;
1592  for (Int_t iSt = 0; iSt < 5; iSt++) sigmaAdd *= de[iSt][TESTBIT(j,iSt)];
1593  spectroEffErr[i] += sigmaAdd;
1594  }
1595  spectroEffErr[i] = TMath::Sqrt(spectroEffErr[i]);
1596 
1597  }
1598 
1599  } else {
1600 
1601  spectroEff = stEff[0] * stEff[1] * stEff[2] * stEff[5];
1602 
1603  for (Int_t i = 0; i < 2; i++) {
1604 
1605  for (Int_t iSt = 0; iSt < 6; iSt++) de[iSt][1] = stEffErr[iSt][i]*stEffErr[iSt][i];
1606 
1607  spectroEffErr[i] = 0.;
1608  for (Int_t j = 1; j < 16; j++) {
1609  Double_t sigmaAdd = de[5][TESTBIT(j,3)];
1610  for (Int_t iSt = 0; iSt < 3; iSt++) sigmaAdd *= de[iSt][TESTBIT(j,iSt)];
1611  spectroEffErr[i] += sigmaAdd;
1612  }
1613  spectroEffErr[i] = TMath::Sqrt(spectroEffErr[i]);
1614 
1615  }
1616 
1617  }
1618 
1619  spectroEffErr[0] = TMath::Min(spectroEff, spectroEffErr[0]);
1620  spectroEffErr[1] = TMath::Min(1.-spectroEff, spectroEffErr[1]);
1621 
1622 }
1623 
1624 
1625 //---------------------------------------------------------------------------
1626 void GetTrackingEfficiency(TArrayD &chEff, TArrayD chEffErr[2], TGraphAsymmErrors *effVsSt[3],
1627  TGraphAsymmErrors *effVsX[3], Int_t ip, Double_t xp, Bool_t print)
1628 {
1630 
1631  Double_t stEff[6];
1632  Double_t stEffErr[6][2];
1633 
1634  // check if unknown chamber efficiency
1635  Bool_t allEffKnown = kTRUE;
1636  for (Int_t iCh = 1; iCh < 11 && allEffKnown; ++iCh) if (chEff[iCh] < 0) allEffKnown = kFALSE;
1637 
1638  if (allEffKnown) {
1639 
1640  // compute station efficiency
1641  for (Int_t iSt = 0; iSt < 5; ++iSt) ComputeStationEfficiency(chEff, chEffErr, iSt, stEff[iSt], stEffErr[iSt]);
1642  ComputeStation45Efficiency(chEff, chEffErr, stEff[5], stEffErr[5]);
1643 
1644  // compute spectrometer efficiency
1645  Double_t spectroEff, spectroEffErr[2];
1646  ComputeTrackingEfficiency(stEff, stEffErr, spectroEff, spectroEffErr);
1647  chEff[0] = spectroEff;
1648  chEffErr[0][0] = spectroEffErr[0];
1649  chEffErr[1][0] = spectroEffErr[1];
1650 
1651  // fill graphs if required
1652  for (Int_t i = 0; i < 3; ++i) {
1653  if (effVsX[i]) {
1654  effVsX[i]->SetPoint(ip,xp,chEff[0]);
1655  effVsX[i]->SetPointError(ip,0.,0.,chEffErr[0][0],chEffErr[1][0]);
1656  }
1657  if (effVsSt[i]) for (Int_t iSt = 0; iSt < 6; ++iSt) {
1658  effVsSt[i]->SetPoint(iSt,iSt+1,stEff[iSt]);
1659  effVsSt[i]->SetPointError(iSt,0.,0.,stEffErr[iSt][0],stEffErr[iSt][1]);
1660  }
1661  }
1662 
1663  } else {
1664 
1665  Double_t edge[2] = {0., 1.};
1666  TArrayD chEffEdge[2];
1667  Double_t stEffEdge[2][6];
1668  Double_t stEffEdgeErr[2][6][2];
1669  Double_t spectroEffEdge[2], spectroEffEdgeErr[2][2];
1670 
1671  for (Int_t i = 0; i < 2; ++i) {
1672 
1673  // set lower(upper) limit of chamber efficiency
1674  chEffEdge[i].Set(11);
1675  for (Int_t iCh = 1; iCh < 11; ++iCh) chEffEdge[i][iCh] = (chEff[iCh] < 0) ? edge[i] : chEff[iCh];
1676 
1677  // compute station efficiency
1678  for (Int_t iSt = 0; iSt < 5; ++iSt) ComputeStationEfficiency(chEffEdge[i], chEffErr, iSt, stEffEdge[i][iSt], stEffEdgeErr[i][iSt]);
1679  ComputeStation45Efficiency(chEffEdge[i], chEffErr, stEffEdge[i][5], stEffEdgeErr[i][5]);
1680 
1681  // compute spectrometer efficiency
1682  ComputeTrackingEfficiency(stEffEdge[i], stEffEdgeErr[i], spectroEffEdge[i], spectroEffEdgeErr[i]);
1683 
1684  // fill graph if required
1685  if (effVsX[i+1]) {
1686  effVsX[i+1]->SetPoint(ip,xp,spectroEffEdge[i]);
1687  effVsX[i+1]->SetPointError(ip,0.,0.,spectroEffEdgeErr[i][0],spectroEffEdgeErr[i][1]);
1688  }
1689  if (effVsSt[i+1]) for (Int_t iSt = 0; iSt < 6; ++iSt) {
1690  effVsSt[i+1]->SetPoint(iSt,iSt+1,stEffEdge[i][iSt]);
1691  effVsSt[i+1]->SetPointError(iSt,0.,0.,stEffEdgeErr[i][iSt][0],stEffEdgeErr[i][iSt][1]);
1692  }
1693 
1694  }
1695 
1696  // combine extreme cases to get station efficiency
1697  for (Int_t iSt = 0; iSt < 6; ++iSt) {
1698  stEff[iSt] = stEffEdge[1][iSt];
1699  stEffErr[iSt][0] = stEffEdge[1][iSt] - stEffEdge[0][iSt] + stEffEdgeErr[0][iSt][0];
1700  stEffErr[iSt][1] = stEffEdgeErr[1][iSt][1];
1701  }
1702 
1703  // combine extreme cases to get spectrometer efficiency
1704  chEff[0] = spectroEffEdge[1];
1705  chEffErr[0][0] = spectroEffEdge[1] - spectroEffEdge[0] + spectroEffEdgeErr[0][0];
1706  chEffErr[1][0] = spectroEffEdgeErr[1][1];
1707 
1708  // fill graph if required
1709  if (effVsX[0]) {
1710  effVsX[0]->SetPoint(ip,xp,chEff[0]);
1711  effVsX[0]->SetPointError(ip,0.,0.,chEffErr[0][0],chEffErr[1][0]);
1712  }
1713  if (effVsSt[0]) for (Int_t iSt = 0; iSt < 6; ++iSt) {
1714  effVsSt[0]->SetPoint(iSt,iSt+1,stEff[iSt]);
1715  effVsSt[0]->SetPointError(iSt,0.,0.,stEffErr[iSt][0],stEffErr[iSt][1]);
1716  }
1717 
1718  }
1719 
1720  // print results
1721  if (print) {
1722  for (Int_t iCh = 1; iCh < 11; ++iCh) {
1723  cout << "Efficiency chamber " << iCh << " : ";
1724  cout << chEff[iCh] << " + " << chEffErr[1][iCh] << " - " << chEffErr[0][iCh] << endl;
1725  }
1726  for (Int_t iSt = 0; iSt < 6; ++iSt) {
1727  if (iSt < 5) cout << "Station " << iSt+1 << " = ";
1728  else cout << "Station 45 = ";
1729  cout << stEff[iSt] << " + " << stEffErr[iSt][1] << " - " << stEffErr[iSt][0] << endl;
1730  }
1731  cout << "Total tracking efficiency : ";
1732  cout << chEff[0] << " + " << chEffErr[1][0] << " - " << chEffErr[0][0] << endl << endl;
1733  }
1734 
1735 }
1736 
1737 
1738 //---------------------------------------------------------------------------
1740  TGraphAsymmErrors &effVsX, Int_t ip, Double_t xp)
1741 {
1744 
1745  if (!runWeights) {
1746  effVsX.SetPoint(ip,xp,-1.);
1747  effVsX.SetPointError(ip,0.,0.,0.,0.);
1748  return;
1749  }
1750 
1751  // initialize
1752  Double_t rec[2] = {0., 0.};
1753  Double_t gen = 0.;
1754  Double_t intEffErr[2] = {0., 0.};
1755  Bool_t ok = kFALSE;
1756 
1757  // loop over runs
1758  Int_t nRuns = effVsRunLow.GetN();
1759  for (Int_t iRun = 0; iRun < nRuns; ++iRun) {
1760 
1761  // get run weight
1762  TString sRun = effVsRunLow.GetXaxis()->GetBinLabel(iRun+1);
1763  TParameter<Double_t> *weight = static_cast<TParameter<Double_t>*>(runWeights->FindObject(sRun.Data()));
1764  if (!weight) {
1765  printf("weight not found for run %s\n", sRun.Data());
1766  continue;
1767  }
1768  Double_t w = weight->GetVal();
1769  Double_t w2 = w*w;
1770 
1771  // get efficiency and error
1772  Double_t eff[2] = {effVsRunLow.GetY()[iRun], effVsRunUp.GetY()[iRun]};
1773  Double_t effErr[2] = {effVsRunLow.GetErrorYlow(iRun), effVsRunUp.GetErrorYhigh(iRun)};
1774  if (eff[0] < 0. || eff[1] < 0.) {
1775  printf("no efficiency measurement --> use 0(1) ± 0 as lower(upper) limit\n");
1776  eff[0] = 0.;
1777  eff[1] = 1.;
1778  effErr[0] = 0.;
1779  effErr[1] = 0.;
1780  } else ok = kTRUE;
1781 
1782  // integrate
1783  gen += w;
1784  for (Int_t i = 0; i < 2; ++i) {
1785  rec[i] += w*eff[i];
1786  intEffErr[i] += w2*effErr[i]*effErr[i];
1787  }
1788 
1789  }
1790 
1791  // compute efficiency
1792  if (gen > 0. && ok) {
1793 
1794  effVsX.SetPoint(ip,xp,rec[1]/gen);
1795  effVsX.SetPointError(ip,0.,0.,(rec[1]-rec[0]+TMath::Sqrt(intEffErr[0]))/gen,TMath::Sqrt(intEffErr[1])/gen);
1796 
1797  } else {
1798 
1799  if (gen <= 0.) printf("impossible to integrate, all weights = 0 or unknown ?!?\n");
1800  else printf("efficiency never measured --> return -1 ± 0\n");
1801 
1802  effVsX.SetPoint(ip,xp,-1.);
1803  effVsX.SetPointError(ip,0.,0.,0.,0.);
1804 
1805  }
1806 
1807 }
1808 
1809 
1810 //---------------------------------------------------------------------------
1812 {
1816 
1817  if (runWeights) return;
1818 
1819  ifstream inFile(fileName.Data());
1820  if (!inFile.is_open()) {
1821  printf("cannot open file %s\n", fileName.Data());
1822  return;
1823  }
1824 
1825  runWeights = new THashList(1000);
1826  runWeights->SetOwner();
1827 
1828  TString line;
1829  while (! inFile.eof() ) {
1830 
1831  line.ReadLine(inFile,kTRUE);
1832  if(line.IsNull()) continue;
1833 
1834  TObjArray *param = line.Tokenize(" ");
1835  if (param->GetEntries() != 2) {
1836  printf("bad input line %s", line.Data());
1837  continue;
1838  }
1839 
1840  Int_t run = ((TObjString*)param->UncheckedAt(0))->String().Atoi();
1841  if (run < 0) {
1842  printf("invalid run number: %d", run);
1843  continue;
1844  }
1845 
1846  Float_t weight = ((TObjString*)param->UncheckedAt(1))->String().Atof();
1847  if (weight <= 0.) {
1848  printf("invalid weight: %g", weight);
1849  continue;
1850  }
1851 
1852  runWeights->Add(new TParameter<Double_t>(Form("%d",run), weight));
1853 
1854  delete param;
1855  }
1856 
1857  inFile.close();
1858 
1859 }
1860 
1861 
1862 //---------------------------------------------------------------------------
1863 void SetCentPtCh(THnSparse& SparseData)
1864 {
1867 
1868  if (ptMax >= 0 && ptMax < ptMin) {
1869  printf("Wrong Pt range, ptMax must be higher than ptMin\n");
1870  exit(-1);
1871  }
1872 
1873  if (charge < -1 || charge > 1) {
1874  printf("Selected charge must be 0, 1 or -1\n");
1875  exit(-1);
1876  }
1877 
1878  // adjust centrality range
1879  centMax = TMath::Max(centMax-1.e-12, centMin);
1880 
1881  // set the centrality range for integration
1882  Int_t lowBin = SparseData.GetAxis(1)->FindBin(centMin);
1883  Int_t upBin = SparseData.GetAxis(1)->FindBin(centMax);
1884  SparseData.GetAxis(1)->SetRange(lowBin, upBin);
1885 
1886  // set the pt range for integration
1887  lowBin = SparseData.GetAxis(2)->FindBin(ptMin);
1888  if (ptMax == -1) { upBin = SparseData.GetAxis(2)->GetNbins()+1; }
1889  else { upBin = SparseData.GetAxis(2)->FindBin(ptMax);}
1890  SparseData.GetAxis(2)->SetRange(lowBin, upBin);
1891 
1892  // set the charge range
1893  if (charge != 0) {
1894  lowBin = SparseData.GetAxis(5)->FindBin(charge);
1895  SparseData.GetAxis(5)->SetRange(lowBin, lowBin);
1896  }
1897 
1898 }
1899 
1900 
1901 //---------------------------------------------------------------------------
1902 TGraphAsymmErrors* CreateGraph(const char* name, const char* title, int value)
1903 {
1905 
1907 
1908  if ( value >= 0 ) {
1909 
1910  g->SetName(Form(name,value));
1911  g->SetTitle(Form(title,value));
1912 
1913  } else {
1914 
1915  g->SetName(name);
1916  g->SetTitle(title);
1917 
1918 
1919  }
1920 
1921  return g;
1922 
1923 }
1924 
1925 
1926 //---------------------------------------------------------------------------
1927 void BeautifyGraph(TGraphAsymmErrors &g, const char* xAxisName, const char* yAxisName)
1928 {
1930 
1931  g.SetLineStyle(1);
1932  g.SetLineColor(1);
1933  g.SetMarkerStyle(20);
1934  g.SetMarkerSize(0.7);
1935  g.SetMarkerColor(2);
1936 
1937  g.GetXaxis()->SetTitle(xAxisName);
1938  g.GetXaxis()->CenterTitle(kTRUE);
1939  g.GetXaxis()->SetLabelFont(22);
1940  g.GetXaxis()->SetTitleFont(22);
1941 
1942  g.GetYaxis()->SetTitle(yAxisName);
1943  g.GetYaxis()->CenterTitle(kTRUE);
1944  g.GetYaxis()->SetLabelFont(22);
1945  g.GetYaxis()->SetTitleFont(22);
1946 
1947 }
1948 
1949 
1950 //---------------------------------------------------------------------------
1951 void BeautifyGraphs(TObjArray& array, const char* xAxisName, const char* yAxisName)
1952 {
1954 
1955  for ( Int_t i = 0; i <= array.GetLast(); ++i)
1956  BeautifyGraph(*static_cast<TGraphAsymmErrors*>(array.UncheckedAt(i)), xAxisName, yAxisName);
1957 
1958 }
1959 
1960 
1961 //---------------------------------------------------------------------------
1962 void SetRunLabel(TGraphAsymmErrors &g, Int_t irun, const TList& runs)
1963 {
1965 
1966  g.GetXaxis()->Set(irun+1, -0.5, irun+0.5);
1967 
1968  TIter nextRun(&runs);
1969  TObjString *srun = 0x0;
1970  Int_t iirun = 1;
1971  while ((srun = static_cast<TObjString*>(nextRun()))) {
1972  g.GetXaxis()->SetBinLabel(iirun, srun->GetName());
1973  iirun++;
1974  }
1975 
1976 }
1977 
1978 
1979 //---------------------------------------------------------------------------
1980 void SetRunLabel(TObjArray& array, Int_t irun, const TList& runs)
1981 {
1983 
1984  for (Int_t i = 0; i <= array.GetLast(); ++i)
1985  SetRunLabel(*static_cast<TGraphAsymmErrors*>(array.UncheckedAt(i)), irun, runs);
1986 
1987 }
1988 
1989 
Int_t charge
void ComputeStationEfficiency(TArrayD &chEff, TArrayD chEffErr[2], Int_t iSt, Double_t &stEff, Double_t stEffErr[2])
void MuonTrackingEfficiency(TString runList="runList.txt", TString fileNameWeights="", TString objNameExtension="", TString fileNameData="AnalysisResults.root", TString fileNameSave="efficiency_new.root")
double Double_t
Definition: External.C:58
void GetStationEfficiency(TArrayD &chEff, TArrayD chEffErr[2], Int_t iSt, TGraphAsymmErrors *effVsX[3], Int_t ip, Double_t xp)
Definition: External.C:236
const char * title
Definition: MakeQAPdf.C:26
void GetTrackingEfficiency(TArrayD &chEff, TArrayD chEffErr[2], TGraphAsymmErrors *effVsSt[3], TGraphAsymmErrors *effVsX[3], Int_t ip, Double_t xp, Bool_t print=kFALSE)
void ComputeTrackingEfficiency(Double_t stEff[6], Double_t stEffErr[6][2], Double_t &spectroEff, Double_t spectroEffErr[2])
void IntegrateMuonEfficiency(TGraphAsymmErrors &effVsRunLow, TGraphAsymmErrors &effVsRunUp, TGraphAsymmErrors &effVsX, Int_t ip, Double_t xp)
TString fileName
void PlotMuonEfficiencyVsRun(TString runList, TString fileNameData, TString fileNameSave, Bool_t print, Bool_t draw)
THashList * runWeights
char Char_t
Definition: External.C:18
void PlotMuonEfficiency(TString fileNameData, TString fileNameSave, Bool_t saveEdges, Bool_t print, Bool_t draw)
TCanvas * c
Definition: TestFitELoss.C:172
void PlotIntegratedMuonEfficiency(TString fileNameWeights, TString fileNameSave, Bool_t print, Bool_t draw)
void PlotMuonEfficiencyVsXY(TString xVar, TString yVar, TString fileNameData, TString fileNameSave, Bool_t draw, Bool_t rap=kFALSE)
Double_t ptMin
Bool_t moreTrackCandidates
Int_t nDE
void LoadRunWeights(TString fileName)
void GetDEEfficiency(THnSparse &TT, THnSparse &TD, TGraphAsymmErrors &effVsDE)
TString gObjNameExtension
int Int_t
Definition: External.C:63
void SetRunLabel(TGraphAsymmErrors &g, Int_t irun, const TList &runs)
float Float_t
Definition: External.C:68
void PlotIntegratedMuonEfficiencyVsX(TString var, TString runList, TString fileNameWeights, TString fileNameData, TString fileNameSave, Bool_t print, Bool_t draw)
void BeautifyGraphs(TObjArray &array, const char *xAxisName, const char *yAxisName)
Definition: External.C:212
TGraphAsymmErrors * CreateGraph(const char *name, const char *title, int value=-1)
void PlotIntegratedMuonEfficiencyPerDE(TString fileNameWeights, TString fileNameSave)
void ComputeStation45Efficiency(TArrayD &chEff, TArrayD chEffErr[2], Double_t &st45Eff, Double_t st45EffErr[2])
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
const Char_t * effErrMode
void BeautifyGraph(TGraphAsymmErrors &g, const char *xAxisName, const char *yAxisName)
Double_t centMax
Bool_t draw[nPtBins]
void PlotMuonEfficiencyVsX(TString var, TString fileNameData, TString fileNameSave, Bool_t saveEdges, Bool_t print, Bool_t draw)
TFile * file
TList with histograms for a given trigger.
bool Bool_t
Definition: External.C:53
void PlotMuonEfficiencyPerDE(TString fileNameData, TString fileNameSave, Bool_t saveEdges)
Double_t ptMax
void PlotMuonEfficiencyPerDEVsRun(TString runList, TString fileNameData, TString fileNameSave)
void GetStation45Efficiency(TArrayD &chEff, TArrayD chEffErr[2], TGraphAsymmErrors *effVsX[3], Int_t ip, Double_t xp)
Bool_t GetChamberEfficiency(THnSparse &TT, THnSparse &TD, TArrayD &chEff, TArrayD chEffErr[2], Bool_t printError=kFALSE)
void SetCentPtCh(THnSparse &SparseData)
Double_t centMin