AliPhysics  a9863a5 (a9863a5)
 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.;
51 Double_t centMax = 999.;
52 Double_t ptMin = 0.;
53 Double_t ptMax = -1.;
54 Int_t charge = 0; // 0 selects + and -, -1 and +1 selects - or + muons respectively
55 
56 Bool_t moreTrackCandidates = kFALSE;
57 
58 THashList *runWeights = 0x0;
59 TString gObjNameExtension = "";
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 
88 void LoadRunWeights(TString fileName);
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  if (saveEdges) {
591 
592  // lower = 0 ± 0
593  effVsCh[1]->SetPoint(0,0.,0.);
594  effVsCh[1]->SetPointError(0,0.,0.,0.,0.);
595 
596  // upper = 1 ± 0
597  effVsCh[2]->SetPoint(0,0.,1.);
598  effVsCh[2]->SetPointError(0,0.,0.,0.,0.);
599 
600  }
601 
602  }
603 
604  // fill graph vs chamber
605  for (Int_t iCh = 1; iCh < 11; iCh++) {
606  for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i) {
607  effVsCh[i]->SetPoint(iCh,iCh,chEff[iCh]);
608  effVsCh[i]->SetPointError(iCh,0.,0.,chEffErr[0][iCh],chEffErr[1][iCh]);
609  }
610  }
611 
612  // close input file
613  file->Close();
614 
615  // display
616  for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i) {
617  effVsCh[i]->GetXaxis()->Set(12, -0.5, 10.5);
618  effVsCh[i]->GetXaxis()->SetNdivisions(11);
619  BeautifyGraph(*effVsCh[i], "chamber", "efficiency");
620  effVsSt[i]->GetXaxis()->Set(7, 0.5, 6.5);
621  effVsSt[i]->GetXaxis()->SetNdivisions(6);
622  BeautifyGraph(*effVsSt[i], "station", "efficiency");
623  }
624 
625  if (draw) {
626  TCanvas *c = new TCanvas("cEfficiency", "Measured tracking efficiency" , 1000, 400);
627  c->Divide(2,1);
628  gROOT->SetSelectedPad(c->cd(1));
629  effVsCh[0]->DrawClone("ap");
630  gROOT->SetSelectedPad(c->cd(2));
631  effVsSt[0]->DrawClone("ap");
632  }
633 
634  // save output
635  file = new TFile(fileNameSave.Data(),"update");
636  for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i) effVsCh[i]->Write(0x0, TObject::kOverwrite);
637  for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i) effVsSt[i]->Write(0x0, TObject::kOverwrite);
638  file->Close();
639 
640  // clean memory
641  for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i) {
642  delete effVsCh[i];
643  delete effVsSt[i];
644  }
645 
646 }
647 
648 
649 //---------------------------------------------------------------------------
650 void PlotMuonEfficiencyVsRun(TString runList, TString fileNameData, TString fileNameSave, Bool_t print, Bool_t draw)
651 {
653 
654  printf("plotting efficiency versus run...\n");
655 
656  // open run list
657  ifstream inFile(runList.Data());
658  if (!inFile.is_open()) {
659  printf("cannot open file %s\n",runList.Data());
660  return;
661  }
662 
663  // output graphs
664  TObjArray chamberVsRunGraphs; // 10 graphs: 1 for each chamber
665  chamberVsRunGraphs.SetOwner(kTRUE);
666  for ( Int_t iCh = 1; iCh < 11; ++iCh)
667  chamberVsRunGraphs.Add(CreateGraph("effCh%dVsRun", "Measured efficiency for chamber %d versus run",iCh));
668 
669  TObjArray stationVsRunGraphs[3]; // 6 graphs: 1 for each station, and 1 for station 4&5
670  TGraphAsymmErrors *trkVsRun[3];
671  TString nameAdd[3] = {"", "Low", "Up"};
672  TString titleAdd[3] = {"", " - lower limit", " - upper limit"};
673  for (Int_t i = 0; i < 3; ++i) {
674 
675  stationVsRunGraphs[i].SetOwner(kTRUE);
676  for ( Int_t iSt = 1; iSt < 6; ++iSt)
677  stationVsRunGraphs[i].Add(CreateGraph(Form("effSt%%dVsRun%s",nameAdd[i].Data()),
678  Form("Measured efficiency for station %%d versus run%s",titleAdd[i].Data()),iSt));
679  stationVsRunGraphs[i].Add(CreateGraph(Form("effSt4&5VsRun%s",nameAdd[i].Data()),
680  Form("Measured efficiency for station 4&5 versus run%s",titleAdd[i].Data())));
681 
682  trkVsRun[i] = CreateGraph(Form("trackingEffVsRun%s",nameAdd[i].Data()),
683  Form("Measured tracking efficiency versus run%s", titleAdd[i].Data()));
684 
685  }
686 
687  Int_t irun = -1;
688  TList runs;
689  runs.SetOwner();
690 
691  // loop over runs
692  while (!inFile.eof()) {
693 
694  // get current run number
695  TString currRun;
696  currRun.ReadLine(inFile,kTRUE);
697  if(currRun.IsNull()) continue;
698  runs.AddLast(new TObjString(currRun));
699  irun++;
700 
701  Int_t run = currRun.Atoi();
702 
703  printf("run %d: ", run);
704 
705  // compute efficiencies for this run
706  TString dataFile = Form("runs/%d/%s", run, fileNameData.Data());
707  TString outFile = Form("runs/%d/%s", run, fileNameSave.Data());
708  PlotMuonEfficiency(dataFile, outFile, kTRUE, print, kFALSE);
709 
710  TFile *file = new TFile(outFile.Data(), "read");
711  if (!file || !file->IsOpen()) {
712  printf("cannot open file\n");
713  continue;
714  }
715 
716  // loop over central value and edges
717  for (Int_t i = 0; i < 3; ++i) {
718 
719  // get results
720  TGraphAsymmErrors *effVsCh;
721  TGraphAsymmErrors *effVsSt;
722  effVsCh = static_cast<TGraphAsymmErrors*>(file->FindObjectAny(Form("chamberEff%s",nameAdd[i].Data())));
723  effVsSt = static_cast<TGraphAsymmErrors*>(file->FindObjectAny(Form("stationEff%s",nameAdd[i].Data())));
724  if (!effVsCh || !effVsSt) {
725  printf("efficiency graph not found\n");
726  continue;
727  }
728 
729  // fill chamber efficiency
730  if (i == 0) for ( Int_t iCh = 0; iCh < 10; ++iCh) {
731  TGraphAsymmErrors *g = static_cast<TGraphAsymmErrors*>(chamberVsRunGraphs.UncheckedAt(iCh));
732  g->SetPoint(irun,irun,effVsCh->GetY()[iCh+1]);
733  g->SetPointError(irun,0.,0.,effVsCh->GetErrorYlow(iCh+1),effVsCh->GetErrorYhigh(iCh+1));
734  }
735 
736  // fill station efficiency
737  for ( Int_t iSt = 0; iSt < 6; ++iSt) {
738  TGraphAsymmErrors *g = static_cast<TGraphAsymmErrors*>(stationVsRunGraphs[i].UncheckedAt(iSt));
739  g->SetPoint(irun,irun,effVsSt->GetY()[iSt]);
740  g->SetPointError(irun,0.,0.,effVsSt->GetErrorYlow(iSt),effVsSt->GetErrorYhigh(iSt));
741  }
742 
743  // fill spectrometer efficiency
744  trkVsRun[i]->SetPoint(irun,irun,effVsCh->GetY()[0]);
745  trkVsRun[i]->SetPointError(irun,0.,0.,effVsCh->GetErrorYlow(0),effVsCh->GetErrorYhigh(0));
746 
747  }
748 
749  file->Close();
750 
751  }
752 
753  inFile.close();
754 
755  // display
756  BeautifyGraphs(chamberVsRunGraphs,"run number","efficiency");
757  SetRunLabel(chamberVsRunGraphs,irun,runs);
758  for (Int_t i = 0; i < 3; ++i) {
759  BeautifyGraph(*trkVsRun[i],"run number","efficiency");
760  SetRunLabel(*trkVsRun[i],irun,runs);
761  BeautifyGraphs(stationVsRunGraphs[i],"run number","efficiency");
762  SetRunLabel(stationVsRunGraphs[i],irun,runs);
763  }
764 
765  if (draw) {
766  new TCanvas("cTrackingEffVsRun", "Tracking efficiency versus run",1000,400);
767  trkVsRun[0]->DrawClone("ap");
768  }
769 
770  // save output
771  TFile* file = new TFile(fileNameSave.Data(),"update");
772  chamberVsRunGraphs.Write("ChamberEffVsRun", TObject::kOverwrite | TObject::kSingleKey);
773  for (Int_t i = 0; i < 3; ++i)
774  stationVsRunGraphs[i].Write(Form("StationEffVsRun%s",nameAdd[i].Data()), TObject::kOverwrite | TObject::kSingleKey);
775  for (Int_t i = 0; i < 3; ++i) trkVsRun[i]->Write(0x0, TObject::kOverwrite);
776  file->Close();
777 
778  // clean memory
779  for (Int_t i = 0; i < 3; ++i) delete trkVsRun[i];
780 
781 }
782 
783 
784 //---------------------------------------------------------------------------
785 void PlotIntegratedMuonEfficiency(TString fileNameWeights, TString fileNameSave, Bool_t print, Bool_t draw)
786 {
788 
789  printf("plotting integrated efficiency...\n");
790 
791  // load run weights
792  if (!fileNameWeights.IsNull()) LoadRunWeights(fileNameWeights);
793  if (!runWeights) {
794  printf("Cannot compute integrated efficiency without run-by-run weights\n");
795  return;
796  }
797 
798  // get input hists
799  TFile *file = new TFile(fileNameSave.Data(), "update");
800  if (!file || !file->IsOpen()) {
801  printf("cannot open file\n");
802  return;
803  }
804  TObjArray *chamberVsRunGraphs = static_cast<TObjArray*>(file->FindObjectAny("ChamberEffVsRun"));
805  TObjArray *stationVsRunGraphs[2];
806  stationVsRunGraphs[0] = static_cast<TObjArray*>(file->FindObjectAny("StationEffVsRunLow"));
807  stationVsRunGraphs[1] = static_cast<TObjArray*>(file->FindObjectAny("StationEffVsRunUp"));
808  TGraphAsymmErrors *trkVsRun[2];
809  trkVsRun[0] = static_cast<TGraphAsymmErrors*>(file->FindObjectAny("trackingEffVsRunLow"));
810  trkVsRun[1] = static_cast<TGraphAsymmErrors*>(file->FindObjectAny("trackingEffVsRunUp"));
811  if (!chamberVsRunGraphs || !stationVsRunGraphs[0] || !stationVsRunGraphs[1] || !trkVsRun[0] || !trkVsRun[1]) {
812  printf("object not found --> you must first plot the efficiency versus run\n");
813  return;
814  }
815 
816  // output graphs
817  TGraphAsymmErrors *effVsCh = CreateGraph("integratedChamberEff", "Integrated efficiency per chamber (0 = spectro)");
818  TGraphAsymmErrors *effVsSt = CreateGraph("integratedStationEff", "Integrated efficiency per station (6 = st4&5)");
819 
820  // integrate spectrometer efficiency
821  IntegrateMuonEfficiency(*trkVsRun[0], *trkVsRun[1], *effVsCh, 0, 0.);
822 
823  // integrate chamber efficiency
824  for ( Int_t iCh = 0; iCh < 10; ++iCh) {
825  TGraphAsymmErrors *g = static_cast<TGraphAsymmErrors*>(chamberVsRunGraphs->UncheckedAt(iCh));
826  IntegrateMuonEfficiency(*g, *g, *effVsCh, iCh+1, iCh+1);
827  }
828 
829  // integrate station efficiency
830  for ( Int_t iSt = 0; iSt < 6; ++iSt) {
831  TGraphAsymmErrors *gLow = static_cast<TGraphAsymmErrors*>(stationVsRunGraphs[0]->UncheckedAt(iSt));
832  TGraphAsymmErrors *gUp = static_cast<TGraphAsymmErrors*>(stationVsRunGraphs[1]->UncheckedAt(iSt));
833  IntegrateMuonEfficiency(*gLow, *gUp, *effVsSt, iSt, iSt+1);
834  }
835 
836  // print results
837  if (print) {
838  for (Int_t iCh = 1; iCh < 11; ++iCh) {
839  cout << "Efficiency chamber " << iCh << " : ";
840  cout << effVsCh->GetY()[iCh] << " + " << effVsCh->GetErrorYhigh(iCh) << " - " << effVsCh->GetErrorYlow(iCh) << endl;
841  }
842  for (Int_t iSt = 0; iSt < 6; ++iSt) {
843  if (iSt < 5) cout << "Station " << iSt+1 << " = ";
844  else cout << "Station 45 = ";
845  cout << effVsSt->GetY()[iSt] << " + " << effVsSt->GetErrorYhigh(iSt) << " - " << effVsSt->GetErrorYlow(iSt) << endl;
846  }
847  cout << "Total tracking efficiency : ";
848  cout << effVsCh->GetY()[0] << " + " << effVsCh->GetErrorYhigh(0) << " - " << effVsCh->GetErrorYlow(0) << endl << endl;
849  }
850 
851  // display
852  effVsCh->GetXaxis()->Set(12, -0.5, 10.5);
853  effVsCh->GetXaxis()->SetNdivisions(11);
854  BeautifyGraph(*effVsCh, "chamber", "efficiency");
855  effVsSt->GetXaxis()->Set(7, 0.5, 6.5);
856  effVsSt->GetXaxis()->SetNdivisions(6);
857  BeautifyGraph(*effVsSt, "station", "efficiency");
858 
859  if (draw) {
860  TCanvas *c = new TCanvas("cIntegratedEfficiency", "Integrated tracking efficiency" , 1000, 400);
861  c->Divide(2,1);
862  gROOT->SetSelectedPad(c->cd(1));
863  effVsCh->DrawClone("ap");
864  gROOT->SetSelectedPad(c->cd(2));
865  effVsSt->DrawClone("ap");
866  }
867 
868  // save output
869  effVsCh->Write(0x0, TObject::kOverwrite);
870  effVsSt->Write(0x0, TObject::kOverwrite);
871  file->Close();
872 
873  // clean memory
874  delete effVsCh;
875  delete effVsSt;
876 
877 }
878 
879 
880 //---------------------------------------------------------------------------
881 void PlotMuonEfficiencyPerDE(TString fileNameData, TString fileNameSave, Bool_t saveEdges)
882 {
884 
885  printf("plotting efficiency per DE...\n");
886 
887  // get input hists
888  TFile *file = new TFile(fileNameData.Data(), "read");
889  if (!file || !file->IsOpen()) {
890  printf("cannot open file %s \n",fileNameData.Data());
891  return;
892  }
893  TList *listTT = static_cast<TList*>(file->FindObjectAny(Form("TotalTracksPerChamber%s", gObjNameExtension.Data())));
894  TList *listTD = static_cast<TList*>(file->FindObjectAny(Form("TracksDetectedPerChamber%s", gObjNameExtension.Data())));
895 
896  // output graph arrays
897  TObjArray chamberVsDEGraphs; // 10 graphs: 1 for each chamber
898  chamberVsDEGraphs.SetOwner(kTRUE);
899  TGraphAsymmErrors *gCh[10]; // shortcut
900 
901  TObjArray stationVsDEGraphs[3]; // 6 graphs: 1 for each station, and 1 for station 4&5
902  TString nameAdd[3] = {"", "Low", "Up"};
903  TString titleAdd[3] = {"", " - lower limit", " - upper limit"};
904  for (Int_t i = 0; i < 3; ++i) stationVsDEGraphs[i].SetOwner(kTRUE);
905  TGraphAsymmErrors *gSt[6][3]; // shortcut
906  for (Int_t iSt = 0; iSt < 6; ++iSt) for (Int_t i = 0; i < 3; ++i) gSt[iSt][i] = 0x0;
907 
908  // loop over chambers
909  for (Int_t iCh = 0; iCh < 10; ++iCh) {
910 
911  // output graphs
912  chamberVsDEGraphs.Add(CreateGraph("effCh%dVsDE","Measured efficiency for chamber %d per DE",iCh+1));
913  gCh[iCh] = static_cast<TGraphAsymmErrors*>(chamberVsDEGraphs.UncheckedAt(iCh));
914 
915  // get input hists
916  THnSparse *TT = static_cast<THnSparse*>(listTT->At(iCh));
917  THnSparse *TD = static_cast<THnSparse*>(listTD->At(iCh));
918 
919  // set the centrality and pT range for integration
920  SetCentPtCh(*TT);
921  SetCentPtCh(*TD);
922 
923  // compute DE efficency and errors
924  GetDEEfficiency(*TT, *TD, *gCh[iCh]);
925 
926  }
927 
928  // close input file
929  file->Close();
930 
931  TArrayD chEff(11);
932  TArrayD chEffErr[2];
933  chEffErr[0].Set(11);
934  chEffErr[1].Set(11);
935 
936  // loop over the first 3 stations
937  for (Int_t iSt = 0; iSt < 3; ++iSt) {
938 
939  // output graphs
940  for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i) {
941  stationVsDEGraphs[i].Add(CreateGraph(Form("effSt%%dVsDE%s",nameAdd[i].Data()),
942  Form("Measured efficiency for station %%d per DE%s",titleAdd[i].Data()),iSt+1));
943  gSt[iSt][i] = static_cast<TGraphAsymmErrors*>(stationVsDEGraphs[i].UncheckedAt(iSt));
944  }
945 
946  // loop over DE
947  Int_t nDE = gCh[2*iSt]->GetN();
948  for (Int_t iDE = 0; iDE < nDE; ++iDE) {
949 
950  // copy efficiency
951  for (Int_t iCh = 0; iCh < 2; ++iCh) {
952  chEff[2*iSt+iCh+1] = gCh[2*iSt+iCh]->GetY()[iDE];
953  chEffErr[0][2*iSt+iCh+1] = gCh[2*iSt+iCh]->GetErrorYlow(iDE);
954  chEffErr[1][2*iSt+iCh+1] = gCh[2*iSt+iCh]->GetErrorYhigh(iDE);
955  }
956 
957  // compute station efficiency
958  GetStationEfficiency(chEff, chEffErr, iSt, gSt[iSt], iDE, iDE);
959 
960  }
961 
962  }
963 
964  // output graphs for last 2 stations
965  for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i) {
966  for (Int_t iSt = 3; iSt < 5; ++iSt) {
967  stationVsDEGraphs[i].Add(CreateGraph(Form("effSt%%dVsDE%s",nameAdd[i].Data()),
968  Form("Measured efficiency for station %%d per DE%s",titleAdd[i].Data()),iSt+1));
969  gSt[iSt][i] = static_cast<TGraphAsymmErrors*>(stationVsDEGraphs[i].UncheckedAt(iSt));
970  }
971  stationVsDEGraphs[i].Add(CreateGraph(Form("effSt4&5VsDE%s",nameAdd[i].Data()),
972  Form("Measured efficiency for station 4&5 per DE%s",titleAdd[i].Data())));
973  gSt[5][i] = static_cast<TGraphAsymmErrors*>(stationVsDEGraphs[i].UncheckedAt(5));
974  }
975 
976  // loop over DE
977  Int_t nDE = gCh[6]->GetN();
978  for (Int_t iDE = 0; iDE < nDE; ++iDE) {
979 
980  // copy efficiency
981  for (Int_t iCh = 6; iCh < 10; ++iCh) {
982  chEff[iCh+1] = gCh[iCh]->GetY()[iDE];
983  chEffErr[0][iCh+1] = gCh[iCh]->GetErrorYlow(iDE);
984  chEffErr[1][iCh+1] = gCh[iCh]->GetErrorYhigh(iDE);
985  }
986 
987  // compute station 4&5 efficiency individually
988  for (Int_t iSt = 3; iSt < 5; ++iSt) GetStationEfficiency(chEff, chEffErr, iSt, gSt[iSt], iDE, iDE);
989 
990  // compute station 4&5 efficiency together
991  GetStation45Efficiency(chEff, chEffErr, gSt[5], iDE, iDE);
992 
993  }
994 
995  // display
996  for (Int_t iCh = 0; iCh < 10; ++iCh) {
997  Int_t nDE = gCh[iCh]->GetN();
998  gCh[iCh]->GetXaxis()->Set(nDE+1, -0.5, nDE-0.5);
999  gCh[iCh]->GetXaxis()->SetNdivisions(nDE);
1000  }
1001  BeautifyGraphs(chamberVsDEGraphs,"Detection Element","efficiency");
1002  for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i) {
1003  for (Int_t iSt = 0; iSt < 6; ++iSt) {
1004  Int_t nDE = gSt[iSt][i]->GetN();
1005  gSt[iSt][i]->GetXaxis()->Set(nDE+1, -0.5, nDE-0.5);
1006  gSt[iSt][i]->GetXaxis()->SetNdivisions(nDE);
1007  }
1008  BeautifyGraphs(stationVsDEGraphs[i],"Detection Element","efficiency");
1009 
1010  }
1011 
1012  // save Output
1013  file = new TFile(fileNameSave.Data(),"update");
1014  chamberVsDEGraphs.Write("ChamberEffVsDE", TObject::kOverwrite | TObject::kSingleKey);
1015  for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i)
1016  stationVsDEGraphs[i].Write(Form("StationEffVsDE%s",nameAdd[i].Data()), TObject::kOverwrite | TObject::kSingleKey);
1017  file->Close();
1018 
1019 }
1020 
1021 
1022 //---------------------------------------------------------------------------
1023 void PlotMuonEfficiencyPerDEVsRun(TString runList, TString fileNameData, TString fileNameSave)
1024 {
1026 
1027  printf("plotting efficiency per DE versus run...\n");
1028 
1029  // open run list
1030  ifstream inFile(runList.Data());
1031  if (!inFile.is_open()) {
1032  printf("cannot open file %s\n",runList.Data());
1033  return;
1034  }
1035 
1036  // output graphs
1037  TObjArray deVsRunGraphs; // 1 graph per DE
1038  deVsRunGraphs.SetOwner(kTRUE);
1039  TObjArray stDEVsRunGraphs[3]; // 1 graph per pair (quartet) of DE in individual stations (stations 4&5)
1040  TString nameAdd[3] = {"", "Low", "Up"};
1041  TString titleAdd[3] = {"", " - lower limit", " - upper limit"};
1042  for (Int_t i = 0; i < 3; ++i) stDEVsRunGraphs[i].SetOwner(kTRUE);
1043  Bool_t createGraph = kTRUE;
1044 
1045  Int_t irun = -1;
1046  TList runs;
1047  runs.SetOwner();
1048 
1049  // loop over runs
1050  while (!inFile.eof()) {
1051 
1052  // get current run number
1053  TString currRun;
1054  currRun.ReadLine(inFile,kTRUE);
1055  if(currRun.IsNull()) continue;
1056  runs.AddLast(new TObjString(currRun));
1057  irun++;
1058 
1059  Int_t run = currRun.Atoi();
1060 
1061  printf("run %d: ", run);
1062 
1063  // compute efficiencies for this run
1064  TString dataFile = Form("runs/%d/%s", run, fileNameData.Data());
1065  TString outFile = Form("runs/%d/%s", run, fileNameSave.Data());
1066  PlotMuonEfficiencyPerDE(dataFile, outFile, kTRUE);
1067 
1068  TFile *file = new TFile(outFile.Data(), "read");
1069  if (!file || !file->IsOpen()) {
1070  printf("cannot open file\n");
1071  continue;
1072  }
1073 
1074  // get results
1075  TObjArray *chamberVsDEGraphs = static_cast<TObjArray*>(file->FindObjectAny("ChamberEffVsDE"));
1076  if (!chamberVsDEGraphs) {
1077  printf("efficiency graph not found\n");
1078  continue;
1079  }
1080 
1081  Int_t currentDE = 0;
1082 
1083  // loop over chambers
1084  for ( Int_t iCh = 0; iCh < 10; ++iCh) {
1085 
1086  // input graph
1087  TGraphAsymmErrors *g = static_cast<TGraphAsymmErrors*>(chamberVsDEGraphs->UncheckedAt(iCh));
1088  Int_t nDE = g->GetN();
1089 
1090  // loop over DE
1091  for (Int_t iDE = 0; iDE < nDE; ++iDE) {
1092 
1093  // output graph
1094  if (createGraph) deVsRunGraphs.Add(CreateGraph("effDE%dVsRun","Measured efficiency for DE %d versus run",100*(iCh+1)+iDE));
1095  TGraphAsymmErrors *gDE= static_cast<TGraphAsymmErrors*>(deVsRunGraphs.UncheckedAt(currentDE++));
1096 
1097  // fill DE efficiency
1098  gDE->SetPoint(irun,irun,g->GetY()[iDE]);
1099  gDE->SetPointError(irun,0.,0.,g->GetErrorYlow(iDE),g->GetErrorYhigh(iDE));
1100 
1101  }
1102 
1103  }
1104 
1105  // loop over central value and edges
1106  for (Int_t i = 0; i < 3; ++i) {
1107 
1108  // get results
1109  TObjArray *stationVsDEGraphs = static_cast<TObjArray*>(file->FindObjectAny(Form("StationEffVsDE%s",nameAdd[i].Data())));
1110  if (!stationVsDEGraphs) {
1111  printf("efficiency graph not found\n");
1112  continue;
1113  }
1114 
1115  Int_t currentStDE = 0;
1116 
1117  // loop over stations
1118  for ( Int_t iSt = 0; iSt < 6; ++iSt) {
1119 
1120  // input graph
1121  TGraphAsymmErrors *g = static_cast<TGraphAsymmErrors*>(stationVsDEGraphs->UncheckedAt(iSt));
1122  Int_t nDE = g->GetN();
1123 
1124  // loop over DE
1125  for (Int_t iDE = 0; iDE < nDE; ++iDE) {
1126 
1127  // output graph
1128  if (createGraph) {
1129  TString sSt = (iSt<5) ? Form("%d",iSt+1) : "4&5";
1130  stDEVsRunGraphs[i].Add(CreateGraph(Form("effSt%sDE%%dVsRun%s",sSt.Data(),nameAdd[i].Data()),
1131  Form("Measured efficiency for DE %%d in station %s versus run%s",sSt.Data(),titleAdd[i].Data()),iDE));
1132  }
1133  TGraphAsymmErrors *gDE= static_cast<TGraphAsymmErrors*>(stDEVsRunGraphs[i].UncheckedAt(currentStDE++));
1134 
1135  // fill DE efficiency
1136  gDE->SetPoint(irun,irun,g->GetY()[iDE]);
1137  gDE->SetPointError(irun,0.,0.,g->GetErrorYlow(iDE),g->GetErrorYhigh(iDE));
1138 
1139  }
1140 
1141  }
1142 
1143  }
1144 
1145  file->Close();
1146 
1147  createGraph = kFALSE;
1148 
1149  }
1150 
1151  inFile.close();
1152 
1153  // display
1154  BeautifyGraphs(deVsRunGraphs,"run number","efficiency");
1155  SetRunLabel(deVsRunGraphs,irun,runs);
1156  for (Int_t i = 0; i < 3; ++i) {
1157  BeautifyGraphs(stDEVsRunGraphs[i],"run number","efficiency");
1158  SetRunLabel(stDEVsRunGraphs[i],irun,runs);
1159  }
1160 
1161  // save output
1162  TFile* file = new TFile(fileNameSave.Data(),"update");
1163  deVsRunGraphs.Write("DEEffVsRun", TObject::kOverwrite | TObject::kSingleKey);
1164  for (Int_t i = 0; i < 3; ++i)
1165  stDEVsRunGraphs[i].Write(Form("DEEffPerStationVsRun%s",nameAdd[i].Data()), TObject::kOverwrite | TObject::kSingleKey);
1166  file->Close();
1167 
1168 }
1169 
1170 
1171 //---------------------------------------------------------------------------
1172 void PlotIntegratedMuonEfficiencyPerDE(TString fileNameWeights, TString fileNameSave)
1173 {
1175 
1176  printf("plotting integrated efficiency per DE...\n");
1177 
1178  // load run weights
1179  if (!fileNameWeights.IsNull()) LoadRunWeights(fileNameWeights);
1180  if (!runWeights) {
1181  printf("Cannot compute integrated efficiency without run-by-run weights\n");
1182  return;
1183  }
1184 
1185  // get input hists
1186  TFile *file = new TFile(fileNameSave.Data(), "update");
1187  if (!file || !file->IsOpen()) {
1188  printf("cannot open file\n");
1189  return;
1190  }
1191  TObjArray *deVsRunGraphs = static_cast<TObjArray*>(file->FindObjectAny("DEEffVsRun"));
1192  TObjArray *stDEVsRunGraphs[2];
1193  stDEVsRunGraphs[0] = static_cast<TObjArray*>(file->FindObjectAny("DEEffPerStationVsRunLow"));
1194  stDEVsRunGraphs[1] = static_cast<TObjArray*>(file->FindObjectAny("DEEffPerStationVsRunUp"));
1195  if (!deVsRunGraphs || !stDEVsRunGraphs[0] || !stDEVsRunGraphs[1]) {
1196  printf("object not found --> you must first plot the efficiency versus run\n");
1197  return;
1198  }
1199 
1200  // output graph arrays
1201  TObjArray chamberVsDEGraphs; // 10 graphs: 1 for each chamber
1202  chamberVsDEGraphs.SetOwner(kTRUE);
1203  for ( Int_t iCh = 0; iCh < 10; ++iCh)
1204  chamberVsDEGraphs.Add(CreateGraph("integratedEffCh%dVsDE","Integrated efficiency for chamber %d per DE",iCh+1));
1205 
1206  TObjArray stationVsDEGraphs; // 6 graphs: 1 for each station, and 1 for station 4&5
1207  stationVsDEGraphs.SetOwner(kTRUE);
1208  for ( Int_t iSt = 0; iSt < 5; ++iSt)
1209  stationVsDEGraphs.Add(CreateGraph("integratedEffSt%dVsDE","Integrated efficiency for station %d per DE",iSt+1));
1210  stationVsDEGraphs.Add(CreateGraph("integratedEffSt4&5VsDE","Integrated efficiency for station 4&5 per DE"));
1211 
1212  // Loop over DE
1213  TIter nextDE(deVsRunGraphs);
1214  TGraphAsymmErrors *gDE = 0x0;
1215  while ((gDE = static_cast<TGraphAsymmErrors*>(nextDE()))) {
1216 
1217  // get chamber and DE indices
1218  Int_t deId;
1219  sscanf(gDE->GetName(), "effDE%dVsRun", &deId);
1220  Int_t iCh = deId/100-1;
1221  Int_t iDE = deId%100;
1222 
1223  // integrate DE efficiency
1224  TGraphAsymmErrors *g = static_cast<TGraphAsymmErrors*>(chamberVsDEGraphs.UncheckedAt(iCh));
1225  IntegrateMuonEfficiency(*gDE, *gDE, *g, iDE, iDE);
1226 
1227  }
1228 
1229  // Loop over DE per station
1230  Int_t ng = stDEVsRunGraphs[0]->GetEntries();
1231  for (Int_t ig = 0; ig < ng; ++ig) {
1232 
1233  TGraphAsymmErrors *gDELow = static_cast<TGraphAsymmErrors*>(stDEVsRunGraphs[0]->UncheckedAt(ig));
1234  TGraphAsymmErrors *gDEUp = static_cast<TGraphAsymmErrors*>(stDEVsRunGraphs[1]->UncheckedAt(ig));
1235 
1236  // get station and DE indices
1237  Int_t iSt, iDE;
1238  if (strstr(gDELow->GetName(), "4&5")) {
1239  iSt = 5;
1240  sscanf(gDELow->GetName(), "effSt4&5DE%dVsRunLow", &iDE);
1241  } else {
1242  sscanf(gDELow->GetName(), "effSt%dDE%dVsRunLow", &iSt, &iDE);
1243  iSt--;
1244  }
1245 
1246  // Integrate DE efficiency per station
1247  TGraphAsymmErrors *g = static_cast<TGraphAsymmErrors*>(stationVsDEGraphs.UncheckedAt(iSt));
1248  IntegrateMuonEfficiency(*gDELow, *gDEUp, *g, iDE, iDE);
1249 
1250  }
1251 
1252  // display
1253  for ( Int_t iCh = 0; iCh < 10; ++iCh) {
1254  TGraphAsymmErrors *g = static_cast<TGraphAsymmErrors*>(chamberVsDEGraphs.UncheckedAt(iCh));
1255  Int_t nDE = g->GetN();
1256  g->GetXaxis()->Set(nDE+1, -0.5, nDE-0.5);
1257  g->GetXaxis()->SetNdivisions(nDE);
1258  }
1259  BeautifyGraphs(chamberVsDEGraphs,"Detection Element","efficiency");
1260  for ( Int_t iSt = 0; iSt < 6; ++iSt) {
1261  TGraphAsymmErrors *g = static_cast<TGraphAsymmErrors*>(stationVsDEGraphs.UncheckedAt(iSt));
1262  Int_t nDE = g->GetN();
1263  g->GetXaxis()->Set(nDE+1, -0.5, nDE-0.5);
1264  g->GetXaxis()->SetNdivisions(nDE);
1265  }
1266  BeautifyGraphs(stationVsDEGraphs,"Detection Element","efficiency");
1267 
1268  // save Output
1269  chamberVsDEGraphs.Write("IntegratedChamberEffVsDE", TObject::kOverwrite | TObject::kSingleKey);
1270  stationVsDEGraphs.Write("IntegratedStationEffVsDE", TObject::kOverwrite | TObject::kSingleKey);
1271  file->Close();
1272 
1273 }
1274 
1275 
1276 //---------------------------------------------------------------------------
1277 Bool_t GetChamberEfficiency(THnSparse &TT, THnSparse &TD, TArrayD &chEff, TArrayD chEffErr[2], Bool_t printError)
1278 {
1281 
1282  // project track hists over the chamber axis
1283  TH1D *TTdraw = TT.Projection(0,"e");
1284  TH1D *TDdraw = TD.Projection(0,"e");
1285 
1286  // compute chamber efficiency and errors
1287  TGraphAsymmErrors *efficiency = (TTdraw->GetEntries() > 0) ? new TGraphAsymmErrors(TDdraw, TTdraw,Form("%se0",effErrMode)) : 0x0;
1288  Bool_t ok = (efficiency);
1289 
1290  // fill arrays
1291  if (ok) {
1292 
1293  Bool_t missingEff = kFALSE;
1294 
1295  for (Int_t i = 0; i < 10; i++) {
1296 
1297  if (TTdraw->GetBinContent(i+1) > 0) {
1298 
1299  chEff[i+1] = efficiency->GetY()[i];
1300  chEffErr[0][i+1] = efficiency->GetErrorYlow(i);
1301  chEffErr[1][i+1] = efficiency->GetErrorYhigh(i);
1302 
1303  } else {
1304 
1305  chEff[i+1] = -1.;
1306  chEffErr[0][i+1] = 0.;
1307  chEffErr[1][i+1] = 0.;
1308 
1309  missingEff = kTRUE;
1310 
1311  }
1312 
1313  }
1314 
1315  if (missingEff && printError) cout << "efficiency partially unknown" << endl;
1316 
1317  } else {
1318 
1319  for (Int_t i = 0; i < 10; i++) {
1320 
1321  chEff[i+1] = -1.;
1322  chEffErr[0][i+1] = 0.;
1323  chEffErr[1][i+1] = 0.;
1324 
1325  }
1326 
1327  if (printError) cout << "efficiency unknown" << endl << endl;
1328 
1329  }
1330 
1331  // clean memory
1332  delete TTdraw;
1333  delete TDdraw;
1334  delete efficiency;
1335 
1336  return ok;
1337 
1338 }
1339 
1340 
1341 //---------------------------------------------------------------------------
1342 void GetDEEfficiency(THnSparse &TT, THnSparse &TD, TGraphAsymmErrors &effVsDE)
1343 {
1345 
1346  // project track hists over the chamber axis
1347  TH1D *TTdraw = TT.Projection(0,"e");
1348  TH1D *TDdraw = TD.Projection(0,"e");
1349 
1350  // compute DE efficiency and errors
1351  TGraphAsymmErrors *efficiency = (TTdraw->GetEntries() > 0) ? new TGraphAsymmErrors(TDdraw, TTdraw,Form("%se0",effErrMode)) : 0x0;
1352  Int_t nDE = TTdraw->GetNbinsX();
1353 
1354  if (efficiency) {
1355 
1356  for (Int_t iDE = 0; iDE < nDE; ++iDE) {
1357 
1358  if (TTdraw->GetBinContent(iDE+1) > 0) {
1359 
1360  effVsDE.SetPoint(iDE,iDE,efficiency->GetY()[iDE]);
1361  effVsDE.SetPointError(iDE,0,0,efficiency->GetErrorYlow(iDE),efficiency->GetErrorYhigh(iDE));
1362 
1363  } else {
1364 
1365  effVsDE.SetPoint(iDE,iDE,-1);
1366  effVsDE.SetPointError(iDE,0,0,0,0);
1367 
1368  }
1369 
1370  }
1371 
1372  } else {
1373 
1374  for (Int_t iDE = 0; iDE < nDE; ++iDE) {
1375 
1376  effVsDE.SetPoint(iDE,iDE,-1);
1377  effVsDE.SetPointError(iDE,0,0,0,0);
1378 
1379  }
1380 
1381  }
1382 
1383  // clean memory
1384  delete TTdraw;
1385  delete TDdraw;
1386  delete efficiency;
1387 
1388 }
1389 
1390 
1391 //---------------------------------------------------------------------------
1392 void ComputeStationEfficiency(TArrayD &chEff, TArrayD chEffErr[2], Int_t iSt, Double_t &stEff, Double_t stEffErr[2])
1393 {
1395 
1396  stEff = 1.-(1.-chEff[2*iSt+1])*(1.-chEff[2*iSt+2]);
1397 
1398  Double_t d1 = (1. - chEff[2*iSt+2]); d1 *= d1;
1399  Double_t d2 = (1. - chEff[2*iSt+1]); d2 *= d2;
1400 
1401  for (Int_t i = 0; i < 2; ++i) {
1402  Double_t s1 = chEffErr[i][2*iSt+1] * chEffErr[i][2*iSt+1];
1403  Double_t s2 = chEffErr[i][2*iSt+2] * chEffErr[i][2*iSt+2];
1404  stEffErr[i] = TMath::Sqrt(d1*s1 + d2*s2 + s1*s2);
1405  }
1406 
1407  stEffErr[0] = TMath::Min(stEff, stEffErr[0]);
1408  stEffErr[1] = TMath::Min(1.-stEff, stEffErr[1]);
1409 
1410 }
1411 
1412 
1413 //---------------------------------------------------------------------------
1414 void GetStationEfficiency(TArrayD &chEff, TArrayD chEffErr[2], Int_t iSt, TGraphAsymmErrors *effVsX[3], Int_t ip, Double_t xp)
1415 {
1417 
1418  if (chEff[2*iSt+1] >= 0 && chEff[2*iSt+2] >= 0) {
1419 
1420  // compute station efficiency from known chamber efficiency
1421  Double_t stEff, stEffErr[2];
1422  ComputeStationEfficiency(chEff, chEffErr, iSt, stEff, stEffErr);
1423 
1424  // fill graphs if required
1425  for (Int_t i = 0; i < 3; ++i) {
1426  if (effVsX[i]) {
1427  effVsX[i]->SetPoint(ip,xp,stEff);
1428  effVsX[i]->SetPointError(ip,0.,0.,stEffErr[0],stEffErr[1]);
1429  }
1430  }
1431 
1432  } else {
1433 
1434  Double_t edge[2] = {0., 1.};
1435  TArrayD chEffEdge[2];
1436  Double_t stEffEdge[2];
1437  Double_t stEffEdgeErr[2][2];
1438 
1439  for (Int_t i = 0; i < 2; ++i) {
1440 
1441  // set lower(upper) limit of chamber efficiency
1442  chEffEdge[i].Set(11);
1443  for (Int_t iCh = 1; iCh < 3; ++iCh) chEffEdge[i][2*iSt+iCh] = (chEff[2*iSt+iCh] < 0) ? edge[i] : chEff[2*iSt+iCh];
1444 
1445  // compute station efficiency
1446  ComputeStationEfficiency(chEffEdge[i], chEffErr, iSt, stEffEdge[i], stEffEdgeErr[i]);
1447 
1448  // fill graph if required
1449  if (effVsX[i+1]) {
1450  effVsX[i+1]->SetPoint(ip,xp,stEffEdge[i]);
1451  effVsX[i+1]->SetPointError(ip,0.,0.,stEffEdgeErr[i][0],stEffEdgeErr[i][1]);
1452  }
1453 
1454  }
1455 
1456  // combine extreme cases to get station efficiency and fill graph if required
1457  if (effVsX[0]) {
1458  effVsX[0]->SetPoint(ip,xp,stEffEdge[1]);
1459  effVsX[0]->SetPointError(ip,0.,0.,stEffEdge[1]-stEffEdge[0]+stEffEdgeErr[0][0],stEffEdgeErr[1][1]);
1460  }
1461 
1462  }
1463 
1464 }
1465 
1466 
1467 //---------------------------------------------------------------------------
1468 void ComputeStation45Efficiency(TArrayD &chEff, TArrayD chEffErr[2], Double_t &st45Eff, Double_t st45EffErr[2])
1469 {
1471 
1472  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];
1473 
1474  Double_t d1 = chEff[8]*chEff[9] + chEff[8]*chEff[10] + chEff[9]*chEff[10] - 3.*chEff[8]*chEff[9]*chEff[10]; d1 *= d1;
1475  Double_t d2 = chEff[7]*chEff[9] + chEff[7]*chEff[10] + chEff[9]*chEff[10] - 3.*chEff[7]*chEff[9]*chEff[10]; d2 *= d2;
1476  Double_t d3 = chEff[7]*chEff[8] + chEff[7]*chEff[10] + chEff[8]*chEff[10] - 3.*chEff[7]*chEff[8]*chEff[10]; d3 *= d3;
1477  Double_t d4 = chEff[7]*chEff[8] + chEff[7]*chEff[9] + chEff[8]*chEff[9] - 3.*chEff[7]*chEff[8]*chEff[9]; d4 *= d4;
1478  Double_t d12 = chEff[9] + chEff[10] - 3.*chEff[9]*chEff[10]; d12 *= d12;
1479  Double_t d13 = chEff[8] + chEff[10] - 3.*chEff[8]*chEff[10]; d13 *= d13;
1480  Double_t d14 = chEff[8] + chEff[9] - 3.*chEff[8]*chEff[9]; d14 *= d14;
1481  Double_t d23 = chEff[7] + chEff[10] - 3.*chEff[7]*chEff[10]; d23 *= d23;
1482  Double_t d24 = chEff[7] + chEff[9] - 3.*chEff[7]*chEff[9]; d24 *= d24;
1483  Double_t d34 = chEff[7] + chEff[8] - 3.*chEff[7]*chEff[8]; d34 *= d34;
1484  Double_t d123 = 1. - 3.*chEff[10]; d123 *= d123;
1485  Double_t d124 = 1. - 3.*chEff[9]; d124 *= d124;
1486  Double_t d134 = 1. - 3.*chEff[8]; d134 *= d134;
1487  Double_t d234 = 1. - 3.*chEff[7]; d234 *= d234;
1488  Double_t d1234 = - 3.; d1234 *= d1234;
1489 
1490  for (Int_t i = 0; i < 2; ++i) {
1491  Double_t s1 = chEffErr[i][7] * chEffErr[i][7];
1492  Double_t s2 = chEffErr[i][8] * chEffErr[i][8];
1493  Double_t s3 = chEffErr[i][9] * chEffErr[i][9];
1494  Double_t s4 = chEffErr[i][10] * chEffErr[i][10];
1495  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);
1496  }
1497 
1498  st45EffErr[0] = TMath::Min(st45Eff, st45EffErr[0]);
1499  st45EffErr[1] = TMath::Min(1.-st45Eff, st45EffErr[1]);
1500 
1501 }
1502 
1503 
1504 //---------------------------------------------------------------------------
1505 void GetStation45Efficiency(TArrayD &chEff, TArrayD chEffErr[2], TGraphAsymmErrors *effVsX[3], Int_t ip, Double_t xp)
1506 {
1508 
1509  if (chEff[7] >= 0 && chEff[8] >= 0 && chEff[9] >= 0 && chEff[10] >= 0) {
1510 
1511  // compute station efficiency from known chamber efficiency
1512  Double_t stEff, stEffErr[2];
1513  ComputeStation45Efficiency(chEff, chEffErr, stEff, stEffErr);
1514 
1515  // fill graphs if required
1516  for (Int_t i = 0; i < 3; ++i) {
1517  if (effVsX[i]) {
1518  effVsX[i]->SetPoint(ip,xp,stEff);
1519  effVsX[i]->SetPointError(ip,0.,0.,stEffErr[0],stEffErr[1]);
1520  }
1521  }
1522 
1523  } else {
1524 
1525  Double_t edge[2] = {0., 1.};
1526  TArrayD chEffEdge[2];
1527  Double_t stEffEdge[2];
1528  Double_t stEffEdgeErr[2][2];
1529 
1530  for (Int_t i = 0; i < 2; ++i) {
1531 
1532  // set lower(upper) limit of chamber efficiency
1533  chEffEdge[i].Set(11);
1534  for (Int_t iCh = 7; iCh < 11; ++iCh) chEffEdge[i][iCh] = (chEff[iCh] < 0) ? edge[i] : chEff[iCh];
1535 
1536  // compute station efficiency
1537  ComputeStation45Efficiency(chEffEdge[i], chEffErr, stEffEdge[i], stEffEdgeErr[i]);
1538 
1539  // fill graph if required
1540  if (effVsX[i+1]) {
1541  effVsX[i+1]->SetPoint(ip,xp,stEffEdge[i]);
1542  effVsX[i+1]->SetPointError(ip,0.,0.,stEffEdgeErr[i][0],stEffEdgeErr[i][1]);
1543  }
1544 
1545  }
1546 
1547  // combine extreme cases to get station efficiency and fill graph if required
1548  if (effVsX[0]) {
1549  effVsX[0]->SetPoint(ip,xp,stEffEdge[1]);
1550  effVsX[0]->SetPointError(ip,0.,0.,stEffEdge[1]-stEffEdge[0]+stEffEdgeErr[0][0],stEffEdgeErr[1][1]);
1551  }
1552 
1553  }
1554 
1555 }
1556 
1557 
1558 //---------------------------------------------------------------------------
1559 void ComputeTrackingEfficiency(Double_t stEff[6], Double_t stEffErr[6][2], Double_t &spectroEff, Double_t spectroEffErr[2])
1560 {
1562 
1563  Double_t de[6][2];
1564  for (Int_t iSt = 0; iSt < 6; iSt++) de[iSt][0] = stEff[iSt]*stEff[iSt];
1565 
1566  if (moreTrackCandidates) {
1567 
1568  spectroEff = stEff[0] * stEff[1] * stEff[2] * stEff[3] * stEff[4];
1569 
1570  for (Int_t i = 0; i < 2; i++) {
1571 
1572  for (Int_t iSt = 0; iSt < 6; iSt++) de[iSt][1] = stEffErr[iSt][i]*stEffErr[iSt][i];
1573 
1574  spectroEffErr[i] = 0.;
1575  for (Int_t j = 1; j < 32; j++) {
1576  Double_t sigmaAdd = 1.;
1577  for (Int_t iSt = 0; iSt < 5; iSt++) sigmaAdd *= de[iSt][TESTBIT(j,iSt)];
1578  spectroEffErr[i] += sigmaAdd;
1579  }
1580  spectroEffErr[i] = TMath::Sqrt(spectroEffErr[i]);
1581 
1582  }
1583 
1584  } else {
1585 
1586  spectroEff = stEff[0] * stEff[1] * stEff[2] * stEff[5];
1587 
1588  for (Int_t i = 0; i < 2; i++) {
1589 
1590  for (Int_t iSt = 0; iSt < 6; iSt++) de[iSt][1] = stEffErr[iSt][i]*stEffErr[iSt][i];
1591 
1592  spectroEffErr[i] = 0.;
1593  for (Int_t j = 1; j < 16; j++) {
1594  Double_t sigmaAdd = de[5][TESTBIT(j,3)];
1595  for (Int_t iSt = 0; iSt < 3; iSt++) sigmaAdd *= de[iSt][TESTBIT(j,iSt)];
1596  spectroEffErr[i] += sigmaAdd;
1597  }
1598  spectroEffErr[i] = TMath::Sqrt(spectroEffErr[i]);
1599 
1600  }
1601 
1602  }
1603 
1604  spectroEffErr[0] = TMath::Min(spectroEff, spectroEffErr[0]);
1605  spectroEffErr[1] = TMath::Min(1.-spectroEff, spectroEffErr[1]);
1606 
1607 }
1608 
1609 
1610 //---------------------------------------------------------------------------
1611 void GetTrackingEfficiency(TArrayD &chEff, TArrayD chEffErr[2], TGraphAsymmErrors *effVsSt[3],
1612  TGraphAsymmErrors *effVsX[3], Int_t ip, Double_t xp, Bool_t print)
1613 {
1615 
1616  Double_t stEff[6];
1617  Double_t stEffErr[6][2];
1618 
1619  // check if unknown chamber efficiency
1620  Bool_t allEffKnown = kTRUE;
1621  for (Int_t iCh = 1; iCh < 11 && allEffKnown; ++iCh) if (chEff[iCh] < 0) allEffKnown = kFALSE;
1622 
1623  if (allEffKnown) {
1624 
1625  // compute station efficiency
1626  for (Int_t iSt = 0; iSt < 5; ++iSt) ComputeStationEfficiency(chEff, chEffErr, iSt, stEff[iSt], stEffErr[iSt]);
1627  ComputeStation45Efficiency(chEff, chEffErr, stEff[5], stEffErr[5]);
1628 
1629  // compute spectrometer efficiency
1630  Double_t spectroEff, spectroEffErr[2];
1631  ComputeTrackingEfficiency(stEff, stEffErr, spectroEff, spectroEffErr);
1632  chEff[0] = spectroEff;
1633  chEffErr[0][0] = spectroEffErr[0];
1634  chEffErr[1][0] = spectroEffErr[1];
1635 
1636  // fill graphs if required
1637  for (Int_t i = 0; i < 3; ++i) {
1638  if (effVsX[i]) {
1639  effVsX[i]->SetPoint(ip,xp,chEff[0]);
1640  effVsX[i]->SetPointError(ip,0.,0.,chEffErr[0][0],chEffErr[1][0]);
1641  }
1642  if (effVsSt[i]) for (Int_t iSt = 0; iSt < 6; ++iSt) {
1643  effVsSt[i]->SetPoint(iSt,iSt+1,stEff[iSt]);
1644  effVsSt[i]->SetPointError(iSt,0.,0.,stEffErr[iSt][0],stEffErr[iSt][1]);
1645  }
1646  }
1647 
1648  } else {
1649 
1650  Double_t edge[2] = {0., 1.};
1651  TArrayD chEffEdge[2];
1652  Double_t stEffEdge[2][6];
1653  Double_t stEffEdgeErr[2][6][2];
1654  Double_t spectroEffEdge[2], spectroEffEdgeErr[2][2];
1655 
1656  for (Int_t i = 0; i < 2; ++i) {
1657 
1658  // set lower(upper) limit of chamber efficiency
1659  chEffEdge[i].Set(11);
1660  for (Int_t iCh = 1; iCh < 11; ++iCh) chEffEdge[i][iCh] = (chEff[iCh] < 0) ? edge[i] : chEff[iCh];
1661 
1662  // compute station efficiency
1663  for (Int_t iSt = 0; iSt < 5; ++iSt) ComputeStationEfficiency(chEffEdge[i], chEffErr, iSt, stEffEdge[i][iSt], stEffEdgeErr[i][iSt]);
1664  ComputeStation45Efficiency(chEffEdge[i], chEffErr, stEffEdge[i][5], stEffEdgeErr[i][5]);
1665 
1666  // compute spectrometer efficiency
1667  ComputeTrackingEfficiency(stEffEdge[i], stEffEdgeErr[i], spectroEffEdge[i], spectroEffEdgeErr[i]);
1668 
1669  // fill graph if required
1670  if (effVsX[i+1]) {
1671  effVsX[i+1]->SetPoint(ip,xp,spectroEffEdge[i]);
1672  effVsX[i+1]->SetPointError(ip,0.,0.,spectroEffEdgeErr[i][0],spectroEffEdgeErr[i][1]);
1673  }
1674  if (effVsSt[i+1]) for (Int_t iSt = 0; iSt < 6; ++iSt) {
1675  effVsSt[i+1]->SetPoint(iSt,iSt+1,stEffEdge[i][iSt]);
1676  effVsSt[i+1]->SetPointError(iSt,0.,0.,stEffEdgeErr[i][iSt][0],stEffEdgeErr[i][iSt][1]);
1677  }
1678 
1679  }
1680 
1681  // combine extreme cases to get station efficiency
1682  for (Int_t iSt = 0; iSt < 6; ++iSt) {
1683  stEff[iSt] = stEffEdge[1][iSt];
1684  stEffErr[iSt][0] = stEffEdge[1][iSt] - stEffEdge[0][iSt] + stEffEdgeErr[0][iSt][0];
1685  stEffErr[iSt][1] = stEffEdgeErr[1][iSt][1];
1686  }
1687 
1688  // combine extreme cases to get spectrometer efficiency
1689  chEff[0] = spectroEffEdge[1];
1690  chEffErr[0][0] = spectroEffEdge[1] - spectroEffEdge[0] + spectroEffEdgeErr[0][0];
1691  chEffErr[1][0] = spectroEffEdgeErr[1][1];
1692 
1693  // fill graph if required
1694  if (effVsX[0]) {
1695  effVsX[0]->SetPoint(ip,xp,chEff[0]);
1696  effVsX[0]->SetPointError(ip,0.,0.,chEffErr[0][0],chEffErr[1][0]);
1697  }
1698  if (effVsSt[0]) for (Int_t iSt = 0; iSt < 6; ++iSt) {
1699  effVsSt[0]->SetPoint(iSt,iSt+1,stEff[iSt]);
1700  effVsSt[0]->SetPointError(iSt,0.,0.,stEffErr[iSt][0],stEffErr[iSt][1]);
1701  }
1702 
1703  }
1704 
1705  // print results
1706  if (print) {
1707  for (Int_t iCh = 1; iCh < 11; ++iCh) {
1708  cout << "Efficiency chamber " << iCh << " : ";
1709  cout << chEff[iCh] << " + " << chEffErr[1][iCh] << " - " << chEffErr[0][iCh] << endl;
1710  }
1711  for (Int_t iSt = 0; iSt < 6; ++iSt) {
1712  if (iSt < 5) cout << "Station " << iSt+1 << " = ";
1713  else cout << "Station 45 = ";
1714  cout << stEff[iSt] << " + " << stEffErr[iSt][1] << " - " << stEffErr[iSt][0] << endl;
1715  }
1716  cout << "Total tracking efficiency : ";
1717  cout << chEff[0] << " + " << chEffErr[1][0] << " - " << chEffErr[0][0] << endl << endl;
1718  }
1719 
1720 }
1721 
1722 
1723 //---------------------------------------------------------------------------
1724 void IntegrateMuonEfficiency(TGraphAsymmErrors &effVsRunLow, TGraphAsymmErrors &effVsRunUp,
1725  TGraphAsymmErrors &effVsX, Int_t ip, Double_t xp)
1726 {
1729 
1730  if (!runWeights) {
1731  effVsX.SetPoint(ip,xp,-1.);
1732  effVsX.SetPointError(ip,0.,0.,0.,0.);
1733  return;
1734  }
1735 
1736  // initialize
1737  Double_t rec[2] = {0., 0.};
1738  Double_t gen = 0.;
1739  Double_t intEffErr[2] = {0., 0.};
1740  Bool_t ok = kFALSE;
1741 
1742  // loop over runs
1743  Int_t nRuns = effVsRunLow.GetN();
1744  for (Int_t iRun = 0; iRun < nRuns; ++iRun) {
1745 
1746  // get run weight
1747  TString sRun = effVsRunLow.GetXaxis()->GetBinLabel(iRun+1);
1748  TParameter<Double_t> *weight = static_cast<TParameter<Double_t>*>(runWeights->FindObject(sRun.Data()));
1749  if (!weight) {
1750  printf("weight not found for run %s\n", sRun.Data());
1751  continue;
1752  }
1753  Double_t w = weight->GetVal();
1754  Double_t w2 = w*w;
1755 
1756  // get efficiency and error
1757  Double_t eff[2] = {effVsRunLow.GetY()[iRun], effVsRunUp.GetY()[iRun]};
1758  Double_t effErr[2] = {effVsRunLow.GetErrorYlow(iRun), effVsRunUp.GetErrorYhigh(iRun)};
1759  if (eff[0] < 0. || eff[1] < 0.) {
1760  printf("no efficiency measurement --> use 0(1) ± 0 as lower(upper) limit\n");
1761  eff[0] = 0.;
1762  eff[1] = 1.;
1763  effErr[0] = 0.;
1764  effErr[1] = 0.;
1765  } else ok = kTRUE;
1766 
1767  // integrate
1768  gen += w;
1769  for (Int_t i = 0; i < 2; ++i) {
1770  rec[i] += w*eff[i];
1771  intEffErr[i] += w2*effErr[i]*effErr[i];
1772  }
1773 
1774  }
1775 
1776  // compute efficiency
1777  if (gen > 0. && ok) {
1778 
1779  effVsX.SetPoint(ip,xp,rec[1]/gen);
1780  effVsX.SetPointError(ip,0.,0.,(rec[1]-rec[0]+TMath::Sqrt(intEffErr[0]))/gen,TMath::Sqrt(intEffErr[1])/gen);
1781 
1782  } else {
1783 
1784  if (gen <= 0.) printf("impossible to integrate, all weights = 0 or unknown ?!?\n");
1785  else printf("efficiency never measured --> return -1 ± 0\n");
1786 
1787  effVsX.SetPoint(ip,xp,-1.);
1788  effVsX.SetPointError(ip,0.,0.,0.,0.);
1789 
1790  }
1791 
1792 }
1793 
1794 
1795 //---------------------------------------------------------------------------
1797 {
1801 
1802  if (runWeights) return;
1803 
1804  ifstream inFile(fileName.Data());
1805  if (!inFile.is_open()) {
1806  printf("cannot open file %s\n", fileName.Data());
1807  return;
1808  }
1809 
1810  runWeights = new THashList(1000);
1811  runWeights->SetOwner();
1812 
1813  TString line;
1814  while (! inFile.eof() ) {
1815 
1816  line.ReadLine(inFile,kTRUE);
1817  if(line.IsNull()) continue;
1818 
1819  TObjArray *param = line.Tokenize(" ");
1820  if (param->GetEntries() != 2) {
1821  printf("bad input line %s", line.Data());
1822  continue;
1823  }
1824 
1825  Int_t run = ((TObjString*)param->UncheckedAt(0))->String().Atoi();
1826  if (run < 0) {
1827  printf("invalid run number: %d", run);
1828  continue;
1829  }
1830 
1831  Float_t weight = ((TObjString*)param->UncheckedAt(1))->String().Atof();
1832  if (weight <= 0.) {
1833  printf("invalid weight: %g", weight);
1834  continue;
1835  }
1836 
1837  runWeights->Add(new TParameter<Double_t>(Form("%d",run), weight));
1838 
1839  delete param;
1840  }
1841 
1842  inFile.close();
1843 
1844 }
1845 
1846 
1847 //---------------------------------------------------------------------------
1848 void SetCentPtCh(THnSparse& SparseData)
1849 {
1852 
1853  if (ptMax >= 0 && ptMax < ptMin) {
1854  printf("Wrong Pt range, ptMax must be higher than ptMin\n");
1855  exit(-1);
1856  }
1857 
1858  if (charge < -1 || charge > 1) {
1859  printf("Selected charge must be 0, 1 or -1\n");
1860  exit(-1);
1861  }
1862 
1863  // adjust centrality range
1864  centMax = TMath::Max(centMax-1.e-12, centMin);
1865 
1866  // set the centrality range for integration
1867  Int_t lowBin = SparseData.GetAxis(1)->FindBin(centMin);
1868  Int_t upBin = SparseData.GetAxis(1)->FindBin(centMax);
1869  SparseData.GetAxis(1)->SetRange(lowBin, upBin);
1870 
1871  // set the pt range for integration
1872  lowBin = SparseData.GetAxis(2)->FindBin(ptMin);
1873  if (ptMax == -1) { upBin = SparseData.GetAxis(2)->GetNbins()+1; }
1874  else { upBin = SparseData.GetAxis(2)->FindBin(ptMax);}
1875  SparseData.GetAxis(2)->SetRange(lowBin, upBin);
1876 
1877  // set the charge range
1878  if (charge != 0) {
1879  lowBin = SparseData.GetAxis(5)->FindBin(charge);
1880  SparseData.GetAxis(5)->SetRange(lowBin, lowBin);
1881  }
1882 
1883 }
1884 
1885 
1886 //---------------------------------------------------------------------------
1887 TGraphAsymmErrors* CreateGraph(const char* name, const char* title, int value)
1888 {
1890 
1891  TGraphAsymmErrors* g = new TGraphAsymmErrors();
1892 
1893  if ( value >= 0 ) {
1894 
1895  g->SetName(Form(name,value));
1896  g->SetTitle(Form(title,value));
1897 
1898  } else {
1899 
1900  g->SetName(name);
1901  g->SetTitle(title);
1902 
1903 
1904  }
1905 
1906  return g;
1907 
1908 }
1909 
1910 
1911 //---------------------------------------------------------------------------
1912 void BeautifyGraph(TGraphAsymmErrors &g, const char* xAxisName, const char* yAxisName)
1913 {
1915 
1916  g.SetLineStyle(1);
1917  g.SetLineColor(1);
1918  g.SetMarkerStyle(20);
1919  g.SetMarkerSize(0.7);
1920  g.SetMarkerColor(2);
1921 
1922  g.GetXaxis()->SetTitle(xAxisName);
1923  g.GetXaxis()->CenterTitle(kTRUE);
1924  g.GetXaxis()->SetLabelFont(22);
1925  g.GetXaxis()->SetTitleFont(22);
1926 
1927  g.GetYaxis()->SetTitle(yAxisName);
1928  g.GetYaxis()->CenterTitle(kTRUE);
1929  g.GetYaxis()->SetLabelFont(22);
1930  g.GetYaxis()->SetTitleFont(22);
1931 
1932 }
1933 
1934 
1935 //---------------------------------------------------------------------------
1936 void BeautifyGraphs(TObjArray& array, const char* xAxisName, const char* yAxisName)
1937 {
1939 
1940  for ( Int_t i = 0; i <= array.GetLast(); ++i)
1941  BeautifyGraph(*static_cast<TGraphAsymmErrors*>(array.UncheckedAt(i)), xAxisName, yAxisName);
1942 
1943 }
1944 
1945 
1946 //---------------------------------------------------------------------------
1947 void SetRunLabel(TGraphAsymmErrors &g, Int_t irun, const TList& runs)
1948 {
1950 
1951  g.GetXaxis()->Set(irun+1, -0.5, irun+0.5);
1952 
1953  TIter nextRun(&runs);
1954  TObjString *srun = 0x0;
1955  Int_t iirun = 1;
1956  while ((srun = static_cast<TObjString*>(nextRun()))) {
1957  g.GetXaxis()->SetBinLabel(iirun, srun->GetName());
1958  iirun++;
1959  }
1960 
1961 }
1962 
1963 
1964 //---------------------------------------------------------------------------
1965 void SetRunLabel(TObjArray& array, Int_t irun, const TList& runs)
1966 {
1968 
1969  for (Int_t i = 0; i <= array.GetLast(); ++i)
1970  SetRunLabel(*static_cast<TGraphAsymmErrors*>(array.UncheckedAt(i)), irun, runs);
1971 
1972 }
1973 
1974 
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")
void GetStationEfficiency(TArrayD &chEff, TArrayD chEffErr[2], Int_t iSt, TGraphAsymmErrors *effVsX[3], Int_t ip, Double_t xp)
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
void PlotMuonEfficiency(TString fileNameData, TString fileNameSave, Bool_t saveEdges, Bool_t print, Bool_t draw)
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
void SetRunLabel(TGraphAsymmErrors &g, Int_t irun, const TList &runs)
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)
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
void PlotMuonEfficiencyVsX(TString var, TString fileNameData, TString fileNameSave, Bool_t saveEdges, Bool_t print, Bool_t draw)
TFile * file
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