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