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