-#include <fstream>
-#include <stdio.h>
-#include <stdlib.h>
-#include <TH1I.h>
-#include <TH2I.h>
-#include <TH1D.h>
-#include <TH2D.h>
+/**
+ *
+ * \file MuonTrackingEfficiency.C
+ * \author Philippe Pillot, Antoine Lardeux, Lizardo Valencia Palomo, Javier Martin Blanco
+ * \brief Compute trk efficiency at DE, chamber, station and spectro levels vs various variables from the output of the efficiency task
+ *
+ * 1) efficiency estimator and error calculation at chamber and DE level:
+ * 2 options:
+ * - using bayesian method with uniform prior
+ * - using Clopper-Pearson or other frequentist method
+ * if n ≠ 0: use above methods
+ * if n = 0: eff = -1 ± 0
+ *
+ * 2) efficiency and error propagation at station and spectrometer level:
+ * if eff = -1 for one or several ch/DE:
+ * - assume eff_ch = 1 ± 0 to compute eff_up and err_up with std error propagation at nth order
+ * - assume eff_ch = 0 ± 0 to compute eff_low and err_low with std error propagation at nth order
+ * - eff_spectro = eff_up + err_up - (eff_up-eff_low + err_low)
+ * otherwise: std efficiency and error propagation at nth order
+ *
+ * 3) efficiency and error integration over runs
+ * - performed from efficiency plots versus run
+ * - compute weighted average and apply standard error propagation to err_up and err_low
+ * - do it for both extreme cases above (with eff_up ± err–up and eff_low ± err_low)
+ * - int_eff = int_eff_up + int_err_up - (int_eff_up-int_eff_low + int_err_low)
+ *
+ */
+
+#include <TROOT.h>
+#include <TH1.h>
+#include <TH2.h>
+#include <THnSparse.h>
#include <TAxis.h>
-#include <TLegend.h>
#include <TString.h>
+#include <TObjString.h>
#include <Riostream.h>
#include <TFile.h>
#include <TList.h>
#include <TCanvas.h>
-#include <TSystem.h>
-#include <TStyle.h>
#include <TGraphAsymmErrors.h>
-#include <vector>
-#include <TGaxis.h>
+#include <TMath.h>
+#include <TArrayD.h>
+#include <TStyle.h>
+#include <THashList.h>
+#include <TParameter.h>
+
+//const Char_t *effErrMode = "cp"; // Clopper-Pearson
+const Char_t *effErrMode = "b(1,1)mode"; // Bayesian with uniform prior
+
+Double_t centMin = -999.;
+Double_t centMax = 999.;
+Double_t ptMin = 0.;
+Double_t ptMax = -1.;
+Int_t charge = 0; // 0 selects + and -, -1 and +1 selects - or + muons respectively
+
+Bool_t moreTrackCandidates = kFALSE;
+
+THashList *runWeights = 0x0;
+
+
+void PlotMuonEfficiencyVsX(TString var, TString fileNameData, TString fileNameSave, Bool_t saveEdges, Bool_t print, Bool_t draw);
+void PlotIntegratedMuonEfficiencyVsX(TString var, TString runList, TString fileNameWeights,
+ TString fileNameData, 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);
+
+void PlotMuonEfficiency(TString fileNameData, TString fileNameSave, Bool_t saveEdges, Bool_t print, Bool_t draw);
+void PlotMuonEfficiencyVsRun(TString runList, TString fileNameData, TString fileNameSave, Bool_t print, Bool_t draw);
+void PlotIntegratedMuonEfficiency(TString fileNameWeights, TString fileNameSave, Bool_t print, Bool_t draw);
+
+void PlotMuonEfficiencyPerDE(TString fileNameData, TString fileNameSave, Bool_t saveEdges);
+void PlotMuonEfficiencyPerDEVsRun(TString runList, TString fileNameData, TString fileNameSave);
+void PlotIntegratedMuonEfficiencyPerDE(TString fileNameWeights, TString fileNameSave);
+
+Bool_t GetChamberEfficiency(THnSparse &TT, THnSparse &TD, TArrayD &chEff, TArrayD chEffErr[2], Bool_t printError = kFALSE);
+void GetDEEfficiency(THnSparse &TT, THnSparse &TD, TGraphAsymmErrors &effVsDE);
+void ComputeStationEfficiency(TArrayD &chEff, TArrayD chEffErr[2], Int_t iSt, Double_t &stEff, Double_t stEffErr[2]);
+void GetStationEfficiency(TArrayD &chEff, TArrayD chEffErr[2], Int_t iSt, TGraphAsymmErrors *effVsX[3], Int_t ip, Double_t xp);
+void ComputeStation45Efficiency(TArrayD &chEff, TArrayD chEffErr[2], Double_t &st45Eff, Double_t st45EffErr[2]);
+void GetStation45Efficiency(TArrayD &chEff, TArrayD chEffErr[2], TGraphAsymmErrors *effVsX[3], Int_t ip, Double_t xp);
+void ComputeTrackingEfficiency(Double_t stEff[6], Double_t stEffErr[6][2], Double_t &spectroEff, Double_t spectroEffErr[2]);
+void GetTrackingEfficiency(TArrayD &chEff, TArrayD chEffErr[2], TGraphAsymmErrors *effVsSt[3],
+ TGraphAsymmErrors *effVsX[3], Int_t ip, Double_t xp, Bool_t print = kFALSE);
+void IntegrateMuonEfficiency(TGraphAsymmErrors &effVsRunLow, TGraphAsymmErrors &effVsRunUp,
+ TGraphAsymmErrors &effVsX, Int_t ip, Double_t xp);
+
+void LoadRunWeights(TString fileName);
+void SetCentPtCh(THnSparse& SparseData);
+TGraphAsymmErrors* CreateGraph(const char* name, const char* title, int value=-1);
+void BeautifyGraph(TGraphAsymmErrors &g, const char* xAxisName, const char* yAxisName);
+void BeautifyGraphs(TObjArray& array, const char* xAxisName, const char* yAxisName);
+void SetRunLabel(TGraphAsymmErrors &g, Int_t irun, const TList& runs);
+void SetRunLabel(TObjArray& array, Int_t irun, const TList& runs);
+
+
+//---------------------------------------------------------------------------
+void MuonTrackingEfficiency(TString runList = "runList.txt",
+ TString fileNameWeights = "",
+ TString fileNameData ="AnalysisResults.root",
+ TString fileNameSave = "efficiency_new.root")
+{
+ /// main function to compute, print and plot efficiencies
+
+ PlotMuonEfficiencyVsX("centrality", fileNameData, fileNameSave, kFALSE, kFALSE, kTRUE);
+ PlotMuonEfficiencyVsX("pt", fileNameData, fileNameSave, kFALSE, kFALSE, kTRUE);
+ PlotMuonEfficiencyVsX("y", fileNameData, fileNameSave, kFALSE, kFALSE, kTRUE);
+ PlotMuonEfficiencyVsX("phi", fileNameData, fileNameSave, kFALSE, kFALSE, kTRUE);
+ PlotMuonEfficiencyVsX("charge", fileNameData, fileNameSave, kFALSE, kFALSE, kTRUE);
+
+ PlotIntegratedMuonEfficiencyVsX("centrality", runList, fileNameWeights, fileNameData, fileNameSave, kFALSE, kTRUE);
+ PlotIntegratedMuonEfficiencyVsX("pt", runList, fileNameWeights, fileNameData, fileNameSave, kFALSE, kTRUE);
+ PlotIntegratedMuonEfficiencyVsX("y", runList, fileNameWeights, fileNameData, fileNameSave, kFALSE, kTRUE);
+ PlotIntegratedMuonEfficiencyVsX("phi", runList, fileNameWeights, fileNameData, fileNameSave, kFALSE, kTRUE);
+ PlotIntegratedMuonEfficiencyVsX("charge", runList, fileNameWeights, fileNameData, fileNameSave, kFALSE, kTRUE);
+
+ PlotMuonEfficiencyVsXY("pt", "centrality", fileNameData, fileNameSave, kTRUE);
+ PlotMuonEfficiencyVsXY("y", "centrality", fileNameData, fileNameSave, kTRUE);
+ PlotMuonEfficiencyVsXY("pt", "y", fileNameData, fileNameSave, kTRUE);
+ PlotMuonEfficiencyVsXY("phi", "y", fileNameData, fileNameSave, kTRUE, kTRUE);
+
+ PlotMuonEfficiency(fileNameData, fileNameSave, kFALSE, kTRUE, kTRUE);
+ PlotMuonEfficiencyVsRun(runList, fileNameData, fileNameSave, kFALSE, kTRUE);
+ PlotIntegratedMuonEfficiency(fileNameWeights, fileNameSave, kTRUE, kTRUE);
+
+ PlotMuonEfficiencyPerDE(fileNameData, fileNameSave, kFALSE);
+ PlotMuonEfficiencyPerDEVsRun(runList, fileNameData, fileNameSave);
+ PlotIntegratedMuonEfficiencyPerDE(fileNameWeights, fileNameSave);
+
+}
+
+
+//---------------------------------------------------------------------------
+void PlotMuonEfficiencyVsX(TString var, TString fileNameData, TString fileNameSave, Bool_t saveEdges, Bool_t print, Bool_t draw)
+{
+ /// plot the tracking efficiency versus X
+
+ printf("plotting efficiency versus %s...\n", var.Data());
+
+ Int_t xDim = -1;
+ if (var == "centrality") xDim = 1;
+ else if (var == "pt") xDim = 2;
+ else if (var == "y") xDim = 3;
+ else if (var == "phi") xDim = 4;
+ else if (var == "charge") xDim = 5;
+ else {
+ printf("incorrect variable. Choices are centrality, pt, y, phi and charge.\n");
+ return;
+ }
+
+ // get input hists
+ TFile *file = new TFile(fileNameData.Data(), "read");
+ if (!file || !file->IsOpen()) {
+ printf("cannot open file %s \n",fileNameData.Data());
+ return;
+ }
+ TList *listTT = static_cast<TList*>(file->FindObjectAny("TotalTracksPerChamber"));
+ TList *listTD = static_cast<TList*>(file->FindObjectAny("TracksDetectedPerChamber"));
+ THnSparse *TT = static_cast<THnSparse*>(listTT->At(10));
+ THnSparse *TD = static_cast<THnSparse*>(listTD->At(10));
+
+ // output graph
+ TGraphAsymmErrors *effVsX[3] = {0x0, 0x0, 0x0};
+ TString nameAdd[3] = {"", "Low", "Up"};
+ TString titleAdd[3] = {"", " - lower limit", " - upper limit"};
+ for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i) {
+ effVsX[i] = CreateGraph(Form("trackingEffVs%s%s",var.Data(),nameAdd[i].Data()),
+ Form("Measured tracking efficiency versus %s%s",var.Data(),titleAdd[i].Data()));
+ }
+
+ // set the centrality and pT range for integration
+ SetCentPtCh(*TT);
+ SetCentPtCh(*TD);
+
+ // loop over X bins
+ TArrayD chEff(11);
+ TArrayD chEffErr[2];
+ chEffErr[0].Set(11);
+ chEffErr[1].Set(11);
+ for (Int_t ix = 1; ix <= TT->GetAxis(xDim)->GetNbins(); ix++) {
+
+ if (print) cout << var.Data() << " " << TT->GetAxis(xDim)->GetBinLowEdge(ix) << "-" << TT->GetAxis(xDim)->GetBinUpEdge(ix) << ":" << endl;
+
+ // set the var range to the current bin
+ TT->GetAxis(xDim)->SetRange(ix, ix);
+ TD->GetAxis(xDim)->SetRange(ix, ix);
+
+ // compute chamber efficency and errors
+ if (GetChamberEfficiency(*TT, *TD, chEff, chEffErr, print)) {
+
+ // compute the overall tracking efficiency
+ TGraphAsymmErrors *dummy[3] = {0x0, 0x0, 0x0};
+ GetTrackingEfficiency(chEff, chEffErr, dummy, effVsX, ix-1, TT->GetAxis(xDim)->GetBinCenter(ix), print);
+
+ } else {
+
+ // fill graph with 1 +0 -1
+ effVsX[0]->SetPoint(ix-1,TT->GetAxis(xDim)->GetBinCenter(ix),1.);
+ effVsX[0]->SetPointError(ix-1,0.,0.,1.,0.);
+
+ if (saveEdges) {
+
+ // lower = 0 ± 0
+ effVsX[1]->SetPoint(ix-1,TT->GetAxis(xDim)->GetBinCenter(ix),0.);
+ effVsX[1]->SetPointError(ix-1,0.,0.,0.,0.);
+
+ // upper = 1 ± 0
+ effVsX[2]->SetPoint(ix-1,TT->GetAxis(xDim)->GetBinCenter(ix),1.);
+ effVsX[2]->SetPointError(ix-1,0.,0.,0.,0.);
+
+ }
+
+ }
+
+ }
+
+ // close input file
+ file->Close();
+
+ // display
+ for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i) BeautifyGraph(*effVsX[i], var.Data(), "efficiency");
+
+ if (draw) {
+ new TCanvas(Form("cTrackingEffVs%s",var.Data()), Form("Measured tracking efficiency versus %s",var.Data()),1000,400);
+ effVsX[0]->DrawClone("ap");
+ }
+
+ // save output
+ file = new TFile(fileNameSave.Data(),"update");
+ for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i) effVsX[i]->Write(0x0, TObject::kOverwrite);
+ file->Close();
+
+ // clean memory
+ for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i) delete effVsX[i];
+
+}
+
+
+//---------------------------------------------------------------------------
+void PlotIntegratedMuonEfficiencyVsX(TString var, TString runList, TString fileNameWeights,
+ TString fileNameData, TString fileNameSave, Bool_t print, Bool_t draw)
+{
+ /// plot the tracking efficiency versus X for each run and integrated
+
+ printf("plotting integrated efficiency versus %s...\n", var.Data());
+
+ // load run weights
+ if (!fileNameWeights.IsNull()) LoadRunWeights(fileNameWeights);
+ if (!runWeights) {
+ printf("Cannot compute integrated efficiency without run-by-run weights\n");
+ return;
+ }
+
+ // open run list
+ ifstream inFile(runList.Data());
+ if (!inFile.is_open()) {
+ printf("cannot open file %s\n",runList.Data());
+ return;
+ }
+
+ // output graph
+ TGraphAsymmErrors *intEffVsX = CreateGraph(Form("integratedTrackingEffVs%s",var.Data()),
+ Form("Integrated tracking efficiency versus %s",var.Data()));
+
+ Int_t n = -1;
+ TArrayD x;
+ TArrayD rec[2];
+ Double_t gen = 0.;
+ TArrayD effErr[2];
+
+ // loop over runs
+ while (!inFile.eof()) {
+
+ // get current run number
+ TString currRun;
+ currRun.ReadLine(inFile,kTRUE);
+ if(currRun.IsNull() || !currRun.IsDec()) continue;
+ Int_t run = currRun.Atoi();
+
+ printf("run %d: ", run);
+
+ // compute efficiency vs var
+ TString dataFile = Form("runs/%d/%s", run, fileNameData.Data());
+ TString outFile = Form("runs/%d/%s", run, fileNameSave.Data());
+ PlotMuonEfficiencyVsX(var, dataFile, outFile, kTRUE, print, kFALSE);
+
+ // get run weight
+ TParameter<Double_t> *weight = static_cast<TParameter<Double_t>*>(runWeights->FindObject(currRun.Data()));
+ if (!weight) {
+ printf("weight not found for run %s\n", currRun.Data());
+ continue;
+ }
+ Double_t w = weight->GetVal();
+ Double_t w2 = w*w;
+
+ // get result
+ TFile *file = new TFile(outFile.Data(), "read");
+ if (!file || !file->IsOpen()) {
+ printf("cannot open file\n");
+ continue;
+ }
+ TGraphAsymmErrors *effVsX[2];
+ effVsX[0] = static_cast<TGraphAsymmErrors*>(file->FindObjectAny(Form("trackingEffVs%sLow",var.Data())));
+ effVsX[1] = static_cast<TGraphAsymmErrors*>(file->FindObjectAny(Form("trackingEffVs%sUp",var.Data())));
+ if (!effVsX[0] || !effVsX[1]) {
+ printf("trackingEffVs%sLow(Up) object not found\n", var.Data());
+ continue;
+ }
+
+ // prepare the arrays if not already done
+ if (n < 0) {
+ n = effVsX[0]->GetN();
+ x.Set(n, effVsX[0]->GetX());
+ for (Int_t i = 0; i < 2; ++i) {
+ rec[i].Set(n);
+ effErr[i].Set(n);
+ }
+ } else if (n != effVsX[0]->GetN()) {
+ printf("number of points in graph trackingEffVs%sLow(Up) for run %d is different than from previous runs\n", var.Data(), run);
+ continue;
+ }
+
+ // integrate for all bins
+ gen += w;
+ for (Int_t ix = 0; ix < n; ++ix) {
+ Double_t ieffErr[2] = {effVsX[0]->GetErrorYlow(ix), effVsX[1]->GetErrorYhigh(ix)};
+ for (Int_t i = 0; i < 2; ++i) {
+ rec[i][ix] += w*effVsX[i]->GetY()[ix];
+ effErr[i][ix] += w2*ieffErr[i]*ieffErr[i];
+ }
+ }
+
+ // close input file
+ file->Close();
+
+ }
+
+ inFile.close();
+
+ // fill output graph
+ if (gen > 0.) {
+
+ for (Int_t ix = 0; ix < n; ++ix) {
+
+ intEffVsX->SetPoint(ix, x[ix], rec[1][ix]/gen);
+ intEffVsX->SetPointError(ix, 0., 0., (rec[1][ix]-rec[0][ix]+TMath::Sqrt(effErr[0][ix]))/gen, TMath::Sqrt(effErr[1][ix])/gen);
+
+ }
+
+ } else {
+
+ for (Int_t ix = 0; ix < n; ++ix) {
+
+ printf("impossible to integrate, all weights = 0 or unknown ?!?\n");
+
+ intEffVsX->SetPoint(ix, x[ix], -1.);
+ intEffVsX->SetPointError(ix, 0., 0., 0., 0.);
+
+ }
+
+ }
+
+ // display
+ BeautifyGraph(*intEffVsX, var.Data(), "efficiency");
+
+ if (draw) {
+ new TCanvas(Form("cIntegratedTrackingEffVs%s",var.Data()), Form("Integrated tracking efficiency versus %s",var.Data()),1000,400);
+ intEffVsX->DrawClone("ap");
+ }
+
+ // save output
+ TFile *file = new TFile(fileNameSave.Data(),"update");
+ intEffVsX->Write(0x0, TObject::kOverwrite);
+ file->Close();
+
+ // clean memory
+ delete intEffVsX;
+
+}
+
+
+//---------------------------------------------------------------------------
+void PlotMuonEfficiencyVsXY(TString xVar, TString yVar, TString fileNameData, TString fileNameSave, Bool_t draw, Bool_t rap)
+{
+ /// plot the tracking efficiency versus X,Y
+
+ printf("plotting efficiency versus %s/%s...\n", xVar.Data(), yVar.Data());
+
+ Int_t xDim = -1;
+ if (xVar == "centrality") xDim = 1;
+ else if (xVar == "pt") xDim = 2;
+ else if (xVar == "y") xDim = 3;
+ else if (xVar == "phi") xDim = 4;
+ else if (xVar == "charge") xDim = 5;
+ else {
+ printf("incorrect variable. Choices are centrality, pt, y, phi and charge.\n");
+ return;
+ }
+ Int_t yDim = -1;
+ if (yVar == "centrality") yDim = 1;
+ else if (yVar == "pt") yDim = 2;
+ else if (yVar == "y") yDim = 3;
+ else if (yVar == "phi") yDim = 4;
+ else if (yVar == "charge") yDim = 5;
+ else {
+ printf("incorrect variable. Choices are centrality, pt, y, phi and charge.\n");
+ return;
+ }
+
+ // get input hists
+ TFile *file = new TFile(fileNameData.Data(), "read");
+ if (!file || !file->IsOpen()) {
+ printf("cannot open file %s\n",fileNameData.Data());
+ return;
+ }
+ TList *listTT = static_cast<TList*>(file->FindObjectAny("TotalTracksPerChamber"));
+ TList *listTD = static_cast<TList*>(file->FindObjectAny("TracksDetectedPerChamber"));
+ THnSparse *TT = static_cast<THnSparse*>(listTT->At(10));
+ THnSparse *TD = static_cast<THnSparse*>(listTD->At(10));
+
+ // output map
+ Int_t nxBins = TT->GetAxis(xDim)->GetNbins();
+ Int_t nyBins = TT->GetAxis(yDim)->GetNbins();
+ TH2F *effVsXY = new TH2F(Form("trackingEffVs%s-%s",xVar.Data(),yVar.Data()),
+ Form("Measured tracking efficiency versus %s and %s",xVar.Data(),yVar.Data()),
+ nxBins, TT->GetAxis(xDim)->GetBinLowEdge(1), TT->GetAxis(xDim)->GetBinUpEdge(nxBins),
+ nyBins, TT->GetAxis(yDim)->GetBinLowEdge(1), TT->GetAxis(yDim)->GetBinUpEdge(nyBins));
+ effVsXY->SetDirectory(0);
+
+ // set the centrality and pT range for integration
+ SetCentPtCh(*TT);
+ SetCentPtCh(*TD);
+
+ // loop over X/Y bins
+ TArrayD chEff(11);
+ TArrayD chEffErr[2];
+ chEffErr[0].Set(11);
+ chEffErr[1].Set(11);
+ for (Int_t ix = 1; ix <= nxBins; ++ix) {
+
+ // set x range
+ TT->GetAxis(xDim)->SetRange(ix, ix);
+ TD->GetAxis(xDim)->SetRange(ix, ix);
+
+ for (Int_t iy = 1; iy <= nyBins; ++iy) {
+
+ // set y range
+ TT->GetAxis(yDim)->SetRange(iy, iy);
+ TD->GetAxis(yDim)->SetRange(iy, iy);
+
+ // compute chamber efficency and errors
+ if (GetChamberEfficiency(*TT, *TD, chEff, chEffErr, kFALSE)) {
+
+ // compute the overall tracking efficiency
+ TGraphAsymmErrors *dummy[3] = {0x0, 0x0, 0x0};
+ GetTrackingEfficiency(chEff, chEffErr, dummy, dummy, 0, 0.);
+
+ // fill histo
+ effVsXY->Fill(TT->GetAxis(xDim)->GetBinCenter(ix),TT->GetAxis(yDim)->GetBinCenter(iy),chEff[0]);
+ effVsXY->SetBinError(ix,iy,TMath::Max(chEffErr[0][0], chEffErr[1][0]));
+
+ } else {
+
+ // fill histo with 0 ± 1
+ effVsXY->Fill(TT->GetAxis(xDim)->GetBinCenter(ix),TT->GetAxis(yDim)->GetBinCenter(iy),0.);
+ effVsXY->SetBinError(ix,iy,1.);
+
+ }
+
+ }
+
+ }
+
+ // close input file
+ file->Close();
+
+ // display
+ effVsXY->GetXaxis()->SetTitle(xVar.Data());
+ effVsXY->GetXaxis()->CenterTitle(kTRUE);
+ effVsXY->GetXaxis()->SetLabelFont(22);
+ effVsXY->GetXaxis()->SetTitleFont(22);
+ effVsXY->GetYaxis()->SetTitle(yVar.Data());
+ effVsXY->GetYaxis()->CenterTitle(kTRUE);
+ effVsXY->GetYaxis()->SetLabelFont(22);
+ effVsXY->GetYaxis()->SetTitleFont(22);
+ effVsXY->GetZaxis()->SetTitle("efficiency");
+ effVsXY->GetZaxis()->SetLabelFont(22);
+ effVsXY->GetZaxis()->SetTitleFont(22);
+
+ if (draw) {
+ new TCanvas(Form("cTrackingEffVs%s-%s",xVar.Data(),yVar.Data()), Form("Measured tracking efficiency versus %s and %s",xVar.Data(),yVar.Data()),700,600);
+ effVsXY->DrawClone("surf1");
+ }
+
+ // save output
+ file = new TFile(fileNameSave.Data(),"update");
+ effVsXY->Write(0x0, TObject::kOverwrite);
+
+ // add an histo with variable size rapidity bins
+ if (yDim == 3 && rap) {
+
+ TH2F* effVsXYrap = new TH2F();
+ TString rapName = Form("trackingEffVs%s-%sRapBins", xVar.Data(), yVar.Data());
+ TString rapTitle = Form("Measured tracking efficiency versus %s and %s", xVar.Data(), yVar.Data());
+ effVsXYrap->SetTitle(rapTitle.Data());
+ effVsXYrap->SetName(rapName.Data());
+
+ Double_t xBinEdge[nxBins+1];
+ for (Int_t xbin = 0; xbin <= nxBins; ++xbin)
+ xBinEdge[xbin] = effVsXY->GetXaxis()->GetBinLowEdge(xbin+1);
+
+ Double_t yBinEdge[nyBins+1];
+ for (Int_t ybin = 0; ybin <= nyBins; ++ybin)
+ yBinEdge[ybin] = 2*TMath::ATan(TMath::Exp((effVsXY->GetYaxis()->GetBinLowEdge(ybin+1))));
+
+ effVsXYrap->SetBins(nxBins, xBinEdge, nyBins, yBinEdge);
+
+ for (Int_t xbin = 1; xbin <= nxBins; ++xbin)
+ for (Int_t ybin = 1; ybin <= nyBins; ++ybin)
+ effVsXYrap->SetBinContent(xbin, ybin, effVsXY->GetBinContent(xbin,ybin));
+
+ effVsXYrap->Write(0x0, TObject::kOverwrite);
+
+ delete effVsXYrap;
+
+ }
+
+ file->Close();
+
+ // clean memory
+ delete effVsXY;
+
+}
+
+
+//---------------------------------------------------------------------------
+void PlotMuonEfficiency(TString fileNameData, TString fileNameSave, Bool_t saveEdges, Bool_t print, Bool_t draw)
+{
+ /// plot chamber, station and overall tracking efficiency
+
+ printf("plotting efficiency...\n");
+
+ // get input hists
+ TFile *file = new TFile(fileNameData.Data(), "read");
+ if (!file || !file->IsOpen()) {
+ printf("cannot open file %s \n",fileNameData.Data());
+ return;
+ }
+ TList *listTT = static_cast<TList*>(file->FindObjectAny("TotalTracksPerChamber"));
+ TList *listTD = static_cast<TList*>(file->FindObjectAny("TracksDetectedPerChamber"));
+ THnSparse *TT = static_cast<THnSparse*>(listTT->At(10));
+ THnSparse *TD = static_cast<THnSparse*>(listTD->At(10));
+
+ // output graphs
+ TGraphAsymmErrors *effVsCh[3] = {0x0, 0x0, 0x0};
+ TGraphAsymmErrors *effVsSt[3] = {0x0, 0x0, 0x0};
+ TString nameAdd[3] = {"", "Low", "Up"};
+ TString titleAdd[3] = {"", " - lower limit", " - upper limit"};
+ for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i) {
+ effVsCh[i] = CreateGraph(Form("chamberEff%s",nameAdd[i].Data()),
+ Form("Measured efficiency per chamber (0 = spectro)%s",titleAdd[i].Data()));
+ effVsSt[i] = CreateGraph(Form("stationEff%s",nameAdd[i].Data()),
+ Form("Measured efficiency per station (6 = st4&5)%s",titleAdd[i].Data()));
+ }
+
+ // set the centrality and pT ranges for integration
+ SetCentPtCh(*TT);
+ SetCentPtCh(*TD);
+
+ TArrayD chEff(11);
+ TArrayD chEffErr[2];
+ chEffErr[0].Set(11);
+ chEffErr[1].Set(11);
+
+ // compute chamber efficency and errors
+ if (GetChamberEfficiency(*TT, *TD, chEff, chEffErr, print)) {
+
+ // compute the overall tracking efficiency
+ GetTrackingEfficiency(chEff, chEffErr, effVsSt, effVsCh, 0, 0., print);
+
+ } else {
+
+ // set tracking efficiency to 1 +0 -1
+ effVsCh[0]->SetPoint(0,0.,1.);
+ effVsCh[0]->SetPointError(0,0.,0.,1.,0.);
+
+ if (saveEdges) {
+
+ // lower = 0 ± 0
+ effVsCh[1]->SetPoint(0,0.,0.);
+ effVsCh[1]->SetPointError(0,0.,0.,0.,0.);
+
+ // upper = 1 ± 0
+ effVsCh[2]->SetPoint(0,0.,1.);
+ effVsCh[2]->SetPointError(0,0.,0.,0.,0.);
+
+ }
+
+ }
+
+ // fill graph vs chamber
+ for (Int_t iCh = 1; iCh < 11; iCh++) {
+ for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i) {
+ effVsCh[i]->SetPoint(iCh,iCh,chEff[iCh]);
+ effVsCh[i]->SetPointError(iCh,0.,0.,chEffErr[0][iCh],chEffErr[1][iCh]);
+ }
+ }
+
+ // close input file
+ file->Close();
+
+ // display
+ for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i) {
+ effVsCh[i]->GetXaxis()->Set(12, -0.5, 10.5);
+ effVsCh[i]->GetXaxis()->SetNdivisions(11);
+ BeautifyGraph(*effVsCh[i], "chamber", "efficiency");
+ effVsSt[i]->GetXaxis()->Set(7, 0.5, 6.5);
+ effVsSt[i]->GetXaxis()->SetNdivisions(6);
+ BeautifyGraph(*effVsSt[i], "station", "efficiency");
+ }
+
+ if (draw) {
+ TCanvas *c = new TCanvas("cEfficiency", "Measured tracking efficiency" , 1000, 400);
+ c->Divide(2,1);
+ gROOT->SetSelectedPad(c->cd(1));
+ effVsCh[0]->DrawClone("ap");
+ gROOT->SetSelectedPad(c->cd(2));
+ effVsSt[0]->DrawClone("ap");
+ }
+
+ // save output
+ file = new TFile(fileNameSave.Data(),"update");
+ for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i) effVsCh[i]->Write(0x0, TObject::kOverwrite);
+ for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i) effVsSt[i]->Write(0x0, TObject::kOverwrite);
+ file->Close();
+
+ // clean memory
+ for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i) {
+ delete effVsCh[i];
+ delete effVsSt[i];
+ }
+
+}
+
+
+//---------------------------------------------------------------------------
+void PlotMuonEfficiencyVsRun(TString runList, TString fileNameData, TString fileNameSave, Bool_t print, Bool_t draw)
+{
+ /// plot chamber, station and overall tracking efficiency versus run
+
+ printf("plotting efficiency versus run...\n");
+
+ // open run list
+ ifstream inFile(runList.Data());
+ if (!inFile.is_open()) {
+ printf("cannot open file %s\n",runList.Data());
+ return;
+ }
+
+ // output graphs
+ TObjArray chamberVsRunGraphs; // 10 graphs: 1 for each chamber
+ chamberVsRunGraphs.SetOwner(kTRUE);
+ for ( Int_t iCh = 1; iCh < 11; ++iCh)
+ chamberVsRunGraphs.Add(CreateGraph("effCh%dVsRun", "Measured efficiency for chamber %d versus run",iCh));
+
+ TObjArray stationVsRunGraphs[3]; // 6 graphs: 1 for each station, and 1 for station 4&5
+ TGraphAsymmErrors *trkVsRun[3];
+ TString nameAdd[3] = {"", "Low", "Up"};
+ TString titleAdd[3] = {"", " - lower limit", " - upper limit"};
+ for (Int_t i = 0; i < 3; ++i) {
+
+ stationVsRunGraphs[i].SetOwner(kTRUE);
+ for ( Int_t iSt = 1; iSt < 6; ++iSt)
+ stationVsRunGraphs[i].Add(CreateGraph(Form("effSt%%dVsRun%s",nameAdd[i].Data()),
+ Form("Measured efficiency for station %%d versus run%s",titleAdd[i].Data()),iSt));
+ stationVsRunGraphs[i].Add(CreateGraph(Form("effSt4&5VsRun%s",nameAdd[i].Data()),
+ Form("Measured efficiency for station 4&5 versus run%s",titleAdd[i].Data())));
+
+ trkVsRun[i] = CreateGraph(Form("trackingEffVsRun%s",nameAdd[i].Data()),
+ Form("Measured tracking efficiency versus run%s", titleAdd[i].Data()));
+
+ }
+
+ Int_t irun = -1;
+ TList runs;
+ runs.SetOwner();
+
+ // loop over runs
+ while (!inFile.eof()) {
+
+ // get current run number
+ TString currRun;
+ currRun.ReadLine(inFile,kTRUE);
+ if(currRun.IsNull()) continue;
+ runs.AddLast(new TObjString(currRun));
+ irun++;
+
+ Int_t run = currRun.Atoi();
+
+ printf("run %d: ", run);
+
+ // compute efficiencies for this run
+ TString dataFile = Form("runs/%d/%s", run, fileNameData.Data());
+ TString outFile = Form("runs/%d/%s", run, fileNameSave.Data());
+ PlotMuonEfficiency(dataFile, outFile, kTRUE, print, kFALSE);
+
+ TFile *file = new TFile(outFile.Data(), "read");
+ if (!file || !file->IsOpen()) {
+ printf("cannot open file\n");
+ continue;
+ }
+
+ // loop over central value and edges
+ for (Int_t i = 0; i < 3; ++i) {
+
+ // get results
+ TGraphAsymmErrors *effVsCh;
+ TGraphAsymmErrors *effVsSt;
+ effVsCh = static_cast<TGraphAsymmErrors*>(file->FindObjectAny(Form("chamberEff%s",nameAdd[i].Data())));
+ effVsSt = static_cast<TGraphAsymmErrors*>(file->FindObjectAny(Form("stationEff%s",nameAdd[i].Data())));
+ if (!effVsCh || !effVsSt) {
+ printf("efficiency graph not found\n");
+ continue;
+ }
+
+ // fill chamber efficiency
+ if (i == 0) for ( Int_t iCh = 0; iCh < 10; ++iCh) {
+ TGraphAsymmErrors *g = static_cast<TGraphAsymmErrors*>(chamberVsRunGraphs.UncheckedAt(iCh));
+ g->SetPoint(irun,irun,effVsCh->GetY()[iCh+1]);
+ g->SetPointError(irun,0.,0.,effVsCh->GetErrorYlow(iCh+1),effVsCh->GetErrorYhigh(iCh+1));
+ }
+
+ // fill station efficiency
+ for ( Int_t iSt = 0; iSt < 6; ++iSt) {
+ TGraphAsymmErrors *g = static_cast<TGraphAsymmErrors*>(stationVsRunGraphs[i].UncheckedAt(iSt));
+ g->SetPoint(irun,irun,effVsSt->GetY()[iSt]);
+ g->SetPointError(irun,0.,0.,effVsSt->GetErrorYlow(iSt),effVsSt->GetErrorYhigh(iSt));
+ }
+
+ // fill spectrometer efficiency
+ trkVsRun[i]->SetPoint(irun,irun,effVsCh->GetY()[0]);
+ trkVsRun[i]->SetPointError(irun,0.,0.,effVsCh->GetErrorYlow(0),effVsCh->GetErrorYhigh(0));
+
+ }
+
+ file->Close();
+
+ }
+
+ inFile.close();
+
+ // display
+ BeautifyGraphs(chamberVsRunGraphs,"run number","efficiency");
+ SetRunLabel(chamberVsRunGraphs,irun,runs);
+ for (Int_t i = 0; i < 3; ++i) {
+ BeautifyGraph(*trkVsRun[i],"run number","efficiency");
+ SetRunLabel(*trkVsRun[i],irun,runs);
+ BeautifyGraphs(stationVsRunGraphs[i],"run number","efficiency");
+ SetRunLabel(stationVsRunGraphs[i],irun,runs);
+ }
+
+ if (draw) {
+ new TCanvas("cTrackingEffVsRun", "Tracking efficiency versus run",1000,400);
+ trkVsRun[0]->DrawClone("ap");
+ }
+
+ // save output
+ TFile* file = new TFile(fileNameSave.Data(),"update");
+ chamberVsRunGraphs.Write("ChamberEffVsRun", TObject::kOverwrite | TObject::kSingleKey);
+ for (Int_t i = 0; i < 3; ++i)
+ stationVsRunGraphs[i].Write(Form("StationEffVsRun%s",nameAdd[i].Data()), TObject::kOverwrite | TObject::kSingleKey);
+ for (Int_t i = 0; i < 3; ++i) trkVsRun[i]->Write(0x0, TObject::kOverwrite);
+ file->Close();
+
+ // clean memory
+ for (Int_t i = 0; i < 3; ++i) delete trkVsRun[i];
+
+}
+
+
+//---------------------------------------------------------------------------
+void PlotIntegratedMuonEfficiency(TString fileNameWeights, TString fileNameSave, Bool_t print, Bool_t draw)
+{
+ /// plot chamber, station and overall tracking efficiency integrated over runs
+
+ printf("plotting integrated efficiency...\n");
+
+ // load run weights
+ if (!fileNameWeights.IsNull()) LoadRunWeights(fileNameWeights);
+ if (!runWeights) {
+ printf("Cannot compute integrated efficiency without run-by-run weights\n");
+ return;
+ }
+
+ // get input hists
+ TFile *file = new TFile(fileNameSave.Data(), "update");
+ if (!file || !file->IsOpen()) {
+ printf("cannot open file\n");
+ return;
+ }
+ TObjArray *chamberVsRunGraphs = static_cast<TObjArray*>(file->FindObjectAny("ChamberEffVsRun"));
+ TObjArray *stationVsRunGraphs[2];
+ stationVsRunGraphs[0] = static_cast<TObjArray*>(file->FindObjectAny("StationEffVsRunLow"));
+ stationVsRunGraphs[1] = static_cast<TObjArray*>(file->FindObjectAny("StationEffVsRunUp"));
+ TGraphAsymmErrors *trkVsRun[2];
+ trkVsRun[0] = static_cast<TGraphAsymmErrors*>(file->FindObjectAny("trackingEffVsRunLow"));
+ trkVsRun[1] = static_cast<TGraphAsymmErrors*>(file->FindObjectAny("trackingEffVsRunUp"));
+ if (!chamberVsRunGraphs || !stationVsRunGraphs[0] || !stationVsRunGraphs[1] || !trkVsRun[0] || !trkVsRun[1]) {
+ printf("object not found --> you must first plot the efficiency versus run\n");
+ return;
+ }
+
+ // output graphs
+ TGraphAsymmErrors *effVsCh = CreateGraph("integratedChamberEff", "Integrated efficiency per chamber (0 = spectro)");
+ TGraphAsymmErrors *effVsSt = CreateGraph("integratedStationEff", "Integrated efficiency per station (6 = st4&5)");
+
+ // integrate spectrometer efficiency
+ IntegrateMuonEfficiency(*trkVsRun[0], *trkVsRun[1], *effVsCh, 0, 0.);
+
+ // integrate chamber efficiency
+ for ( Int_t iCh = 0; iCh < 10; ++iCh) {
+ TGraphAsymmErrors *g = static_cast<TGraphAsymmErrors*>(chamberVsRunGraphs->UncheckedAt(iCh));
+ IntegrateMuonEfficiency(*g, *g, *effVsCh, iCh+1, iCh+1);
+ }
+
+ // integrate station efficiency
+ for ( Int_t iSt = 0; iSt < 6; ++iSt) {
+ TGraphAsymmErrors *gLow = static_cast<TGraphAsymmErrors*>(stationVsRunGraphs[0]->UncheckedAt(iSt));
+ TGraphAsymmErrors *gUp = static_cast<TGraphAsymmErrors*>(stationVsRunGraphs[1]->UncheckedAt(iSt));
+ IntegrateMuonEfficiency(*gLow, *gUp, *effVsSt, iSt, iSt+1);
+ }
+
+ // print results
+ if (print) {
+ for (Int_t iCh = 1; iCh < 11; ++iCh) {
+ cout << "Efficiency chamber " << iCh << " : ";
+ cout << effVsCh->GetY()[iCh] << " + " << effVsCh->GetErrorYhigh(iCh) << " - " << effVsCh->GetErrorYlow(iCh) << endl;
+ }
+ for (Int_t iSt = 0; iSt < 6; ++iSt) {
+ if (iSt < 5) cout << "Station " << iSt+1 << " = ";
+ else cout << "Station 45 = ";
+ cout << effVsSt->GetY()[iSt] << " + " << effVsSt->GetErrorYhigh(iSt) << " - " << effVsSt->GetErrorYlow(iSt) << endl;
+ }
+ cout << "Total tracking efficiency : ";
+ cout << effVsCh->GetY()[0] << " + " << effVsCh->GetErrorYhigh(0) << " - " << effVsCh->GetErrorYlow(0) << endl << endl;
+ }
+
+ // display
+ effVsCh->GetXaxis()->Set(12, -0.5, 10.5);
+ effVsCh->GetXaxis()->SetNdivisions(11);
+ BeautifyGraph(*effVsCh, "chamber", "efficiency");
+ effVsSt->GetXaxis()->Set(7, 0.5, 6.5);
+ effVsSt->GetXaxis()->SetNdivisions(6);
+ BeautifyGraph(*effVsSt, "station", "efficiency");
+
+ if (draw) {
+ TCanvas *c = new TCanvas("cIntegratedEfficiency", "Integrated tracking efficiency" , 1000, 400);
+ c->Divide(2,1);
+ gROOT->SetSelectedPad(c->cd(1));
+ effVsCh->DrawClone("ap");
+ gROOT->SetSelectedPad(c->cd(2));
+ effVsSt->DrawClone("ap");
+ }
+
+ // save output
+ effVsCh->Write(0x0, TObject::kOverwrite);
+ effVsSt->Write(0x0, TObject::kOverwrite);
+ file->Close();
+
+ // clean memory
+ delete effVsCh;
+ delete effVsSt;
+
+}
+
+
+//---------------------------------------------------------------------------
+void PlotMuonEfficiencyPerDE(TString fileNameData, TString fileNameSave, Bool_t saveEdges)
+{
+ /// plot chamber and station efficiency per DE
+
+ printf("plotting efficiency per DE...\n");
+
+ // get input hists
+ TFile *file = new TFile(fileNameData.Data(), "read");
+ if (!file || !file->IsOpen()) {
+ printf("cannot open file %s \n",fileNameData.Data());
+ return;
+ }
+ TList *listTT = static_cast<TList*>(file->FindObjectAny("TotalTracksPerChamber"));
+ TList *listTD = static_cast<TList*>(file->FindObjectAny("TracksDetectedPerChamber"));
+
+ // output graph arrays
+ TObjArray chamberVsDEGraphs; // 10 graphs: 1 for each chamber
+ chamberVsDEGraphs.SetOwner(kTRUE);
+ TGraphAsymmErrors *gCh[10]; // shortcut
+
+ TObjArray stationVsDEGraphs[3]; // 6 graphs: 1 for each station, and 1 for station 4&5
+ TString nameAdd[3] = {"", "Low", "Up"};
+ TString titleAdd[3] = {"", " - lower limit", " - upper limit"};
+ for (Int_t i = 0; i < 3; ++i) stationVsDEGraphs[i].SetOwner(kTRUE);
+ TGraphAsymmErrors *gSt[6][3]; // shortcut
+ for (Int_t iSt = 0; iSt < 6; ++iSt) for (Int_t i = 0; i < 3; ++i) gSt[iSt][i] = 0x0;
+
+ // loop over chambers
+ for (Int_t iCh = 0; iCh < 10; ++iCh) {
+
+ // output graphs
+ chamberVsDEGraphs.Add(CreateGraph("effCh%dVsDE","Measured efficiency for chamber %d per DE",iCh+1));
+ gCh[iCh] = static_cast<TGraphAsymmErrors*>(chamberVsDEGraphs.UncheckedAt(iCh));
+
+ // get input hists
+ THnSparse *TT = static_cast<THnSparse*>(listTT->At(iCh));
+ THnSparse *TD = static_cast<THnSparse*>(listTD->At(iCh));
+
+ // set the centrality and pT range for integration
+ SetCentPtCh(*TT);
+ SetCentPtCh(*TD);
+
+ // compute DE efficency and errors
+ GetDEEfficiency(*TT, *TD, *gCh[iCh]);
+
+ }
+
+ // close input file
+ file->Close();
+
+ TArrayD chEff(11);
+ TArrayD chEffErr[2];
+ chEffErr[0].Set(11);
+ chEffErr[1].Set(11);
+
+ // loop over the first 3 stations
+ for (Int_t iSt = 0; iSt < 3; ++iSt) {
+
+ // output graphs
+ for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i) {
+ stationVsDEGraphs[i].Add(CreateGraph(Form("effSt%%dVsDE%s",nameAdd[i].Data()),
+ Form("Measured efficiency for station %%d per DE%s",titleAdd[i].Data()),iSt+1));
+ gSt[iSt][i] = static_cast<TGraphAsymmErrors*>(stationVsDEGraphs[i].UncheckedAt(iSt));
+ }
+
+ // loop over DE
+ Int_t nDE = gCh[2*iSt]->GetN();
+ for (Int_t iDE = 0; iDE < nDE; ++iDE) {
+
+ // copy efficiency
+ for (Int_t iCh = 0; iCh < 2; ++iCh) {
+ chEff[2*iSt+iCh+1] = gCh[2*iSt+iCh]->GetY()[iDE];
+ chEffErr[0][2*iSt+iCh+1] = gCh[2*iSt+iCh]->GetErrorYlow(iDE);
+ chEffErr[1][2*iSt+iCh+1] = gCh[2*iSt+iCh]->GetErrorYhigh(iDE);
+ }
+
+ // compute station efficiency
+ GetStationEfficiency(chEff, chEffErr, iSt, gSt[iSt], iDE, iDE);
+
+ }
+
+ }
+
+ // output graphs for last 2 stations
+ for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i) {
+ for (Int_t iSt = 3; iSt < 5; ++iSt) {
+ stationVsDEGraphs[i].Add(CreateGraph(Form("effSt%%dVsDE%s",nameAdd[i].Data()),
+ Form("Measured efficiency for station %%d per DE%s",titleAdd[i].Data()),iSt+1));
+ gSt[iSt][i] = static_cast<TGraphAsymmErrors*>(stationVsDEGraphs[i].UncheckedAt(iSt));
+ }
+ stationVsDEGraphs[i].Add(CreateGraph(Form("effSt4&5VsDE%s",nameAdd[i].Data()),
+ Form("Measured efficiency for station 4&5 per DE%s",titleAdd[i].Data())));
+ gSt[5][i] = static_cast<TGraphAsymmErrors*>(stationVsDEGraphs[i].UncheckedAt(5));
+ }
+
+ // loop over DE
+ Int_t nDE = gCh[6]->GetN();
+ for (Int_t iDE = 0; iDE < nDE; ++iDE) {
+
+ // copy efficiency
+ for (Int_t iCh = 6; iCh < 10; ++iCh) {
+ chEff[iCh+1] = gCh[iCh]->GetY()[iDE];
+ chEffErr[0][iCh+1] = gCh[iCh]->GetErrorYlow(iDE);
+ chEffErr[1][iCh+1] = gCh[iCh]->GetErrorYhigh(iDE);
+ }
+
+ // compute station 4&5 efficiency individually
+ for (Int_t iSt = 3; iSt < 5; ++iSt) GetStationEfficiency(chEff, chEffErr, iSt, gSt[iSt], iDE, iDE);
+
+ // compute station 4&5 efficiency together
+ GetStation45Efficiency(chEff, chEffErr, gSt[5], iDE, iDE);
+
+ }
+
+ // display
+ for (Int_t iCh = 0; iCh < 10; ++iCh) {
+ Int_t nDE = gCh[iCh]->GetN();
+ gCh[iCh]->GetXaxis()->Set(nDE+1, -0.5, nDE-0.5);
+ gCh[iCh]->GetXaxis()->SetNdivisions(nDE);
+ }
+ BeautifyGraphs(chamberVsDEGraphs,"Detection Element","efficiency");
+ for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i) {
+ for (Int_t iSt = 0; iSt < 6; ++iSt) {
+ Int_t nDE = gSt[iSt][i]->GetN();
+ gSt[iSt][i]->GetXaxis()->Set(nDE+1, -0.5, nDE-0.5);
+ gSt[iSt][i]->GetXaxis()->SetNdivisions(nDE);
+ }
+ BeautifyGraphs(stationVsDEGraphs[i],"Detection Element","efficiency");
+
+ }
+
+ // save Output
+ file = new TFile(fileNameSave.Data(),"update");
+ chamberVsDEGraphs.Write("ChamberEffVsDE", TObject::kOverwrite | TObject::kSingleKey);
+ for (Int_t i = 0; i < 1 || (saveEdges && i < 3); ++i)
+ stationVsDEGraphs[i].Write(Form("StationEffVsDE%s",nameAdd[i].Data()), TObject::kOverwrite | TObject::kSingleKey);
+ file->Close();
+
+}
+
+
+//---------------------------------------------------------------------------
+void PlotMuonEfficiencyPerDEVsRun(TString runList, TString fileNameData, TString fileNameSave)
+{
+ /// plot chamber and station efficiency per DE versus run
+
+ printf("plotting efficiency per DE versus run...\n");
+
+ // open run list
+ ifstream inFile(runList.Data());
+ if (!inFile.is_open()) {
+ printf("cannot open file %s\n",runList.Data());
+ return;
+ }
+
+ // output graphs
+ TObjArray deVsRunGraphs; // 1 graph per DE
+ deVsRunGraphs.SetOwner(kTRUE);
+ TObjArray stDEVsRunGraphs[3]; // 1 graph per pair (quartet) of DE in individual stations (stations 4&5)
+ TString nameAdd[3] = {"", "Low", "Up"};
+ TString titleAdd[3] = {"", " - lower limit", " - upper limit"};
+ for (Int_t i = 0; i < 3; ++i) stDEVsRunGraphs[i].SetOwner(kTRUE);
+ Bool_t createGraph = kTRUE;
+
+ Int_t irun = -1;
+ TList runs;
+ runs.SetOwner();
+
+ // loop over runs
+ while (!inFile.eof()) {
+
+ // get current run number
+ TString currRun;
+ currRun.ReadLine(inFile,kTRUE);
+ if(currRun.IsNull()) continue;
+ runs.AddLast(new TObjString(currRun));
+ irun++;
+
+ Int_t run = currRun.Atoi();
+
+ printf("run %d: ", run);
+
+ // compute efficiencies for this run
+ TString dataFile = Form("runs/%d/%s", run, fileNameData.Data());
+ TString outFile = Form("runs/%d/%s", run, fileNameSave.Data());
+ PlotMuonEfficiencyPerDE(dataFile, outFile, kTRUE);
+
+ TFile *file = new TFile(outFile.Data(), "read");
+ if (!file || !file->IsOpen()) {
+ printf("cannot open file\n");
+ continue;
+ }
+
+ // get results
+ TObjArray *chamberVsDEGraphs = static_cast<TObjArray*>(file->FindObjectAny("ChamberEffVsDE"));
+ if (!chamberVsDEGraphs) {
+ printf("efficiency graph not found\n");
+ continue;
+ }
+
+ Int_t currentDE = 0;
+
+ // loop over chambers
+ for ( Int_t iCh = 0; iCh < 10; ++iCh) {
+
+ // input graph
+ TGraphAsymmErrors *g = static_cast<TGraphAsymmErrors*>(chamberVsDEGraphs->UncheckedAt(iCh));
+ Int_t nDE = g->GetN();
+
+ // loop over DE
+ for (Int_t iDE = 0; iDE < nDE; ++iDE) {
+
+ // output graph
+ if (createGraph) deVsRunGraphs.Add(CreateGraph("effDE%dVsRun","Measured efficiency for DE %d versus run",100*(iCh+1)+iDE));
+ TGraphAsymmErrors *gDE= static_cast<TGraphAsymmErrors*>(deVsRunGraphs.UncheckedAt(currentDE++));
+
+ // fill DE efficiency
+ gDE->SetPoint(irun,irun,g->GetY()[iDE]);
+ gDE->SetPointError(irun,0.,0.,g->GetErrorYlow(iDE),g->GetErrorYhigh(iDE));
+
+ }
+
+ }
+
+ // loop over central value and edges
+ for (Int_t i = 0; i < 3; ++i) {
+
+ // get results
+ TObjArray *stationVsDEGraphs = static_cast<TObjArray*>(file->FindObjectAny(Form("StationEffVsDE%s",nameAdd[i].Data())));
+ if (!stationVsDEGraphs) {
+ printf("efficiency graph not found\n");
+ continue;
+ }
+
+ Int_t currentStDE = 0;
+
+ // loop over stations
+ for ( Int_t iSt = 0; iSt < 6; ++iSt) {
+
+ // input graph
+ TGraphAsymmErrors *g = static_cast<TGraphAsymmErrors*>(stationVsDEGraphs->UncheckedAt(iSt));
+ Int_t nDE = g->GetN();
+
+ // loop over DE
+ for (Int_t iDE = 0; iDE < nDE; ++iDE) {
+
+ // output graph
+ if (createGraph) {
+ TString sSt = (iSt<5) ? Form("%d",iSt+1) : "4&5";
+ stDEVsRunGraphs[i].Add(CreateGraph(Form("effSt%sDE%%dVsRun%s",sSt.Data(),nameAdd[i].Data()),
+ Form("Measured efficiency for DE %%d in station %s versus run%s",sSt.Data(),titleAdd[i].Data()),iDE));
+ }
+ TGraphAsymmErrors *gDE= static_cast<TGraphAsymmErrors*>(stDEVsRunGraphs[i].UncheckedAt(currentStDE++));
+
+ // fill DE efficiency
+ gDE->SetPoint(irun,irun,g->GetY()[iDE]);
+ gDE->SetPointError(irun,0.,0.,g->GetErrorYlow(iDE),g->GetErrorYhigh(iDE));
+
+ }
+
+ }
+
+ }
+
+ file->Close();
+
+ createGraph = kFALSE;
+
+ }
+
+ inFile.close();
+
+ // display
+ BeautifyGraphs(deVsRunGraphs,"run number","efficiency");
+ SetRunLabel(deVsRunGraphs,irun,runs);
+ for (Int_t i = 0; i < 3; ++i) {
+ BeautifyGraphs(stDEVsRunGraphs[i],"run number","efficiency");
+ SetRunLabel(stDEVsRunGraphs[i],irun,runs);
+ }
+
+ // save output
+ TFile* file = new TFile(fileNameSave.Data(),"update");
+ deVsRunGraphs.Write("DEEffVsRun", TObject::kOverwrite | TObject::kSingleKey);
+ for (Int_t i = 0; i < 3; ++i)
+ stDEVsRunGraphs[i].Write(Form("DEEffPerStationVsRun%s",nameAdd[i].Data()), TObject::kOverwrite | TObject::kSingleKey);
+ file->Close();
+
+}
+
+
+//---------------------------------------------------------------------------
+void PlotIntegratedMuonEfficiencyPerDE(TString fileNameWeights, TString fileNameSave)
+{
+ /// plot chamber and station efficiency per DE integrated over runs
+
+ printf("plotting integrated efficiency per DE...\n");
+
+ // load run weights
+ if (!fileNameWeights.IsNull()) LoadRunWeights(fileNameWeights);
+ if (!runWeights) {
+ printf("Cannot compute integrated efficiency without run-by-run weights\n");
+ return;
+ }
+
+ // get input hists
+ TFile *file = new TFile(fileNameSave.Data(), "update");
+ if (!file || !file->IsOpen()) {
+ printf("cannot open file\n");
+ return;
+ }
+ TObjArray *deVsRunGraphs = static_cast<TObjArray*>(file->FindObjectAny("DEEffVsRun"));
+ TObjArray *stDEVsRunGraphs[2];
+ stDEVsRunGraphs[0] = static_cast<TObjArray*>(file->FindObjectAny("DEEffPerStationVsRunLow"));
+ stDEVsRunGraphs[1] = static_cast<TObjArray*>(file->FindObjectAny("DEEffPerStationVsRunUp"));
+ if (!deVsRunGraphs || !stDEVsRunGraphs[0] || !stDEVsRunGraphs[1]) {
+ printf("object not found --> you must first plot the efficiency versus run\n");
+ return;
+ }
+
+ // output graph arrays
+ TObjArray chamberVsDEGraphs; // 10 graphs: 1 for each chamber
+ chamberVsDEGraphs.SetOwner(kTRUE);
+ for ( Int_t iCh = 0; iCh < 10; ++iCh)
+ chamberVsDEGraphs.Add(CreateGraph("integratedEffCh%dVsDE","Integrated efficiency for chamber %d per DE",iCh+1));
+
+ TObjArray stationVsDEGraphs; // 6 graphs: 1 for each station, and 1 for station 4&5
+ stationVsDEGraphs.SetOwner(kTRUE);
+ for ( Int_t iSt = 0; iSt < 5; ++iSt)
+ stationVsDEGraphs.Add(CreateGraph("integratedEffSt%dVsDE","Integrated efficiency for station %d per DE",iSt+1));
+ stationVsDEGraphs.Add(CreateGraph("integratedEffSt4&5VsDE","Integrated efficiency for station 4&5 per DE"));
+
+ // Loop over DE
+ TIter nextDE(deVsRunGraphs);
+ TGraphAsymmErrors *gDE = 0x0;
+ while ((gDE = static_cast<TGraphAsymmErrors*>(nextDE()))) {
+
+ // get chamber and DE indices
+ Int_t deId;
+ sscanf(gDE->GetName(), "effDE%dVsRun", &deId);
+ Int_t iCh = deId/100-1;
+ Int_t iDE = deId%100;
+
+ // integrate DE efficiency
+ TGraphAsymmErrors *g = static_cast<TGraphAsymmErrors*>(chamberVsDEGraphs.UncheckedAt(iCh));
+ IntegrateMuonEfficiency(*gDE, *gDE, *g, iDE, iDE);
+
+ }
+
+ // Loop over DE per station
+ Int_t ng = stDEVsRunGraphs[0]->GetEntries();
+ for (Int_t ig = 0; ig < ng; ++ig) {
+
+ TGraphAsymmErrors *gDELow = static_cast<TGraphAsymmErrors*>(stDEVsRunGraphs[0]->UncheckedAt(ig));
+ TGraphAsymmErrors *gDEUp = static_cast<TGraphAsymmErrors*>(stDEVsRunGraphs[1]->UncheckedAt(ig));
+
+ // get station and DE indices
+ Int_t iSt, iDE;
+ if (strstr(gDELow->GetName(), "4&5")) {
+ iSt = 5;
+ sscanf(gDELow->GetName(), "effSt4&5DE%dVsRunLow", &iDE);
+ } else {
+ sscanf(gDELow->GetName(), "effSt%dDE%dVsRunLow", &iSt, &iDE);
+ iSt--;
+ }
+
+ // Integrate DE efficiency per station
+ TGraphAsymmErrors *g = static_cast<TGraphAsymmErrors*>(stationVsDEGraphs.UncheckedAt(iSt));
+ IntegrateMuonEfficiency(*gDELow, *gDEUp, *g, iDE, iDE);
+
+ }
+
+ // display
+ for ( Int_t iCh = 0; iCh < 10; ++iCh) {
+ TGraphAsymmErrors *g = static_cast<TGraphAsymmErrors*>(chamberVsDEGraphs.UncheckedAt(iCh));
+ Int_t nDE = g->GetN();
+ g->GetXaxis()->Set(nDE+1, -0.5, nDE-0.5);
+ g->GetXaxis()->SetNdivisions(nDE);
+ }
+ BeautifyGraphs(chamberVsDEGraphs,"Detection Element","efficiency");
+ for ( Int_t iSt = 0; iSt < 6; ++iSt) {
+ TGraphAsymmErrors *g = static_cast<TGraphAsymmErrors*>(stationVsDEGraphs.UncheckedAt(iSt));
+ Int_t nDE = g->GetN();
+ g->GetXaxis()->Set(nDE+1, -0.5, nDE-0.5);
+ g->GetXaxis()->SetNdivisions(nDE);
+ }
+ BeautifyGraphs(stationVsDEGraphs,"Detection Element","efficiency");
-/*
+ // save Output
+ chamberVsDEGraphs.Write("IntegratedChamberEffVsDE", TObject::kOverwrite | TObject::kSingleKey);
+ stationVsDEGraphs.Write("IntegratedStationEffVsDE", TObject::kOverwrite | TObject::kSingleKey);
+ file->Close();
+
+}
-This macro is meant to be used with the output of the Efficiency Task. It computes the uncorrelated efficiencies per chammber, station and globally. It then looks for correlated dead arrays and corrects the previous calculated efficiency.
-Macro created by Lizardo Valencia Palomo
+//---------------------------------------------------------------------------
+Bool_t GetChamberEfficiency(THnSparse &TT, THnSparse &TD, TArrayD &chEff, TArrayD chEffErr[2], Bool_t printError)
+{
+ /// compute chamber efficiency and errors
+ /// return kFALSE if efficiency unknown for all chambers
+
+ // project track hists over the chamber axis
+ TH1D *TTdraw = TT.Projection(0,"e");
+ TH1D *TDdraw = TD.Projection(0,"e");
+
+ // compute chamber efficiency and errors
+ TGraphAsymmErrors *efficiency = (TTdraw->GetEntries() > 0) ? new TGraphAsymmErrors(TDdraw, TTdraw,Form("%se0",effErrMode)) : 0x0;
+ Bool_t ok = (efficiency);
+
+ // fill arrays
+ if (ok) {
+
+ Bool_t missingEff = kFALSE;
+
+ for (Int_t i = 0; i < 10; i++) {
+
+ if (TTdraw->GetBinContent(i+1) > 0) {
+
+ chEff[i+1] = efficiency->GetY()[i];
+ chEffErr[0][i+1] = efficiency->GetErrorYlow(i);
+ chEffErr[1][i+1] = efficiency->GetErrorYhigh(i);
+
+ } else {
+
+ chEff[i+1] = -1.;
+ chEffErr[0][i+1] = 0.;
+ chEffErr[1][i+1] = 0.;
+
+ missingEff = kTRUE;
+
+ }
+
+ }
+
+ if (missingEff && printError) cout << "efficiency partially unknown" << endl;
+
+ } else {
+
+ for (Int_t i = 0; i < 10; i++) {
+
+ chEff[i+1] = -1.;
+ chEffErr[0][i+1] = 0.;
+ chEffErr[1][i+1] = 0.;
+
+ }
+
+ if (printError) cout << "efficiency unknown" << endl << endl;
+
+ }
+
+ // clean memory
+ delete TTdraw;
+ delete TDdraw;
+ delete efficiency;
+
+ return ok;
+
+}
-lizardo.valencia.palomo@cern.ch
-*/
+//---------------------------------------------------------------------------
+void GetDEEfficiency(THnSparse &TT, THnSparse &TD, TGraphAsymmErrors &effVsDE)
+{
+ /// compute DE efficiency and errors
+
+ // project track hists over the chamber axis
+ TH1D *TTdraw = TT.Projection(0,"e");
+ TH1D *TDdraw = TD.Projection(0,"e");
+
+ // compute DE efficiency and errors
+ TGraphAsymmErrors *efficiency = (TTdraw->GetEntries() > 0) ? new TGraphAsymmErrors(TDdraw, TTdraw,Form("%se0",effErrMode)) : 0x0;
+ Int_t nDE = TTdraw->GetNbinsX();
+
+ if (efficiency) {
+
+ for (Int_t iDE = 0; iDE < nDE; ++iDE) {
+
+ if (TTdraw->GetBinContent(iDE+1) > 0) {
+
+ effVsDE.SetPoint(iDE,iDE,efficiency->GetY()[iDE]);
+ effVsDE.SetPointError(iDE,0,0,efficiency->GetErrorYlow(iDE),efficiency->GetErrorYhigh(iDE));
+
+ } else {
+
+ effVsDE.SetPoint(iDE,iDE,-1);
+ effVsDE.SetPointError(iDE,0,0,0,0);
+
+ }
+
+ }
+
+ } else {
+
+ for (Int_t iDE = 0; iDE < nDE; ++iDE) {
+
+ effVsDE.SetPoint(iDE,iDE,-1);
+ effVsDE.SetPointError(iDE,0,0,0,0);
+
+ }
+
+ }
+
+ // clean memory
+ delete TTdraw;
+ delete TDdraw;
+ delete efficiency;
+
+}
+//---------------------------------------------------------------------------
+void ComputeStationEfficiency(TArrayD &chEff, TArrayD chEffErr[2], Int_t iSt, Double_t &stEff, Double_t stEffErr[2])
+{
+ /// compute the station iSt (0...4) efficiency and errors from the individual chamber efficiencies and errors
+
+ stEff = 1.-(1.-chEff[2*iSt+1])*(1.-chEff[2*iSt+2]);
+
+ Double_t d1 = (1. - chEff[2*iSt+2]); d1 *= d1;
+ Double_t d2 = (1. - chEff[2*iSt+1]); d2 *= d2;
+
+ for (Int_t i = 0; i < 2; ++i) {
+ Double_t s1 = chEffErr[i][2*iSt+1] * chEffErr[i][2*iSt+1];
+ Double_t s2 = chEffErr[i][2*iSt+2] * chEffErr[i][2*iSt+2];
+ stEffErr[i] = TMath::Sqrt(d1*s1 + d2*s2 + s1*s2);
+ }
+
+ stEffErr[0] = TMath::Min(stEff, stEffErr[0]);
+ stEffErr[1] = TMath::Min(1.-stEff, stEffErr[1]);
+
+}
+
+
+//---------------------------------------------------------------------------
+void GetStationEfficiency(TArrayD &chEff, TArrayD chEffErr[2], Int_t iSt, TGraphAsymmErrors *effVsX[3], Int_t ip, Double_t xp)
+{
+ /// compute the station iSt (0...4) efficiency and errors from the individual chamber efficiencies and errors
+
+ if (chEff[2*iSt+1] >= 0 && chEff[2*iSt+2] >= 0) {
+
+ // compute station efficiency from known chamber efficiency
+ Double_t stEff, stEffErr[2];
+ ComputeStationEfficiency(chEff, chEffErr, iSt, stEff, stEffErr);
+
+ // fill graphs if required
+ for (Int_t i = 0; i < 3; ++i) {
+ if (effVsX[i]) {
+ effVsX[i]->SetPoint(ip,xp,stEff);
+ effVsX[i]->SetPointError(ip,0.,0.,stEffErr[0],stEffErr[1]);
+ }
+ }
+
+ } else {
+
+ Double_t edge[2] = {0., 1.};
+ TArrayD chEffEdge[2];
+ Double_t stEffEdge[2];
+ Double_t stEffEdgeErr[2][2];
+
+ for (Int_t i = 0; i < 2; ++i) {
+
+ // set lower(upper) limit of chamber efficiency
+ chEffEdge[i].Set(11);
+ for (Int_t iCh = 1; iCh < 3; ++iCh) chEffEdge[i][2*iSt+iCh] = (chEff[2*iSt+iCh] < 0) ? edge[i] : chEff[2*iSt+iCh];
+
+ // compute station efficiency
+ ComputeStationEfficiency(chEffEdge[i], chEffErr, iSt, stEffEdge[i], stEffEdgeErr[i]);
+
+ // fill graph if required
+ if (effVsX[i+1]) {
+ effVsX[i+1]->SetPoint(ip,xp,stEffEdge[i]);
+ effVsX[i+1]->SetPointError(ip,0.,0.,stEffEdgeErr[i][0],stEffEdgeErr[i][1]);
+ }
+
+ }
+
+ // combine extreme cases to get station efficiency and fill graph if required
+ if (effVsX[0]) {
+ effVsX[0]->SetPoint(ip,xp,stEffEdge[1]);
+ effVsX[0]->SetPointError(ip,0.,0.,stEffEdge[1]-stEffEdge[0]+stEffEdgeErr[0][0],stEffEdgeErr[1][1]);
+ }
+
+ }
+
+}
+
+
+//---------------------------------------------------------------------------
+void ComputeStation45Efficiency(TArrayD &chEff, TArrayD chEffErr[2], Double_t &st45Eff, Double_t st45EffErr[2])
+{
+ /// compute the station 4-5 efficiency and errors from chamber efficiencies and errors
+
+ 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];
+
+ Double_t d1 = chEff[8]*chEff[9] + chEff[8]*chEff[10] + chEff[9]*chEff[10] - 3.*chEff[8]*chEff[9]*chEff[10]; d1 *= d1;
+ Double_t d2 = chEff[7]*chEff[9] + chEff[7]*chEff[10] + chEff[9]*chEff[10] - 3.*chEff[7]*chEff[9]*chEff[10]; d2 *= d2;
+ Double_t d3 = chEff[7]*chEff[8] + chEff[7]*chEff[10] + chEff[8]*chEff[10] - 3.*chEff[7]*chEff[8]*chEff[10]; d3 *= d3;
+ Double_t d4 = chEff[7]*chEff[8] + chEff[7]*chEff[9] + chEff[8]*chEff[9] - 3.*chEff[7]*chEff[8]*chEff[9]; d4 *= d4;
+ Double_t d12 = chEff[9] + chEff[10] - 3.*chEff[9]*chEff[10]; d12 *= d12;
+ Double_t d13 = chEff[8] + chEff[10] - 3.*chEff[8]*chEff[10]; d13 *= d13;
+ Double_t d14 = chEff[8] + chEff[9] - 3.*chEff[8]*chEff[9]; d14 *= d14;
+ Double_t d23 = chEff[7] + chEff[10] - 3.*chEff[7]*chEff[10]; d23 *= d23;
+ Double_t d24 = chEff[7] + chEff[9] - 3.*chEff[7]*chEff[9]; d24 *= d24;
+ Double_t d34 = chEff[7] + chEff[8] - 3.*chEff[7]*chEff[8]; d34 *= d34;
+ Double_t d123 = 1. - 3.*chEff[10]; d123 *= d123;
+ Double_t d124 = 1. - 3.*chEff[9]; d124 *= d124;
+ Double_t d134 = 1. - 3.*chEff[8]; d134 *= d134;
+ Double_t d234 = 1. - 3.*chEff[7]; d234 *= d234;
+ Double_t d1234 = - 3.; d1234 *= d1234;
+
+ for (Int_t i = 0; i < 2; ++i) {
+ Double_t s1 = chEffErr[i][7] * chEffErr[i][7];
+ Double_t s2 = chEffErr[i][8] * chEffErr[i][8];
+ Double_t s3 = chEffErr[i][9] * chEffErr[i][9];
+ Double_t s4 = chEffErr[i][10] * chEffErr[i][10];
+ 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);
+ }
+
+ st45EffErr[0] = TMath::Min(st45Eff, st45EffErr[0]);
+ st45EffErr[1] = TMath::Min(1.-st45Eff, st45EffErr[1]);
+
+}
+
-TH1D* GraphToHist(Bool_t ToBePlotted, const TGraphAsymmErrors* Graph, Int_t Chamber);
-std::vector<Double_t> FillVector(TH1D* Chamber, Bool_t Errors);
-std::vector<Double_t> CorrelatedDeadArray(std::vector<Double_t> Chj, std::vector<Double_t> Chi, Int_t NumberOfDEs, Int_t SymmetricGroups);
-void Drawing(std::vector<Double_t> Chi, std::vector<Double_t> Chj, std::vector<Double_t> ChiErr, std::vector<Double_t> ChjErr, Int_t Bins, Int_t SymGroups, Int_t Station);
-Double_t N00(Double_t Average, Double_t CDA, Int_t SimGroup, Int_t SimGroupsPerCh);
+//---------------------------------------------------------------------------
+void GetStation45Efficiency(TArrayD &chEff, TArrayD chEffErr[2], TGraphAsymmErrors *effVsX[3], Int_t ip, Double_t xp)
+{
+ /// compute the station 4-5 efficiency and errors from chamber efficiencies and errors
+
+ if (chEff[7] >= 0 && chEff[8] >= 0 && chEff[9] >= 0 && chEff[10] >= 0) {
+
+ // compute station efficiency from known chamber efficiency
+ Double_t stEff, stEffErr[2];
+ ComputeStation45Efficiency(chEff, chEffErr, stEff, stEffErr);
+
+ // fill graphs if required
+ for (Int_t i = 0; i < 3; ++i) {
+ if (effVsX[i]) {
+ effVsX[i]->SetPoint(ip,xp,stEff);
+ effVsX[i]->SetPointError(ip,0.,0.,stEffErr[0],stEffErr[1]);
+ }
+ }
+
+ } else {
+
+ Double_t edge[2] = {0., 1.};
+ TArrayD chEffEdge[2];
+ Double_t stEffEdge[2];
+ Double_t stEffEdgeErr[2][2];
+
+ for (Int_t i = 0; i < 2; ++i) {
+
+ // set lower(upper) limit of chamber efficiency
+ chEffEdge[i].Set(11);
+ for (Int_t iCh = 7; iCh < 11; ++iCh) chEffEdge[i][iCh] = (chEff[iCh] < 0) ? edge[i] : chEff[iCh];
+
+ // compute station efficiency
+ ComputeStation45Efficiency(chEffEdge[i], chEffErr, stEffEdge[i], stEffEdgeErr[i]);
+
+ // fill graph if required
+ if (effVsX[i+1]) {
+ effVsX[i+1]->SetPoint(ip,xp,stEffEdge[i]);
+ effVsX[i+1]->SetPointError(ip,0.,0.,stEffEdgeErr[i][0],stEffEdgeErr[i][1]);
+ }
+
+ }
+
+ // combine extreme cases to get station efficiency and fill graph if required
+ if (effVsX[0]) {
+ effVsX[0]->SetPoint(ip,xp,stEffEdge[1]);
+ effVsX[0]->SetPointError(ip,0.,0.,stEffEdge[1]-stEffEdge[0]+stEffEdgeErr[0][0],stEffEdgeErr[1][1]);
+ }
+
+ }
+
+}
-void MuonTrackingEfficiency(TString fileName = "./AnalysisResults.root")
+//---------------------------------------------------------------------------
+void ComputeTrackingEfficiency(Double_t stEff[6], Double_t stEffErr[6][2], Double_t &spectroEff, Double_t spectroEffErr[2])
{
-
- TFile *file = new TFile(fileName, "read");
- TDirectoryFile *Directory = (TDirectoryFile*)((file->GetDirectory("MUON_Efficiency")));
-
- std::vector<TH1D*> CorrectedHistos;
- std::vector<TH1D*> CorrectedChiHisto;
- std::vector<TH1D*> CorrectedChjHisto;
+ /// compute the spectrometer efficiency and errors from the station efficiencies and errors
- TList *listTT = (TList *)Directory->Get("TotalTracksPerChamber;1");
- TList *listTD = (TList *)Directory->Get("TracksDetectedPerChamber;1");
- TList *listSD = (TList *)Directory->Get("SingleDetectedPerChamber;1");
-
- TH1D *EfficiencyHisto = 0x0;
- TH1D *TD = 0x0;
- TH1D *TT = 0x0;
+ Double_t de[6][2];
+ for (Int_t iSt = 0; iSt < 6; iSt++) de[iSt][0] = stEff[iSt]*stEff[iSt];
- for (Int_t Chamber = 1; Chamber < 12; Chamber++) //Loop over all chambers + 1
- {
-
- TH2F *TD2Dim = (TH2F *) listTD->At(Chamber - 1);
- TH2F *TT2Dim = (TH2F *) listTT->At(Chamber - 1);
- TH2F *SD2Dim = (TH2F *) listSD->At(Chamber - 1);
-
- Int_t nDEch = 0;
- nDEch = TD2Dim->GetNbinsX();
-
- TH1D *Ncluster = 0x0;
- TH1D *Corrected = 0x0;
- TH1D *SD = 0x0;
-
- TH1D *EfficiencyHistoCh11 = 0x0;
- TH1D *TDch11 = 0x0;
- TH1D *TTch11 = 0x0;
-
- if (Chamber == 11)
- {
-
- TDch11 = new TH1D(Form("TDch11 %d",Chamber), "", nDEch, - 0.5, nDEch - 0.5);
- TTch11 = new TH1D(Form("TTch11 %d",Chamber), "", nDEch, - 0.5, nDEch - 0.5);
-
+ if (moreTrackCandidates) {
+
+ spectroEff = stEff[0] * stEff[1] * stEff[2] * stEff[3] * stEff[4];
+
+ for (Int_t i = 0; i < 2; i++) {
+
+ for (Int_t iSt = 0; iSt < 6; iSt++) de[iSt][1] = stEffErr[iSt][i]*stEffErr[iSt][i];
+
+ spectroEffErr[i] = 0.;
+ for (Int_t j = 1; j < 32; j++) {
+ Double_t sigmaAdd = 1.;
+ for (Int_t iSt = 0; iSt < 5; iSt++) sigmaAdd *= de[iSt][TESTBIT(j,iSt)];
+ spectroEffErr[i] += sigmaAdd;
}
-
- TD = new TH1D(Form("TD %d",Chamber), "", nDEch, - 0.5, nDEch - 0.5);
- TT = new TH1D(Form("TT %d",Chamber), "", nDEch, - 0.5, nDEch - 0.5);
- SD = new TH1D(Form("SD %d",Chamber), "", nDEch, - 0.5, nDEch - 0.5);
-
- std::vector<Double_t> TDentries (nDEch, 0.0);
- std::vector<Double_t> TTentries (nDEch, 0.0);
- std::vector<Double_t> SDentries (nDEch, 0.0);
-
- for (Int_t ix = 1; ix <= nDEch; ix++)
- {
- for (Int_t iy = 1; iy <= TD2Dim->GetNbinsY(); iy++)
- {
-
- TDentries[ix-1] += TD2Dim->GetBinContent(ix,iy);
- TTentries[ix-1] += TT2Dim->GetBinContent(ix,iy);
- SDentries[ix-1] += SD2Dim->GetBinContent(ix,iy);
-
- }
-
- TD->SetBinContent(ix,TDentries[ix-1]);
- TD->SetBinError(ix,sqrt(TDentries[ix-1]));
- TT->SetBinContent(ix,TTentries[ix-1]);
- TT->SetBinError(ix,sqrt(TTentries[ix-1]));
- SD->SetBinContent(ix,SDentries[ix-1]);
- SD->SetBinError(ix,sqrt(SDentries[ix-1]));
-
- if (Chamber != 11) continue;
-
- TDch11->SetBinContent(ix,TDentries[ix-1]);
- TDch11->SetBinError(ix,sqrt(TDentries[ix-1]));
- TTch11->SetBinContent(ix,TTentries[ix-1]);
- TTch11->SetBinError(ix,sqrt(TTentries[ix-1]));
-
- }
+ spectroEffErr[i] = TMath::Sqrt(spectroEffErr[i]);
- TGraphAsymmErrors *EfficiencyGraph = new TGraphAsymmErrors(TD, TT);
-
- EfficiencyHisto = GraphToHist(kFALSE,EfficiencyGraph,Chamber);
-
- Ncluster = new TH1D(Form("Ncluster %d",Chamber), "N_{Cluster}", nDEch, - 0.5, nDEch - 0.5);
- Ncluster->Sumw2();
- Ncluster->Add(TD,SD,1,1);
-
- Corrected = new TH1D(Form("Corrected %d",Chamber), Form("Corrected N_{Cluster} Ch %d",Chamber), nDEch, - 0.5, nDEch - 0.5);
- Corrected->Sumw2();
- Corrected->Divide(Ncluster,EfficiencyHisto,1,1);
-
- CorrectedHistos.push_back(Corrected);
-
- if (Chamber != 11) continue;
-
- TGraphAsymmErrors *EfficiencyGraphCh11 = new TGraphAsymmErrors(TDch11, TTch11);
-
- EfficiencyHistoCh11 = GraphToHist(kTRUE,EfficiencyGraphCh11,Chamber);
-
- TCanvas *CanvasEffPerCh = new TCanvas("Efficiency per chambers", "Efficiency per chambers");
- gStyle->SetOptStat(0);
- CanvasEffPerCh->SetFillColor(0);
- CanvasEffPerCh->SetBorderSize(1);
- CanvasEffPerCh->SetFrameLineWidth(2);
-
- EfficiencyHistoCh11->SetTitle(kFALSE);
- EfficiencyHistoCh11->SetMarkerStyle(20);
- EfficiencyHistoCh11->SetMarkerColor(1);
- EfficiencyHistoCh11->SetMarkerSize(1.3);
- EfficiencyHistoCh11->GetXaxis()->SetTitle("Chamber");
- EfficiencyHistoCh11->GetXaxis()->SetNdivisions(Chamber);
- EfficiencyHistoCh11->GetXaxis()->SetTitleSize(0.05);
- EfficiencyHistoCh11->GetXaxis()->SetTitleOffset(0.91);
- EfficiencyHistoCh11->GetXaxis()->SetLabelSize(0.05);
- EfficiencyHistoCh11->GetYaxis()->SetTitle("Efficiency");
- EfficiencyHistoCh11->GetYaxis()->CenterTitle(kTRUE);
- EfficiencyHistoCh11->GetYaxis()->SetTitleSize(0.06);
- EfficiencyHistoCh11->GetYaxis()->SetTitleOffset(0.79);
- EfficiencyHistoCh11->GetYaxis()->SetLabelSize(0.05);
- EfficiencyHistoCh11->Draw("e");
-
}
-
- /*************************Correlated Dead Arrays check****************************/
-
-
- Double_t MissingTracks [5] = {0.0};//Station one is in slot cero!!!
-
- for (Int_t CDAperSt = 0; CDAperSt < 5; CDAperSt++)
- {
-
- std::vector<Double_t> FilledVectorChi = FillVector(CorrectedHistos[CDAperSt*2], kFALSE);
- std::vector<Double_t> FilledVectorChj = FillVector(CorrectedHistos[(CDAperSt*2)+1], kFALSE);
- std::vector<Double_t> FilledVectorChiErr = FillVector(CorrectedHistos[CDAperSt*2], kTRUE);
- std::vector<Double_t> FilledVectorChjErr = FillVector(CorrectedHistos[(CDAperSt*2)+1], kTRUE);
-
- Int_t DEs = FilledVectorChi.size();
- Int_t SimGroups = 0;
- Double_t ToFillMissingTracks = 0.0;
+
+ } else {
+
+ spectroEff = stEff[0] * stEff[1] * stEff[2] * stEff[5];
+
+ for (Int_t i = 0; i < 2; i++) {
+
+ for (Int_t iSt = 0; iSt < 6; iSt++) de[iSt][1] = stEffErr[iSt][i]*stEffErr[iSt][i];
+
+ spectroEffErr[i] = 0.;
+ for (Int_t j = 1; j < 16; j++) {
+ Double_t sigmaAdd = de[5][TESTBIT(j,3)];
+ for (Int_t iSt = 0; iSt < 3; iSt++) sigmaAdd *= de[iSt][TESTBIT(j,iSt)];
+ spectroEffErr[i] += sigmaAdd;
+ }
+ spectroEffErr[i] = TMath::Sqrt(spectroEffErr[i]);
- if ( CDAperSt < 2 ) SimGroups = static_cast<Int_t>(DEs/2);
- else SimGroups = static_cast<Int_t>(DEs/2) + 1;
-
- std::vector<Double_t> CDAs = CorrelatedDeadArray(FilledVectorChi, FilledVectorChj, DEs, SimGroups);
- Drawing(FilledVectorChi, FilledVectorChj, FilledVectorChiErr, FilledVectorChjErr, DEs, SimGroups, CDAperSt);
-
- for (Int_t DEst = 0; DEst < DEs; DEst++)
- {
- if ( CDAs[DEst] != 0.0 )
- {
-
- printf ("Looks like there is a CDA in DE %d of Station %d: N00 = %.2f \n", DEst, CDAperSt+1, CDAs[DEst]);
- ToFillMissingTracks += CDAs[DEst];
-
- }
- }
-
- MissingTracks[CDAperSt] = ToFillMissingTracks;
-
}
-
-
-
- /**********************************Efficiency from Efficiency Task************************************/
+
+ }
+ spectroEffErr[0] = TMath::Min(spectroEff, spectroEffErr[0]);
+ spectroEffErr[1] = TMath::Min(1.-spectroEff, spectroEffErr[1]);
- vector<Double_t> chambersEff;
- vector<Double_t> chambersEffErr;
- vector<Double_t> StEfficiency;
- vector<Double_t> StEfficiencyErr;
-
- chambersEff.push_back(0.0);
- chambersEffErr.push_back(0.0);
- StEfficiency.push_back(0.0);
- StEfficiencyErr.push_back(0.0);
-
- for (Int_t EfCh = 0; EfCh < 10; EfCh++)
- {
-
- chambersEff.push_back(EfficiencyHisto->GetBinContent(EfCh+1));
- chambersEffErr.push_back(EfficiencyHisto->GetBinError(EfCh+1));
-
- cout << "Efficiency chamber " << EfCh+1 << " = " << chambersEff[EfCh+1] << " +- " << chambersEffErr[EfCh+1] << endl;
-
- }
-
- for (Int_t EffSt = 0; EffSt < 4; EffSt++)
- {
-
- if (EffSt < 3)
- {
-
- StEfficiency.push_back(1.0 - (1.0-chambersEff[(2*EffSt)+1])*(1.0-chambersEff[(2*EffSt)+2]));
- StEfficiencyErr.push_back(StEfficiency[EffSt+1]*(chambersEffErr[(2*EffSt)+1]/chambersEff[(2*EffSt)+1] + chambersEffErr[(2*EffSt)+2]/chambersEff[(2*EffSt)+2]));
-
- cout << "Efficiency station " << EffSt+1 <<" = " <<StEfficiency[EffSt+1] << " +- " <<StEfficiencyErr[EffSt+1]<<endl;
-
- }
-
- else
- {
-
- StEfficiency.push_back(chambersEff[(2*EffSt)+1] * chambersEff[(2*EffSt)+2] * chambersEff[(2*EffSt)+3] * chambersEff[(2*EffSt)+4]
- + (1.0 - chambersEff[(2*EffSt)+1]) * chambersEff[(2*EffSt)+2] * chambersEff[(2*EffSt)+3] * chambersEff[(2*EffSt)+4]
- + chambersEff[(2*EffSt)+1] * (1.0 - chambersEff[(2*EffSt)+2]) * chambersEff[(2*EffSt)+3] * chambersEff[(2*EffSt)+4]
- + chambersEff[(2*EffSt)+1] * chambersEff[(2*EffSt)+2] * (1.0 - chambersEff[(2*EffSt)+3]) * chambersEff[(2*EffSt)+4]
- + chambersEff[(2*EffSt)+1] * chambersEff[(2*EffSt)+2] * chambersEff[(2*EffSt)+3] * (1.0 - chambersEff[(2*EffSt)+4]));
-
- StEfficiencyErr.push_back(StEfficiency[EffSt+1]*(chambersEffErr[(2*EffSt)+1]/chambersEff[(2*EffSt)+1]
- + chambersEffErr[(2*EffSt)+2]/chambersEff[(2*EffSt)+2] + chambersEffErr[(2*EffSt)+3]/chambersEff[(2*EffSt)+3]
- + chambersEffErr[(2*EffSt)+4]/chambersEff[(2*EffSt)+4]));
-
- cout << "Efficiency station 45"<<" = " <<StEfficiency[EffSt+1] << " +- " <<StEfficiencyErr[EffSt+1]<<endl;
- //Notice efficiency from st 45 is in slot 4 of the vector!!!
-
- }
-
- }
-
- StEfficiency.push_back(StEfficiency[1]*StEfficiency[2]*StEfficiency[3]*StEfficiency[4]);
- StEfficiencyErr.push_back(StEfficiencyErr[1]+StEfficiencyErr[2]+StEfficiencyErr[3]+StEfficiencyErr[4]);
-
- cout << "Total efficiency = " <<StEfficiency[5]<< " +- " <<StEfficiencyErr[5]<< endl;
-
-
+}
- /*************************************Correcting efficiencies********************************************/
+//---------------------------------------------------------------------------
+void GetTrackingEfficiency(TArrayD &chEff, TArrayD chEffErr[2], TGraphAsymmErrors *effVsSt[3],
+ TGraphAsymmErrors *effVsX[3], Int_t ip, Double_t xp, Bool_t print)
+{
+ /// compute the tracking efficiency from the individual chamber efficiencies
- TH1D *TTcorrected = new TH1D("TTcorrected", "TDcorrected", 10, -0.5, 9.5);
-
- for (Int_t St = 0; St < 5; St++)
- {
-
- Double_t NewBinContentChi = TT->GetBinContent((2*St)+1) + MissingTracks[St];
- TTcorrected->SetBinContent((2*St)+1,NewBinContentChi);
- TTcorrected->SetBinContent((2*St)+2,TT->GetBinContent((2*St)+2));
-
+ Double_t stEff[6];
+ Double_t stEffErr[6][2];
+
+ // check if unknown chamber efficiency
+ Bool_t allEffKnown = kTRUE;
+ for (Int_t iCh = 1; iCh < 11 && allEffKnown; ++iCh) if (chEff[iCh] < 0) allEffKnown = kFALSE;
+
+ if (allEffKnown) {
+
+ // compute station efficiency
+ for (Int_t iSt = 0; iSt < 5; ++iSt) ComputeStationEfficiency(chEff, chEffErr, iSt, stEff[iSt], stEffErr[iSt]);
+ ComputeStation45Efficiency(chEff, chEffErr, stEff[5], stEffErr[5]);
+
+ // compute spectrometer efficiency
+ Double_t spectroEff, spectroEffErr[2];
+ ComputeTrackingEfficiency(stEff, stEffErr, spectroEff, spectroEffErr);
+ chEff[0] = spectroEff;
+ chEffErr[0][0] = spectroEffErr[0];
+ chEffErr[1][0] = spectroEffErr[1];
+
+ // fill graphs if required
+ for (Int_t i = 0; i < 3; ++i) {
+ if (effVsX[i]) {
+ effVsX[i]->SetPoint(ip,xp,chEff[0]);
+ effVsX[i]->SetPointError(ip,0.,0.,chEffErr[0][0],chEffErr[1][0]);
+ }
+ if (effVsSt[i]) for (Int_t iSt = 0; iSt < 6; ++iSt) {
+ effVsSt[i]->SetPoint(iSt,iSt+1,stEff[iSt]);
+ effVsSt[i]->SetPointError(iSt,0.,0.,stEffErr[iSt][0],stEffErr[iSt][1]);
+ }
}
-
- TGraphAsymmErrors *EfficiencyGraphCorrected = new TGraphAsymmErrors(TD,TTcorrected);
-
- TH1D *EfficiencyHistoCorrected = GraphToHist(1,EfficiencyGraphCorrected,11);//In this case doesnt matter which value is used in the first slot.
-
- vector<Double_t> chambersEffCorrected;
- vector<Double_t> chambersEffErrCorrected;
-
- chambersEffCorrected.push_back(0.0);
- chambersEffErrCorrected.push_back(0.0);
-
- for (Int_t EfChCorr = 0; EfChCorr < 10; EfChCorr++)
- {
-
- chambersEffCorrected.push_back(EfficiencyHistoCorrected->GetBinContent(EfChCorr+1));
- chambersEffErrCorrected.push_back(EfficiencyHistoCorrected->GetBinError(EfChCorr+1));
-
- if ( (chambersEffCorrected[EfChCorr+1]/chambersEff[EfChCorr+1]) == 1 ) continue;
-
- cout << "Corrected efficiency chamber " << EfChCorr+1 << " = " << chambersEffCorrected[EfChCorr+1] << " +- " << chambersEffErrCorrected[EfChCorr+1] << endl;
-
+
+ } else {
+
+ Double_t edge[2] = {0., 1.};
+ TArrayD chEffEdge[2];
+ Double_t stEffEdge[2][6];
+ Double_t stEffEdgeErr[2][6][2];
+ Double_t spectroEffEdge[2], spectroEffEdgeErr[2][2];
+
+ for (Int_t i = 0; i < 2; ++i) {
+
+ // set lower(upper) limit of chamber efficiency
+ chEffEdge[i].Set(11);
+ for (Int_t iCh = 1; iCh < 11; ++iCh) chEffEdge[i][iCh] = (chEff[iCh] < 0) ? edge[i] : chEff[iCh];
+
+ // compute station efficiency
+ for (Int_t iSt = 0; iSt < 5; ++iSt) ComputeStationEfficiency(chEffEdge[i], chEffErr, iSt, stEffEdge[i][iSt], stEffEdgeErr[i][iSt]);
+ ComputeStation45Efficiency(chEffEdge[i], chEffErr, stEffEdge[i][5], stEffEdgeErr[i][5]);
+
+ // compute spectrometer efficiency
+ ComputeTrackingEfficiency(stEffEdge[i], stEffEdgeErr[i], spectroEffEdge[i], spectroEffEdgeErr[i]);
+
+ // fill graph if required
+ if (effVsX[i+1]) {
+ effVsX[i+1]->SetPoint(ip,xp,spectroEffEdge[i]);
+ effVsX[i+1]->SetPointError(ip,0.,0.,spectroEffEdgeErr[i][0],spectroEffEdgeErr[i][1]);
+ }
+ if (effVsSt[i+1]) for (Int_t iSt = 0; iSt < 6; ++iSt) {
+ effVsSt[i+1]->SetPoint(iSt,iSt+1,stEffEdge[i][iSt]);
+ effVsSt[i+1]->SetPointError(iSt,0.,0.,stEffEdgeErr[i][iSt][0],stEffEdgeErr[i][iSt][1]);
+ }
+
}
-
- vector<Double_t> StEfficiencyCorrected;
- vector<Double_t> StEfficiencyErrCorrected;
-
- StEfficiencyCorrected.push_back(0.0);
- StEfficiencyErrCorrected.push_back(0.0);
-
- //By default CDA in Station 4 is not corrected.
+
+ // combine extreme cases to get station efficiency
+ for (Int_t iSt = 0; iSt < 6; ++iSt) {
+ stEff[iSt] = stEffEdge[1][iSt];
+ stEffErr[iSt][0] = stEffEdge[1][iSt] - stEffEdge[0][iSt] + stEffEdgeErr[0][iSt][0];
+ stEffErr[iSt][1] = stEffEdgeErr[1][iSt][1];
+ }
+
+ // combine extreme cases to get spectrometer efficiency
+ chEff[0] = spectroEffEdge[1];
+ chEffErr[0][0] = spectroEffEdge[1] - spectroEffEdge[0] + spectroEffEdgeErr[0][0];
+ chEffErr[1][0] = spectroEffEdgeErr[1][1];
+
+ // fill graph if required
+ if (effVsX[0]) {
+ effVsX[0]->SetPoint(ip,xp,chEff[0]);
+ effVsX[0]->SetPointError(ip,0.,0.,chEffErr[0][0],chEffErr[1][0]);
+ }
+ if (effVsSt[0]) for (Int_t iSt = 0; iSt < 6; ++iSt) {
+ effVsSt[0]->SetPoint(iSt,iSt+1,stEff[iSt]);
+ effVsSt[0]->SetPointError(iSt,0.,0.,stEffErr[iSt][0],stEffErr[iSt][1]);
+ }
+
+ }
- for (Int_t Station = 0; Station < 4; Station++)
- {
-
- if (MissingTracks[Station] == 0.0)
- {
-
- StEfficiencyCorrected.push_back(StEfficiency[Station+1]);
- StEfficiencyErrCorrected.push_back(StEfficiencyErr[Station+1]);
-
- }
-
- else
- {
-
- if (Station < 3)
- {
-
- StEfficiencyCorrected.push_back(StEfficiency[Station+1]*(chambersEffCorrected[(2*Station)+1]/chambersEff[(2*Station)+1]));
- StEfficiencyErrCorrected.push_back(StEfficiencyCorrected[Station+1]*sqrt((StEfficiencyErr[Station+1]/StEfficiency[Station+1])*(StEfficiencyErr[Station+1]/StEfficiency[Station+1])+(chambersEffErrCorrected[(2*Station)+1]/chambersEffCorrected[(2*Station)+1])*(chambersEffErrCorrected[(2*Station)+1]/chambersEffCorrected[(2*Station)+1]) + (chambersEffErr[(2*Station)+1]/chambersEff[(2*Station)+1])*(chambersEffErr[(2*Station)+1]/chambersEff[(2*Station)+1])));
-
- cout<<"Corrected efficiency station "<<Station+1<<" = "<<StEfficiencyCorrected[Station+1]<<" +- "<<StEfficiencyErrCorrected[Station+1]<<endl;
-
- }
-
- if (Station == 3)
- {
-
- StEfficiencyCorrected.push_back(StEfficiency[Station+1]*(chambersEffCorrected[(2*Station)+3]/chambersEff[(2*Station)+3]));
- StEfficiencyErrCorrected.push_back(StEfficiencyCorrected[Station+1]*sqrt((StEfficiencyErr[Station+1]/StEfficiency[Station+1])*(StEfficiencyErr[Station+1]/StEfficiency[Station+1])+(chambersEffErrCorrected[(2*Station)+3]/chambersEffCorrected[(2*Station)+3])*(chambersEffErrCorrected[(2*Station)+3]/chambersEffCorrected[(2*Station)+3]) + (chambersEffErr[(2*Station)+3]/chambersEff[(2*Station)+3])*(chambersEffErr[(2*Station)+3]/chambersEff[(2*Station)+3])));
-
- cout<<"Corrected efficiency station 45"<<" = "<<StEfficiencyCorrected[Station+1]<<" +- "<<StEfficiencyErrCorrected[Station+1]<<endl;
-
- }
-
-
- }
-
+ // print results
+ if (print) {
+ for (Int_t iCh = 1; iCh < 11; ++iCh) {
+ cout << "Efficiency chamber " << iCh << " : ";
+ cout << chEff[iCh] << " + " << chEffErr[1][iCh] << " - " << chEffErr[0][iCh] << endl;
}
-
- StEfficiencyCorrected.push_back(StEfficiencyCorrected[1]*StEfficiencyCorrected[2]*StEfficiencyCorrected[3]*StEfficiencyCorrected[4]);
- StEfficiencyErrCorrected.push_back(StEfficiencyErrCorrected[1]+StEfficiencyErrCorrected[2]+StEfficiencyErrCorrected[3]+StEfficiencyErrCorrected[4]);
-
- cout<<"Corrected total efficiency"<<" = "<<StEfficiencyCorrected[5]<< " +- " <<StEfficiencyErrCorrected[5]<< endl;
-
+ for (Int_t iSt = 0; iSt < 6; ++iSt) {
+ if (iSt < 5) cout << "Station " << iSt+1 << " = ";
+ else cout << "Station 45 = ";
+ cout << stEff[iSt] << " + " << stEffErr[iSt][1] << " - " << stEffErr[iSt][0] << endl;
+ }
+ cout << "Total tracking efficiency : ";
+ cout << chEff[0] << " + " << chEffErr[1][0] << " - " << chEffErr[0][0] << endl << endl;
+ }
+
}
-/*************************************Functions************************************/
-
-
-std::vector<Double_t> FillVector(TH1D* Chamber, Bool_t Errors)
+//---------------------------------------------------------------------------
+void IntegrateMuonEfficiency(TGraphAsymmErrors &effVsRunLow, TGraphAsymmErrors &effVsRunUp,
+ TGraphAsymmErrors &effVsX, Int_t ip, Double_t xp)
{
-
- Int_t NumOfDEs = Chamber->GetNbinsX();
- std::vector<Double_t> ChVector (NumOfDEs, 0);
- for (Int_t i = 1; i <= NumOfDEs; i++)
- {
- if (Errors == kFALSE) ChVector.at(i-1) = Chamber->GetBinContent(i);
- if (Errors == kTRUE) ChVector.at(i-1) = Chamber->GetBinError(i);
- }
-
- return ChVector;
-
+ /// integrate efficiency over runs
+ /// return kFALSE if efficiency unknown in all runs
+
+ if (!runWeights) {
+ effVsX.SetPoint(ip,xp,-1.);
+ effVsX.SetPointError(ip,0.,0.,0.,0.);
+ return;
+ }
+
+ // initialize
+ Double_t rec[2] = {0., 0.};
+ Double_t gen = 0.;
+ Double_t intEffErr[2] = {0., 0.};
+ Bool_t ok = kFALSE;
+
+ // loop over runs
+ Int_t nRuns = effVsRunLow.GetN();
+ for (Int_t iRun = 0; iRun < nRuns; ++iRun) {
+
+ // get run weight
+ TString sRun = effVsRunLow.GetXaxis()->GetBinLabel(iRun+1);
+ TParameter<Double_t> *weight = static_cast<TParameter<Double_t>*>(runWeights->FindObject(sRun.Data()));
+ if (!weight) {
+ printf("weight not found for run %s\n", sRun.Data());
+ continue;
+ }
+ Double_t w = weight->GetVal();
+ Double_t w2 = w*w;
+
+ // get efficiency and error
+ Double_t eff[2] = {effVsRunLow.GetY()[iRun], effVsRunUp.GetY()[iRun]};
+ Double_t effErr[2] = {effVsRunLow.GetErrorYlow(iRun), effVsRunUp.GetErrorYhigh(iRun)};
+ if (eff[0] < 0. || eff[1] < 0.) {
+ printf("no efficiency measurement --> use 0(1) ± 0 as lower(upper) limit\n");
+ eff[0] = 0.;
+ eff[1] = 1.;
+ effErr[0] = 0.;
+ effErr[1] = 0.;
+ } else ok = kTRUE;
+
+ // integrate
+ gen += w;
+ for (Int_t i = 0; i < 2; ++i) {
+ rec[i] += w*eff[i];
+ intEffErr[i] += w2*effErr[i]*effErr[i];
+ }
+
+ }
+
+ // compute efficiency
+ if (gen > 0. && ok) {
+
+ effVsX.SetPoint(ip,xp,rec[1]/gen);
+ effVsX.SetPointError(ip,0.,0.,(rec[1]-rec[0]+TMath::Sqrt(intEffErr[0]))/gen,TMath::Sqrt(intEffErr[1])/gen);
+
+ } else {
+
+ if (gen <= 0.) printf("impossible to integrate, all weights = 0 or unknown ?!?\n");
+ else printf("efficiency never measured --> return -1 ± 0\n");
+
+ effVsX.SetPoint(ip,xp,-1.);
+ effVsX.SetPointError(ip,0.,0.,0.,0.);
+
+ }
+
}
-std::vector<Double_t> CorrelatedDeadArray(std::vector<Double_t> Chi, std::vector<Double_t> Chj, Int_t NumberOfDEs, Int_t SymmetricGroups)
-{
-
- std::vector<Double_t> Average (SymmetricGroups, 0.0);
- std::vector<Double_t> CDA (NumberOfDEs, 0.0);
-
- for (Int_t Sim = 0; Sim < SymmetricGroups/2; Sim++)
- {
- if (Sim == 0)
- {
-
- if ( (Sim == 0) && (SymmetricGroups/2 == 1) )
- {
-
- Average[Sim] = (Chi[Sim]+Chi[Sim+1]+Chi[Sim+2]+Chi[Sim+3])/4;
- Average[Sim+(SymmetricGroups/2)] = (Chj[Sim]+Chj[Sim+1]+Chj[Sim+2]+Chj[Sim+3])/4;
-
- for (Int_t i = 1; i < 5; i++)
- {
-
- if ( (Chi[i-1] < (Average[Sim]-(Average[Sim]*0.2))) && (Chj[i-1] < (Average[Sim+(SymmetricGroups/2)]-(Average[Sim+(SymmetricGroups/2)]*0.2))) )
- {
- CDA[i-1] = N00(Average[Sim], Chi[i-1], Sim, SymmetricGroups/2);
- }
-
- }
-
- }
-
- else
- {
-
- Average[Sim] = (Chi[Sim] + Chi[Sim+(NumberOfDEs/2)])/2;
- Average[Sim+(SymmetricGroups/2)] = (Chj[Sim] + Chj[Sim+(NumberOfDEs/2)])/2;
-
-
- if ( (Chi[Sim] < (Average[Sim]-(Average[Sim]*0.2))) && (Chj[Sim] < (Average[Sim+(SymmetricGroups/2)]-(Average[Sim+(SymmetricGroups/2)]*0.2))) )
- {
- CDA[Sim] = N00(Average[Sim], Chi[Sim], Sim, SymmetricGroups/2);
- }
- if ( (Chi[Sim+(NumberOfDEs/2)] < (Average[Sim]-(Average[Sim]*0.2))) && (Chj[Sim+(NumberOfDEs/2)] < (Average[Sim+(SymmetricGroups/2)]-(Average[Sim+(SymmetricGroups/2)]*0.2))) )
- {
- CDA[Sim+(NumberOfDEs/2)] = N00(Average[Sim], Chi[Sim+(NumberOfDEs/2)], Sim, SymmetricGroups/2);
- }
-
- }
-
- }
-
- else
- {
-
- Average[Sim] = (Chi[Sim] + Chi[(NumberOfDEs/2)-Sim] + Chi[Sim+(NumberOfDEs/2)] + Chi[NumberOfDEs-Sim])/4;
- Average[Sim+(SymmetricGroups/2)] = (Chj[Sim] + Chj[(NumberOfDEs/2)-Sim] + Chj[Sim+(NumberOfDEs/2)] + Chj[NumberOfDEs-Sim])/4;
- if ( (Chi[Sim] < (Average[Sim]-(Average[Sim]*0.2))) && (Chj[Sim] < (Average[Sim+(SymmetricGroups/2)]-(Average[Sim+(SymmetricGroups/2)]*0.2))) )
- {
- CDA[Sim] = N00(Average[Sim], Chi[Sim], Sim, SymmetricGroups/2);
- }
- if ( (Chi[(NumberOfDEs/2)-Sim] < (Average[Sim]-(Average[Sim]*0.2))) && (Chj[(NumberOfDEs/2)-Sim] < (Average[Sim+(SymmetricGroups/2)]-(Average[Sim+(SymmetricGroups/2)]*0.2))) )
- {
- CDA[(NumberOfDEs/2)-Sim] = N00(Average[Sim], Chi[(NumberOfDEs/2)-Sim], Sim, SymmetricGroups/2);
- }
- if ( (Chi[Sim+(NumberOfDEs/2)] < (Average[Sim]-(Average[Sim]*0.2))) && (Chj[Sim+(NumberOfDEs/2)] < (Average[Sim+(SymmetricGroups/2)]-(Average[Sim+(SymmetricGroups/2)]*0.2))) )
- {
- CDA[Sim+(NumberOfDEs/2)] = N00(Average[Sim], Chi[Sim+(NumberOfDEs/2)], Sim, SymmetricGroups/2);
- }
-
- if ( (Chi[NumberOfDEs-Sim] < (Average[Sim]-(Average[Sim]*0.2))) && (Chj[NumberOfDEs-Sim] < (Average[Sim+(SymmetricGroups/2)]-(Average[Sim+(SymmetricGroups/2)]*0.2))) )
- {
- CDA[NumberOfDEs-Sim] = N00(Average[Sim], Chi[NumberOfDEs-Sim], Sim, SymmetricGroups/2);
-
- }
-
- }
-
+//---------------------------------------------------------------------------
+void LoadRunWeights(TString fileName)
+{
+ /// Set the number of interested events per run
+ /// (used to weight the acc*eff correction integrated
+ /// over run for any pt/y/centrality bins)
+
+ if (runWeights) return;
+
+ ifstream inFile(fileName.Data());
+ if (!inFile.is_open()) {
+ printf("cannot open file %s", fileName.Data());
+ return;
+ }
+
+ runWeights = new THashList(1000);
+ runWeights->SetOwner();
+
+ TString line;
+ while (! inFile.eof() ) {
+
+ line.ReadLine(inFile,kTRUE);
+ if(line.IsNull()) continue;
+
+ TObjArray *param = line.Tokenize(" ");
+ if (param->GetEntries() != 2) {
+ printf("bad input line %s", line.Data());
+ continue;
}
-
- return CDA;
-
+
+ Int_t run = ((TObjString*)param->UncheckedAt(0))->String().Atoi();
+ if (run < 0) {
+ printf("invalid run number: %d", run);
+ continue;
+ }
+
+ Float_t weight = ((TObjString*)param->UncheckedAt(1))->String().Atof();
+ if (weight <= 0.) {
+ printf("invalid weight: %g", weight);
+ continue;
+ }
+
+ runWeights->Add(new TParameter<Double_t>(Form("%d",run), weight));
+
+ delete param;
+ }
+
+ inFile.close();
+
}
-Double_t N00(Double_t Average, Double_t CDA, Int_t SimGroup, Int_t SimGroupsPerCh)
+//---------------------------------------------------------------------------
+void SetCentPtCh(THnSparse& SparseData)
{
-
- Double_t RecomputedAverage = 0.0;
- Double_t ToBeAdded = 0.0;
-
- if ( (SimGroup == 0) && (SimGroupsPerCh != 1) ) RecomputedAverage = (2*Average) - CDA;
- else RecomputedAverage = ((4*Average) - CDA)/3;
-
- ToBeAdded = RecomputedAverage - CDA;
-
- return ToBeAdded;
-
+ /// Sets the centrality and pT range and the charge for integration
+ /// If maxPt = -1, it sets maxPt as the maximum Pt value on the THnSparse
+
+ if (ptMax >= 0 && ptMax < ptMin) {
+ printf("Wrong Pt range, ptMax must be higher than ptMin\n");
+ exit(-1);
+ }
+
+ if (charge < -1 || charge > 1) {
+ printf("Selected charge must be 0, 1 or -1\n");
+ exit(-1);
+ }
+
+ // adjust centrality range
+ centMax = TMath::Max(centMax-1.e-12, centMin);
+
+ // set the centrality range for integration
+ Int_t lowBin = SparseData.GetAxis(1)->FindBin(centMin);
+ Int_t upBin = SparseData.GetAxis(1)->FindBin(centMax);
+ SparseData.GetAxis(1)->SetRange(lowBin, upBin);
+
+ // set the pt range for integration
+ lowBin = SparseData.GetAxis(2)->FindBin(ptMin);
+ if (ptMax == -1) { upBin = SparseData.GetAxis(2)->GetNbins()+1; }
+ else { upBin = SparseData.GetAxis(2)->FindBin(ptMax);}
+ SparseData.GetAxis(2)->SetRange(lowBin, upBin);
+
+ // set the charge range
+ if (charge != 0) {
+ lowBin = SparseData.GetAxis(5)->FindBin(charge);
+ SparseData.GetAxis(5)->SetRange(lowBin, lowBin);
+ }
+
}
-void Drawing(std::vector<Double_t> Chi, std::vector<Double_t> Chj, std::vector<Double_t> ChiErr, std::vector<Double_t> ChjErr, Int_t Bins, Int_t SymGroups, Int_t Station)
+//---------------------------------------------------------------------------
+TGraphAsymmErrors* CreateGraph(const char* name, const char* title, int value)
{
+ /// create a graph of a given name and title
+
+ TGraphAsymmErrors* g = new TGraphAsymmErrors();
+
+ if ( value >= 0 ) {
+
+ g->SetName(Form(name,value));
+ g->SetTitle(Form(title,value));
+
+ } else {
+
+ g->SetName(name);
+ g->SetTitle(title);
+
+
+ }
+
+ return g;
+
+}
- TCanvas *Canvas = new TCanvas(Form("Station %d",Station+1), Form("Station %d",Station+1));
- gStyle->SetOptStat(0);
- Canvas->SetFillColor(0);
- if ( (Station+1) < 3) Canvas->Divide(SymGroups/2);
- if ( (Station+1) == 3) Canvas->Divide(Station+1,Station);
- if ( (Station+1) > 3) Canvas->Divide(3,3);
- Canvas->SetBorderSize(1);
- Canvas->SetFrameLineWidth(2);
-
- for (Int_t Sim = 0; Sim < SymGroups/2; Sim++)
- {
-
- TH1D *Chamberi = new TH1D(Form("Chi SimGroup %d St %d",Sim,Station+1), "", Bins, - 0.5, Bins - 0.5);
- TGaxis::SetMaxDigits(2);
- Chamberi->SetTitle(kFALSE);
- Chamberi->SetMarkerStyle(20);
- Chamberi->SetMarkerColor(1);
- Chamberi->SetMarkerSize(1.3);
- Chamberi->GetXaxis()->SetTitle("Detection Element");
- Chamberi->GetXaxis()->SetNdivisions(Bins+1);
- Chamberi->GetXaxis()->SetTitleSize(0.05);
- Chamberi->GetXaxis()->SetTitleOffset(0.91);
- Chamberi->GetXaxis()->SetLabelSize(0.05);
- Chamberi->GetYaxis()->SetTitle("Corrected N_{Cluster}");
- Chamberi->GetYaxis()->CenterTitle(kTRUE);
- Chamberi->GetYaxis()->SetTitleSize(0.06);
- Chamberi->GetYaxis()->SetTitleOffset(0.79);
- Chamberi->GetYaxis()->SetLabelSize(0.05);
- TH1D *Chamberj = new TH1D(Form("Chj SimGroup %d St %d",Sim,Station+1), "", Bins, - 0.5, Bins - 0.5);
- Chamberj->SetTitle(kFALSE);
- Chamberj->SetMarkerStyle(20);
- Chamberj->SetMarkerColor(2);
- Chamberj->SetMarkerSize(1.3);
-
- if (Sim == 0)
- {
-
- if ( (Sim == 0) && (SymGroups/2 == 1) )
- {
-
- for (Int_t i = 1; i < 5; i++)
- {
-
- Chamberi->SetBinContent(i, Chi[i-1]);
- Chamberi->SetBinError(i, ChiErr[i-1]);
- Chamberj->SetBinContent(i, Chj[i-1]);
- Chamberj->SetBinError(i, ChjErr[i-1]);
-
- }
-
- }
-
- else
- {
-
- Chamberi->SetBinContent(Sim+1, Chi[Sim]);
- Chamberi->SetBinContent((Sim+1)+(Bins/2), Chi[Sim+(Bins/2)]);
- Chamberi->SetBinError(Sim+1, ChiErr[Sim]);
- Chamberi->SetBinError((Sim+1)+(Bins/2), ChiErr[Sim+(Bins/2)]);
- Chamberj->SetBinContent(Sim+1, Chj[Sim]);
- Chamberj->SetBinContent((Sim+1)+(Bins/2), Chj[Sim+(Bins/2)]);
- Chamberj->SetBinError(Sim+1, ChjErr[Sim]);
- Chamberj->SetBinError((Sim+1)+(Bins/2), ChjErr[Sim+(Bins/2)]);
-
- }
-
- }
-
- else
- {
-
- Chamberi->SetBinContent(Sim+1, Chi[Sim]);
- Chamberi->SetBinContent((Bins/2)-(Sim-1), Chi[(Bins/2)-Sim]);
- Chamberi->SetBinContent((Sim+1)+(Bins/2), Chi[Sim+(Bins/2)]);
- Chamberi->SetBinContent(Bins-(Sim-1), Chi[Bins-Sim]);
-
- Chamberi->SetBinError(Sim+1, ChiErr[Sim]);
- Chamberi->SetBinError((Bins/2)-(Sim-1), ChiErr[(Bins/2)-Sim]);
- Chamberi->SetBinError((Sim+1)+(Bins/2), ChiErr[Sim+(Bins/2)]);
- Chamberi->SetBinError(Bins-(Sim-1), ChiErr[Bins-Sim]);
-
- Chamberj->SetBinContent(Sim+1, Chj[Sim]);
- Chamberj->SetBinContent((Bins/2)-(Sim-1), Chj[(Bins/2)-Sim]);
- Chamberj->SetBinContent((Sim+1)+(Bins/2), Chj[Sim+(Bins/2)]);
- Chamberj->SetBinContent(Bins-(Sim-1), Chj[Bins-Sim]);
- Chamberj->SetBinError(Sim+1, ChjErr[Sim]);
- Chamberj->SetBinError((Bins/2)-(Sim-1), ChjErr[(Bins/2)-Sim]);
- Chamberj->SetBinError((Sim+1)+(Bins/2), ChjErr[Sim+(Bins/2)]);
- Chamberj->SetBinError(Bins-(Sim-1), ChjErr[Bins-Sim]);
+//---------------------------------------------------------------------------
+void BeautifyGraph(TGraphAsymmErrors &g, const char* xAxisName, const char* yAxisName)
+{
+ /// beautify this graph
+
+ g.SetLineStyle(1);
+ g.SetLineColor(1);
+ g.SetMarkerStyle(20);
+ g.SetMarkerSize(0.7);
+ g.SetMarkerColor(2);
+
+ g.GetXaxis()->SetTitle(xAxisName);
+ g.GetXaxis()->CenterTitle(kTRUE);
+ g.GetXaxis()->SetLabelFont(22);
+ g.GetXaxis()->SetTitleFont(22);
+
+ g.GetYaxis()->SetTitle(yAxisName);
+ g.GetYaxis()->CenterTitle(kTRUE);
+ g.GetYaxis()->SetLabelFont(22);
+ g.GetYaxis()->SetTitleFont(22);
+
+}
- }
- Canvas->cd(Sim+1);
- Chamberi->Draw("e");
- Chamberj->Draw("esame");
+//---------------------------------------------------------------------------
+void BeautifyGraphs(TObjArray& array, const char* xAxisName, const char* yAxisName)
+{
+ /// beautify all the graphs in this array
+
+ for ( Int_t i = 0; i <= array.GetLast(); ++i)
+ BeautifyGraph(*static_cast<TGraphAsymmErrors*>(array.UncheckedAt(i)), xAxisName, yAxisName);
+
+}
- TLegend *Legend = new TLegend (0.75, 0.8, 1.0, 1.0);
- Legend->SetBorderSize(0);
- Legend->SetFillColor(0);
- if ( (Station+1) < 3 ) Legend->SetTextSize(0.05);
- if ( (Station+1) > 2 ) Legend->SetTextSize(0.05);
- Legend->AddEntry(Chamberi,Form("Chamber %d",2*Station+1),"p");
- Legend->AddEntry(Chamberj,Form("Chamber %d",2*Station+2),"p");
- Legend->Draw();
-
- }
+//---------------------------------------------------------------------------
+void SetRunLabel(TGraphAsymmErrors &g, Int_t irun, const TList& runs)
+{
+ /// set the run label for all the graphs
+
+ g.GetXaxis()->Set(irun+1, -0.5, irun+0.5);
+
+ TIter nextRun(&runs);
+ TObjString *srun = 0x0;
+ Int_t iirun = 1;
+ while ((srun = static_cast<TObjString*>(nextRun()))) {
+ g.GetXaxis()->SetBinLabel(iirun, srun->GetName());
+ iirun++;
+ }
+
}
-TH1D* GraphToHist(Bool_t ToBePlotted,const TGraphAsymmErrors* Graph,Int_t Chamber)
+//---------------------------------------------------------------------------
+void SetRunLabel(TObjArray& array, Int_t irun, const TList& runs)
{
-
- Int_t nBins = 0;
- if ( Chamber < 5 ) nBins = 4;
- if ( (Chamber == 5) || (Chamber == 6) ) nBins = 18;
- if ( (Chamber > 6) && (Chamber < 11) ) nBins = 26;
- if ( Chamber == 11 ) nBins = 10;
- TArrayF bins(nBins+1);
- TArrayF y(nBins);
- TArrayF ey(nBins);
- Double_t dx = 1;
- Double_t xmin = 10000;
- Double_t xmax = -10000;
- for (Int_t i = 0; i < nBins; i++) {
- Double_t x = Graph->GetX()[i];
- Double_t exl = Graph->GetEXlow()[i];
- Double_t exh = Graph->GetEXhigh()[i];
- xmin = TMath::Min(x-exl, xmin);
- xmax = TMath::Max(x+exh, xmax);
- bins.fArray[i] = x-exl;
- bins.fArray[i+1] = x+exh;
- Double_t dxi = exh+exl;
- if (dxi == 0 && i != 0) dxi = bins.fArray[i]-bins.fArray[i-1];
- if (dx == 0) dx = dxi;
- else if (dxi != dx) dx = 0;
-
- y.fArray[i] = Graph->GetY()[i];
- ey.fArray[i] = TMath::Max(Graph->GetEYlow()[i],Graph->GetEYhigh()[i]);
-
- }
- TString name(Graph->GetName());
- TString title(Graph->GetTitle());
- TH1D* h = 0;
- if (ToBePlotted == kTRUE) {
- h = new TH1D(name.Data(), title.Data(), nBins, bins[0]-dx/2 + 1.5, bins[nBins]-dx/2 + 1.5);
- }
- else {
- if ( Chamber < 5 ) h = new TH1D(name.Data(), title.Data(), 4, -0.5, 3.5);
- if ( (Chamber == 5) || (Chamber == 6) ) h = new TH1D(name.Data(), title.Data(), 18, -0.5, 17.5);
- if ( (Chamber > 6) && (Chamber < 11) ) h = new TH1D(name.Data(), title.Data(), 26, -0.5, 25.5);
- if ( Chamber == 11 ) h = new TH1D(name.Data(), title.Data(), 10, -0.5, 9.5);
- }
- for (Int_t i = 1; i <= nBins; i++) {
- h->SetBinContent(i, y.fArray[i-1]);
- h->SetBinError(i, ey.fArray[i-1]);
- }
- h->SetMarkerStyle(Graph->GetMarkerStyle());
- h->SetMarkerColor(Graph->GetMarkerColor());
- h->SetMarkerSize(Graph->GetMarkerSize());
- h->SetDirectory(0);
-
- return h;
+ /// set the run label for all the graphs
+
+ for (Int_t i = 0; i <= array.GetLast(); ++i)
+ SetRunLabel(*static_cast<TGraphAsymmErrors*>(array.UncheckedAt(i)), irun, runs);
+
}
+
+
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff