[u/mrichter/AliRoot.git] / PWG3 / muon / trigEffQA.C

#if !defined(__CINT__) || defined(__MAKECINT__)
// ROOT includes
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TGraphAsymmErrors.h"
#include "TSystem.h"
#include "Riostream.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TROOT.h"
#include "TLine.h"
#include "TLegend.h"
#include "TPaveLabel.h"
#include "TMath.h"
#include "TObjArray.h"
#include "TObjString.h"
#include "TString.h"
#include "TGrid.h"

#include "AliCDBManager.h"
#include "AliCDBEntry.h"
#include "AliCDBPath.h"
#include "AliMUONTriggerEfficiencyCells.h"
#include "AliMUONTriggerChamberEfficiency.h"
#endif

const Int_t kNch = 4;
const Double_t kZero = 1.e-7; // Avoid problems when comparing to 0.
void SetMyStyle();
Int_t GetRunNumber(TString);
//Double_t GetBinomial(Double_t*, Double_t* eff2 = 0x0, Double_t* effBoth = 0x0);
Double_t* GetBinomial(Double_t*, Double_t* effErr2 = 0x0, Double_t* effErrBoth = 0x0);
Double_t* GetProdErr(Double_t*, Int_t, Int_t nFactors=kNch);
Double_t* GetConditionalEffErr(Double_t*, Double_t*, Double_t*, Int_t exclude=-1);
TH1* GetHisto(TString, TFile*, TList*);

void trigEffQA(TString fileListName, TString outFilename = "")
{
  ifstream inFile(fileListName.Data());
  TObjArray fileNameArray, runNumArray;
  fileNameArray.SetOwner();
  runNumArray.SetOwner();
  TString currString = "";
  if (inFile.is_open()) {
    while (! inFile.eof() ) {
      currString.ReadLine(inFile,kTRUE); // Read line
      if ( currString.IsNull() || ! currString.Contains(".root") || 
           currString.BeginsWith("#") ) continue;
      fileNameArray.AddLast(new TObjString(currString.Data()));
      Int_t runNum = GetRunNumber(currString);
      runNumArray.AddLast(new TObjString(Form("%i",runNum)));
    }
    inFile.close();
  }
  else {
    printf("Fatal: cannot open input file\n");
    return;
  }

  runNumArray.Sort();

  //TString elementName[2] = {"Board", "RPC"};
  TString countTypeName[3] = {"bendPlane", "nonBendPlane","bothPlanes"};

  TObjArray chEffList(13);

  TCanvas* can = 0x0;
  TString filename = "", histoName = "", histoTitle = "";

  histoName = "totalEffEvolution";
  histoTitle = "Multinomial probability of firing at least 3/4 chambers";
  TH1* totalEff = new TH1D(histoName.Data(), histoTitle.Data(), 0, 0., 1.);
  chEffList.AddAtAndExpand(totalEff, 12);

  SetMyStyle();
  Double_t effValues[3][2*kNch];

  TString runNumString = "";
  for ( Int_t irun=0; irun<runNumArray.GetEntries(); irun++ ) {
    runNumString = runNumArray.At(irun)->GetName();

    // Search corresponding file (for sorting)
    for ( Int_t ifile=0; ifile<fileNameArray.GetEntries(); ifile++ ) {
      filename = fileNameArray.At(ifile)->GetName();
      if ( filename.Contains(runNumString.Data()) ) break;
    }

    if ( filename.Contains("alien://") && ! gGrid )
      gGrid->Connect("alien://");

    TFile* file = TFile::Open(filename.Data());
    if ( ! file ) {
      printf("Warning: cannot find %s\n", filename.Data());
      continue;
    }
    TList* trigEffList = (TList*)file->FindObjectAny("triggerChamberEff");
    if ( ! trigEffList )
      printf("Warning: histo list not found in %s. Check directly in file\n", filename.Data());

    TH1* histoDen = GetHisto("allTracksCountChamber", file, trigEffList);
    if ( ! histoDen ) {
      printf("Error: cannot find histograms in file %s. Skip to next\n", filename.Data());
      continue;
    }
    for(Int_t icount=0; icount<AliMUONTriggerEfficiencyCells::kNcounts-1; icount++) {
      histoName = countTypeName[icount] + "CountChamber";
      TH1* histoNum = GetHisto(histoName, file, trigEffList);
      TGraphAsymmErrors* graph = new TGraphAsymmErrors(histoNum, histoDen);
      Double_t xpt, ypt;
      for ( Int_t ich=0; ich<kNch; ich++ ) {
        graph->GetPoint(ich, xpt, ypt);
        effValues[icount][ich] = ypt;
        Int_t ihisto = kNch*icount+ich;
        TH1* chEff = (TH1*)chEffList.At(ihisto);
        if ( ! chEff ) {
          histoName = Form("effEvolution%sCh%i", countTypeName[icount].Data(), 11+ich);
          histoTitle = Form("Efficiency of trigger chamber %i", 11+ich);
          chEff = new TH1D(histoName.Data(), histoTitle.Data(), 0, 0., 1.);
          chEffList.AddAtAndExpand(chEff, ihisto);
        }
        Int_t currBin = chEff->Fill(runNumString.Data(), ypt);
        Double_t err = 0.5*(graph->GetErrorYlow(ich) + graph->GetErrorYhigh(ich));
        chEff->SetBinError(currBin, err);
        effValues[icount][ich+kNch] = err;
      } // loop on chambers
    } // loop on counts
    Double_t* binomialEff = GetBinomial(effValues[0], effValues[1], effValues[2]);
    Int_t currBin = totalEff->Fill(runNumString, binomialEff[0]);
    // CAVEAT!!!!
    // Well, error calculation of the binomial efficiency is a mess...
    // Sometimes it happens that efficiency + error > 1.
    // In that case reduce the error.
    totalEff->SetBinError(currBin, TMath::Min(binomialEff[1], 1.-binomialEff[0]));
    delete binomialEff;
  } // loop on runs

  TString outCanName = outFilename;
  outCanName.ReplaceAll(".root",".pdf");
  TString canName = "";
  canName = "totalEff";
  can = new TCanvas(canName.Data(), canName.Data(), 200, 10, 600, 600);
  if ( ! outCanName.IsNull() )
    can->Print(Form("%s[", outCanName.Data())); // Open pdf file

  // Set correct range (do not show last empty bins)
  for ( Int_t ihisto=0; ihisto<chEffList.GetEntries(); ihisto++ ) {
    TH1* histo = (TH1*)chEffList.At(ihisto);
    histo->GetXaxis()->SetRange(1,fileNameArray.GetEntries());
    histo->GetXaxis()->SetLabelSize(0.04);
  }

  TH1* totEff = (TH1*)chEffList.At(12);
  totEff->GetYaxis()->SetRangeUser(0.,1.1);
  totEff->GetYaxis()->SetTitle("Probability to satisfy trigger conditions (3/4)");
  totEff->SetStats(kFALSE);
  totEff->DrawCopy();
  if ( ! outCanName.IsNull() )
    can->Print(outCanName.Data());

  Int_t color[3] = {kBlack, kRed, kBlue};
  Int_t markStyle[3] = {20, 24, 26};
  TLegend* leg = 0x0;

  for ( Int_t ich=0; ich<kNch; ich++ ) {
    canName = Form("trigEffCh%i", 11+ich);
    can = new TCanvas(canName.Data(), canName.Data(), 200, 10, 600, 600);
    can->SetGridy();
    leg = new TLegend(0.6, 0.2, 0.9, 0.4);
    leg->SetBorderSize(1);
    //can->Divide(2,2);
    TString drawOpt = "e";
    for(Int_t icount=0; icount<AliMUONTriggerEfficiencyCells::kNcounts-1; icount++) {
      //can->cd(icount+1);
      TH1* chEff = (TH1*)chEffList.At(4*icount+ich);
      chEff->GetYaxis()->SetRangeUser(0.48,1.05);
      chEff->SetStats(kFALSE);
      chEff->GetYaxis()->SetTitle("Trigger chamber efficiency");
      TH1* copyEff = chEff->DrawCopy(drawOpt.Data());
      copyEff->SetLineColor(color[icount]);
      copyEff->SetMarkerColor(color[icount]);
      copyEff->SetMarkerStyle(markStyle[icount]);
      leg->AddEntry(copyEff, countTypeName[icount].Data(), "lp");
      drawOpt = "esame";
    } // loop on counts
    leg->Draw("same");
    if ( ! outCanName.IsNull() )
      can->Print(outCanName.Data());
  } // loop on chambers
  if ( ! outCanName.IsNull() ) {
    can->Print(Form("%s]", outCanName.Data())); // close pdf file
    // There is a bug when creating a pdf
    // So create a ps and then convert via epstopdf
    if ( outCanName.Contains(".ps") || outCanName.Contains(".eps") ) {
      gSystem->Exec(Form("epstopdf %s", outCanName.Data()));
      gSystem->Exec(Form("rm %s", outCanName.Data()));
    }
  }

  TFile* outFile = 0x0;
  if ( ! outFilename.IsNull() ) {
    outFile = new TFile(outFilename.Data(), "recreate");
    chEffList.Write();
    outFile->Close();
  }
}


//_____________________________________________________________________________
void SetMyStyle()
{
    gStyle->SetCanvasColor(10);
    gStyle->SetFrameFillColor(10);
    gStyle->SetStatColor(10);
    gStyle->SetFillColor(10);
    gStyle->SetTitleFillColor(10);

    gStyle->SetTitleXSize(0.03);
    gStyle->SetTitleXOffset(1.1);
    gStyle->SetTitleYSize(0.03);
    gStyle->SetTitleYOffset(1.9);

    gStyle->SetMarkerSize(0.7);
    gStyle->SetHistLineWidth(2);

    gStyle->SetPadLeftMargin(0.12);
    gStyle->SetPadRightMargin(0.04);
    gStyle->SetPadBottomMargin(0.08);
    gStyle->SetPadTopMargin(0.08);

    gROOT->ForceStyle();    
}

//_____________________________________________________________________________
Int_t GetRunNumber(TString filePath)
{
  TObjArray* array = filePath.Tokenize("/");
  array->SetOwner();
  TString auxString = "";
  Int_t runNum = -1;
  for ( Int_t ientry=0; ientry<array->GetEntries(); ientry++ ) {
    auxString = array->At(ientry)->GetName();
    if ( auxString.Contains("000") ) {
      runNum = auxString.Atoi();
      break;
    }
  }
  delete array;

  if ( runNum < 0 ) {
    array = auxString.Tokenize("_");
    array->SetOwner();
    auxString = array->Last()->GetName();
    auxString.ReplaceAll(".root","");
    runNum = auxString.Atoi();
    delete array;
  }

  return runNum;
}


// //_____________________________________________________________________________
// Double_t GetBinomial(Double_t* eff1, Double_t* eff2, Double_t* effBoth)
// {
//   Double_t binomialEff = 0.;
//   Double_t eff44 = 1.;
//   Double_t auxEff[4];
//   Double_t auxBinomial = 1.;

//   // All chamber efficient
//   for ( Int_t ich=0; ich<4; ich++ ) {
//     eff44 *= eff1[ich];
//   }
//   if ( eff2 ) {
//     for ( Int_t ich=0; ich<4; ich++ ) {
//       auxEff[ich] = ( eff1[ich] > kZero ) ? effBoth[ich]/eff1[ich] : 0.;
//     }
//     auxBinomial = GetBinomial(auxEff);
//     eff44 *= auxBinomial;
//   }

//   binomialEff += eff44;

//   // 1 chamber inefficient
//   for ( Int_t ich=0; ich<4; ich++ ) {
//     Double_t eff34 = 1.;
//     for ( Int_t jch=0; jch<4; jch++ ) {
//       eff34 *= ( ich == jch ) ? ( 1. - eff1[jch] ) : eff1[jch];
//     }
//     if ( eff2 ) {
//       for ( Int_t jch=0; jch<4; jch++ ) {
//         if ( ich == jch ) {
//           auxEff[jch] = ( eff1[jch] < 1. ) ? ( eff2[jch] - effBoth[jch] ) / ( 1. - eff1[jch]) : 0.;
//         }
//         else {
//           auxEff[jch] = ( eff1[ich] > kZero ) ? effBoth[ich]/eff1[ich] : 0.;
//         }
//         auxBinomial = GetBinomial(auxEff);
//         eff34 *= auxBinomial;
//       }
//     }

//     binomialEff += eff34;
//   }

//   return binomialEff;
// }


//_____________________________________________________________________________
Double_t* GetBinomial(Double_t* effErr1, Double_t* effErr2, Double_t* effErrBoth)
{
  Double_t effProd[4];
  Double_t defaultEffErr[2] = {1.,0.};
  Double_t* auxBinomial = 0x0;
  Double_t* currEffErr44 = 0x0;
  Double_t* effErrBinomial = new Double_t[2];
  effErrBinomial[0] = 0.;
  effErrBinomial[1] = 0.;

  for ( Int_t ich = -1; ich<kNch; ich++ ) {
    Double_t* currEffErr = GetProdErr(effErr1, ich);
    if ( ich >= 0 ) {
      currEffErr[0] = currEffErr[0] - currEffErr44[0];
      currEffErr[1] = TMath::Sqrt(currEffErr[1]*currEffErr[1] + currEffErr44[1]*currEffErr44[1]);
    }
    if ( effErr2 ) {
      Double_t* auxEffErr = GetConditionalEffErr(effErr1, effErr2, effErrBoth, ich);
      auxBinomial = GetBinomial(auxEffErr);
      delete auxEffErr;
    }
    for ( Int_t ival=0; ival<2; ival++ ) {
      effProd[2*ival] = currEffErr[ival];
      effProd[2*ival+1] = ( effErr2 ) ? auxBinomial[ival] : defaultEffErr[ival];
    }
    if ( ich < 0 ) currEffErr44 = currEffErr;
    else delete currEffErr;
    delete auxBinomial;

    Double_t* effErr = GetProdErr(effProd, -1, 2);
    //printf("%f * %f = %f\n", effProd[0], effProd[1], effErr[0]); // REMEMBER TO CUT
    effErrBinomial[0] += effErr[0];
    effErrBinomial[1] += effErr[1]*effErr[1];
    delete effErr;
  } // loop on chambers

  delete currEffErr44;

  effErrBinomial[1] = TMath::Sqrt(effErrBinomial[1]);

  return effErrBinomial;
}

//_____________________________________________________________________________
Double_t* GetProdErr(Double_t* effErr, Int_t exclude, Int_t nFactors)
{
  Double_t prod = 1.;
  Double_t relErr = 0., relProdErrSquare = 0.;
  for ( Int_t iprod=0; iprod<nFactors; iprod++ ) {
    if ( iprod == exclude ) continue;
    prod *= effErr[iprod];
    relErr = ( effErr[iprod] > kZero ) ? effErr[iprod+nFactors]/effErr[iprod] : 0.;
    relProdErrSquare += relErr*relErr;
    //printf("%f +- %f  ", effErr[iprod], effErr[iprod+nFactors]); // REMEMBER TO CUT
  }
  Double_t* prodErr = new Double_t[2];
  prodErr[0] = prod;
  prodErr[1] = prod*TMath::Sqrt(relProdErrSquare);
  //printf("-> %f %f\n", prodErr[0], prodErr[1]); // REMEMBER TO CUT
  return prodErr;
}


//_____________________________________________________________________________
Double_t* GetConditionalEffErr(Double_t* effErr1, Double_t* effErr2, Double_t* effErrBoth, Int_t exclude)
{
  Double_t* effErr = new Double_t[2*kNch];
  for ( Int_t ich=0; ich<kNch; ich++ ) {
    if ( ich == exclude ) {
      effErr[ich] = ( effErr1[ich] < 1. ) ? ( effErr2[ich] - effErrBoth[ich] ) / ( 1. - effErr1[ich] ) : 0.;
      effErr[ich+kNch] = 0;
      if ( effErr1[ich] < 1. ) {
        Double_t err2 = effErr2[ich+kNch] / ( 1. - effErr1[ich] );
        Double_t errBoth = effErrBoth[ich+kNch] / ( 1. - effErr1[ich] );
        Double_t err1 = effErr1[ich+kNch] * effErr[ich] / ( 1. - effErr1[ich] );
        effErr[ich+kNch] = TMath::Sqrt(err2*err2 + errBoth*errBoth + err1*err1);
      }
    }
    else {
      effErr[ich] = ( effErr1[ich] > kZero ) ? effErrBoth[ich]/effErr1[ich] : 0.;
      Double_t relErr1 = ( effErr1[ich] > kZero ) ? effErr1[ich+kNch]/effErr1[ich] : 0.;
      Double_t relErrBoth = ( effErrBoth[ich] > kZero ) ? effErrBoth[ich+kNch]/effErrBoth[ich] : 0.;
      effErr[ich+kNch] = effErr[ich] * TMath::Sqrt(relErr1*relErr1 + relErrBoth*relErrBoth);
    }
    //printf("%f  %f  %f -> %f\n", effErr1[ich], effErr2[ich], effErrBoth[ich], effErr[ich]); // REMEMBER TO CUT
  } // loop on chambers
  return effErr;
}


//_____________________________________________________________________________
TH1* GetHisto(TString histoName, TFile* file, TList* histoList)
{
  TH1* histo = 0x0;
  if ( histoList )
    histo = (TH1*)histoList->FindObject(histoName.Data());
  else 
    histo = (TH1*)file->FindObjectAny(histoName.Data());
  
  return histo;
}
Commit	Line	Data
352cd18c	1	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	2	// ROOT includes
	3	#include "TFile.h"
	4	#include "TH1.h"
	5	#include "TH2.h"
	6	#include "TGraphAsymmErrors.h"
	7	#include "TSystem.h"
	8	#include "Riostream.h"
	9	#include "TCanvas.h"
	10	#include "TStyle.h"
	11	#include "TROOT.h"
	12	#include "TLine.h"
	13	#include "TLegend.h"
	14	#include "TPaveLabel.h"
	15	#include "TMath.h"
	16	#include "TObjArray.h"
	17	#include "TObjString.h"
	18	#include "TString.h"
	19	#include "TGrid.h"
	20
	21	#include "AliCDBManager.h"
	22	#include "AliCDBEntry.h"
	23	#include "AliCDBPath.h"
	24	#include "AliMUONTriggerEfficiencyCells.h"
	25	#include "AliMUONTriggerChamberEfficiency.h"
	26	#endif
	27
	28	const Int_t kNch = 4;
8415c156	29	const Double_t kZero = 1.e-7; // Avoid problems when comparing to 0.
352cd18c	30	void SetMyStyle();
	31	Int_t GetRunNumber(TString);
	32	//Double_t GetBinomial(Double_t, Double_t eff2 = 0x0, Double_t* effBoth = 0x0);
	33	Double_t* GetBinomial(Double_t, Double_t effErr2 = 0x0, Double_t* effErrBoth = 0x0);
	34	Double_t* GetProdErr(Double_t*, Int_t, Int_t nFactors=kNch);
	35	Double_t* GetConditionalEffErr(Double_t, Double_t, Double_t*, Int_t exclude=-1);
	36	TH1* GetHisto(TString, TFile, TList);
	37
	38	void trigEffQA(TString fileListName, TString outFilename = "")
	39	{
	40	ifstream inFile(fileListName.Data());
	41	TObjArray fileNameArray, runNumArray;
	42	fileNameArray.SetOwner();
	43	runNumArray.SetOwner();
	44	TString currString = "";
	45	if (inFile.is_open()) {
	46	while (! inFile.eof() ) {
	47	currString.ReadLine(inFile,kTRUE); // Read line
	48	if ( currString.IsNull() \|\| ! currString.Contains(".root") \|\|
	49	currString.BeginsWith("#") ) continue;
	50	fileNameArray.AddLast(new TObjString(currString.Data()));
	51	Int_t runNum = GetRunNumber(currString);
	52	runNumArray.AddLast(new TObjString(Form("%i",runNum)));
	53	}
	54	inFile.close();
	55	}
	56	else {
	57	printf("Fatal: cannot open input file\n");
	58	return;
	59	}
	60
	61	runNumArray.Sort();
	62
	63	//TString elementName[2] = {"Board", "RPC"};
	64	TString countTypeName[3] = {"bendPlane", "nonBendPlane","bothPlanes"};
	65
	66	TObjArray chEffList(13);
	67
	68	TCanvas* can = 0x0;
	69	TString filename = "", histoName = "", histoTitle = "";
	70
	71	histoName = "totalEffEvolution";
	72	histoTitle = "Multinomial probability of firing at least 3/4 chambers";
	73	TH1* totalEff = new TH1D(histoName.Data(), histoTitle.Data(), 0, 0., 1.);
	74	chEffList.AddAtAndExpand(totalEff, 12);
	75
	76	SetMyStyle();
	77	Double_t effValues[3][2*kNch];
	78
	79	TString runNumString = "";
	80	for ( Int_t irun=0; irun<runNumArray.GetEntries(); irun++ ) {
	81	runNumString = runNumArray.At(irun)->GetName();
	82
	83	// Search corresponding file (for sorting)
	84	for ( Int_t ifile=0; ifile<fileNameArray.GetEntries(); ifile++ ) {
	85	filename = fileNameArray.At(ifile)->GetName();
	86	if ( filename.Contains(runNumString.Data()) ) break;
	87	}
	88
	89	if ( filename.Contains("alien://") && ! gGrid )
	90	gGrid->Connect("alien://");
	91
	92	TFile* file = TFile::Open(filename.Data());
	93	if ( ! file ) {
94	printf("Warning: cannot find %s\n", filename.Data());
95	continue;
96	}
97	TList* trigEffList = (TList*)file->FindObjectAny("triggerChamberEff");
98	if ( ! trigEffList )
99	printf("Warning: histo list not found in %s. Check directly in file\n", filename.Data());
100
101	TH1* histoDen = GetHisto("allTracksCountChamber", file, trigEffList);
102	if ( ! histoDen ) {
103	printf("Error: cannot find histograms in file %s. Skip to next\n", filename.Data());
104	continue;
105	}
106	for(Int_t icount=0; icount<AliMUONTriggerEfficiencyCells::kNcounts-1; icount++) {
107	histoName = countTypeName[icount] + "CountChamber";
108	TH1* histoNum = GetHisto(histoName, file, trigEffList);
109	TGraphAsymmErrors* graph = new TGraphAsymmErrors(histoNum, histoDen);
110	Double_t xpt, ypt;
111	for ( Int_t ich=0; ich<kNch; ich++ ) {
112	graph->GetPoint(ich, xpt, ypt);
113	effValues[icount][ich] = ypt;
114	Int_t ihisto = kNch*icount+ich;
115	TH1* chEff = (TH1*)chEffList.At(ihisto);
116	if ( ! chEff ) {
117	histoName = Form("effEvolution%sCh%i", countTypeName[icount].Data(), 11+ich);
118	histoTitle = Form("Efficiency of trigger chamber %i", 11+ich);
119	chEff = new TH1D(histoName.Data(), histoTitle.Data(), 0, 0., 1.);
120	chEffList.AddAtAndExpand(chEff, ihisto);
121	}
122	Int_t currBin = chEff->Fill(runNumString.Data(), ypt);
123	Double_t err = 0.5*(graph->GetErrorYlow(ich) + graph->GetErrorYhigh(ich));
124	chEff->SetBinError(currBin, err);
125	effValues[icount][ich+kNch] = err;
126	} // loop on chambers
127	} // loop on counts
128	Double_t* binomialEff = GetBinomial(effValues[0], effValues[1], effValues[2]);
129	Int_t currBin = totalEff->Fill(runNumString, binomialEff[0]);
130	// CAVEAT!!!!
131	// Well, error calculation of the binomial efficiency is a mess...
132	// Sometimes it happens that efficiency + error > 1.
133	// In that case reduce the error.
134	totalEff->SetBinError(currBin, TMath::Min(binomialEff[1], 1.-binomialEff[0]));
135	delete binomialEff;
136	} // loop on runs
137
138	TString outCanName = outFilename;
139	outCanName.ReplaceAll(".root",".pdf");
140	TString canName = "";
141	canName = "totalEff";
142	can = new TCanvas(canName.Data(), canName.Data(), 200, 10, 600, 600);
143	if ( ! outCanName.IsNull() )
144	can->Print(Form("%s[", outCanName.Data())); // Open pdf file
145
146	// Set correct range (do not show last empty bins)
147	for ( Int_t ihisto=0; ihisto<chEffList.GetEntries(); ihisto++ ) {
148	TH1* histo = (TH1*)chEffList.At(ihisto);
149	histo->GetXaxis()->SetRange(1,fileNameArray.GetEntries());
b8e17107	150	histo->GetXaxis()->SetLabelSize(0.04);
352cd18c	151	}
	152
	153	TH1* totEff = (TH1*)chEffList.At(12);
	154	totEff->GetYaxis()->SetRangeUser(0.,1.1);
	155	totEff->GetYaxis()->SetTitle("Probability to satisfy trigger conditions (3/4)");
	156	totEff->SetStats(kFALSE);
	157	totEff->DrawCopy();
	158	if ( ! outCanName.IsNull() )
	159	can->Print(outCanName.Data());
	160
	161	Int_t color[3] = {kBlack, kRed, kBlue};
	162	Int_t markStyle[3] = {20, 24, 26};
	163	TLegend* leg = 0x0;
	164
	165	for ( Int_t ich=0; ich<kNch; ich++ ) {
	166	canName = Form("trigEffCh%i", 11+ich);
	167	can = new TCanvas(canName.Data(), canName.Data(), 200, 10, 600, 600);
	168	can->SetGridy();
	169	leg = new TLegend(0.6, 0.2, 0.9, 0.4);
	170	leg->SetBorderSize(1);
	171	//can->Divide(2,2);
	172	TString drawOpt = "e";
	173	for(Int_t icount=0; icount<AliMUONTriggerEfficiencyCells::kNcounts-1; icount++) {
	174	//can->cd(icount+1);
	175	TH1* chEff = (TH1)chEffList.At(4icount+ich);
	176	chEff->GetYaxis()->SetRangeUser(0.48,1.05);
	177	chEff->SetStats(kFALSE);
	178	chEff->GetYaxis()->SetTitle("Trigger chamber efficiency");
	179	TH1* copyEff = chEff->DrawCopy(drawOpt.Data());
	180	copyEff->SetLineColor(color[icount]);
	181	copyEff->SetMarkerColor(color[icount]);
	182	copyEff->SetMarkerStyle(markStyle[icount]);
	183	leg->AddEntry(copyEff, countTypeName[icount].Data(), "lp");
	184	drawOpt = "esame";
	185	} // loop on counts
	186	leg->Draw("same");
	187	if ( ! outCanName.IsNull() )
	188	can->Print(outCanName.Data());
	189	} // loop on chambers
	190	if ( ! outCanName.IsNull() ) {
	191	can->Print(Form("%s]", outCanName.Data())); // close pdf file
	192	// There is a bug when creating a pdf
	193	// So create a ps and then convert via epstopdf
	194	if ( outCanName.Contains(".ps") \|\| outCanName.Contains(".eps") ) {
	195	gSystem->Exec(Form("epstopdf %s", outCanName.Data()));
	196	gSystem->Exec(Form("rm %s", outCanName.Data()));
	197	}
	198	}
	199
	200	TFile* outFile = 0x0;
	201	if ( ! outFilename.IsNull() ) {
	202	outFile = new TFile(outFilename.Data(), "recreate");
	203	chEffList.Write();
	204	outFile->Close();
	205	}
	206	}
	207
	208
	209	//_____________________________________________________________________________
	210	void SetMyStyle()
	211	{
	212	gStyle->SetCanvasColor(10);
	213	gStyle->SetFrameFillColor(10);
	214	gStyle->SetStatColor(10);
215	gStyle->SetFillColor(10);
216	gStyle->SetTitleFillColor(10);
217
218	gStyle->SetTitleXSize(0.03);
219	gStyle->SetTitleXOffset(1.1);
220	gStyle->SetTitleYSize(0.03);
221	gStyle->SetTitleYOffset(1.9);
222
223	gStyle->SetMarkerSize(0.7);
224	gStyle->SetHistLineWidth(2);
225
226	gStyle->SetPadLeftMargin(0.12);
227	gStyle->SetPadRightMargin(0.04);
228	gStyle->SetPadBottomMargin(0.08);
229	gStyle->SetPadTopMargin(0.08);
230
231	gROOT->ForceStyle();
232	}
233
234	//_____________________________________________________________________________
235	Int_t GetRunNumber(TString filePath)
236	{
237	TObjArray* array = filePath.Tokenize("/");
238	array->SetOwner();
239	TString auxString = "";
240	Int_t runNum = -1;
241	for ( Int_t ientry=0; ientry<array->GetEntries(); ientry++ ) {
242	auxString = array->At(ientry)->GetName();
243	if ( auxString.Contains("000") ) {
244	runNum = auxString.Atoi();
245	break;
246	}
247	}
248	delete array;
249
250	if ( runNum < 0 ) {
251	array = auxString.Tokenize("_");
252	array->SetOwner();
253	auxString = array->Last()->GetName();
254	auxString.ReplaceAll(".root","");
255	runNum = auxString.Atoi();
256	delete array;
257	}
258
259	return runNum;
260	}
261
262
263	// //_____________________________________________________________________________
264	// Double_t GetBinomial(Double_t* eff1, Double_t* eff2, Double_t* effBoth)
265	// {
266	// Double_t binomialEff = 0.;
267	// Double_t eff44 = 1.;
268	// Double_t auxEff[4];
269	// Double_t auxBinomial = 1.;
270
271	// // All chamber efficient
272	// for ( Int_t ich=0; ich<4; ich++ ) {
273	// eff44 *= eff1[ich];
274	// }
275	// if ( eff2 ) {
276	// for ( Int_t ich=0; ich<4; ich++ ) {
8415c156	277	// auxEff[ich] = ( eff1[ich] > kZero ) ? effBoth[ich]/eff1[ich] : 0.;
352cd18c	278	// }
	279	// auxBinomial = GetBinomial(auxEff);
	280	// eff44 *= auxBinomial;
	281	// }
	282
	283	// binomialEff += eff44;
	284
	285	// // 1 chamber inefficient
	286	// for ( Int_t ich=0; ich<4; ich++ ) {
	287	// Double_t eff34 = 1.;
	288	// for ( Int_t jch=0; jch<4; jch++ ) {
	289	// eff34 *= ( ich == jch ) ? ( 1. - eff1[jch] ) : eff1[jch];
	290	// }
	291	// if ( eff2 ) {
	292	// for ( Int_t jch=0; jch<4; jch++ ) {
	293	// if ( ich == jch ) {
	294	// auxEff[jch] = ( eff1[jch] < 1. ) ? ( eff2[jch] - effBoth[jch] ) / ( 1. - eff1[jch]) : 0.;
	295	// }
	296	// else {
8415c156	297	// auxEff[jch] = ( eff1[ich] > kZero ) ? effBoth[ich]/eff1[ich] : 0.;
352cd18c	298	// }
	299	// auxBinomial = GetBinomial(auxEff);
	300	// eff34 *= auxBinomial;
	301	// }
	302	// }
	303
	304	// binomialEff += eff34;
	305	// }
	306
	307	// return binomialEff;
	308	// }
	309
	310
	311	//_____________________________________________________________________________
	312	Double_t* GetBinomial(Double_t* effErr1, Double_t* effErr2, Double_t* effErrBoth)
	313	{
	314	Double_t effProd[4];
	315	Double_t defaultEffErr[2] = {1.,0.};
	316	Double_t* auxBinomial = 0x0;
	317	Double_t* currEffErr44 = 0x0;
	318	Double_t* effErrBinomial = new Double_t[2];
	319	effErrBinomial[0] = 0.;
	320	effErrBinomial[1] = 0.;
	321
	322	for ( Int_t ich = -1; ich<kNch; ich++ ) {
	323	Double_t* currEffErr = GetProdErr(effErr1, ich);
	324	if ( ich >= 0 ) {
	325	currEffErr[0] = currEffErr[0] - currEffErr44[0];
	326	currEffErr[1] = TMath::Sqrt(currEffErr[1]currEffErr[1] + currEffErr44[1]currEffErr44[1]);
	327	}
	328	if ( effErr2 ) {
	329	Double_t* auxEffErr = GetConditionalEffErr(effErr1, effErr2, effErrBoth, ich);
	330	auxBinomial = GetBinomial(auxEffErr);
	331	delete auxEffErr;
	332	}
	333	for ( Int_t ival=0; ival<2; ival++ ) {
	334	effProd[2*ival] = currEffErr[ival];
	335	effProd[2*ival+1] = ( effErr2 ) ? auxBinomial[ival] : defaultEffErr[ival];
	336	}
	337	if ( ich < 0 ) currEffErr44 = currEffErr;
	338	else delete currEffErr;
	339	delete auxBinomial;
	340
	341	Double_t* effErr = GetProdErr(effProd, -1, 2);
	342	//printf("%f * %f = %f\n", effProd[0], effProd[1], effErr[0]); // REMEMBER TO CUT
	343	effErrBinomial[0] += effErr[0];
	344	effErrBinomial[1] += effErr[1]*effErr[1];
	345	delete effErr;
	346	} // loop on chambers
	347
	348	delete currEffErr44;
	349
	350	effErrBinomial[1] = TMath::Sqrt(effErrBinomial[1]);
	351
	352	return effErrBinomial;
	353	}
	354
	355	//_____________________________________________________________________________
	356	Double_t* GetProdErr(Double_t* effErr, Int_t exclude, Int_t nFactors)
	357	{
	358	Double_t prod = 1.;
	359	Double_t relErr = 0., relProdErrSquare = 0.;
	360	for ( Int_t iprod=0; iprod<nFactors; iprod++ ) {
	361	if ( iprod == exclude ) continue;
362	prod *= effErr[iprod];
8415c156	363	relErr = ( effErr[iprod] > kZero ) ? effErr[iprod+nFactors]/effErr[iprod] : 0.;
352cd18c	364	relProdErrSquare += relErr*relErr;
	365	//printf("%f +- %f ", effErr[iprod], effErr[iprod+nFactors]); // REMEMBER TO CUT
	366	}
	367	Double_t* prodErr = new Double_t[2];
	368	prodErr[0] = prod;
	369	prodErr[1] = prod*TMath::Sqrt(relProdErrSquare);
	370	//printf("-> %f %f\n", prodErr[0], prodErr[1]); // REMEMBER TO CUT
	371	return prodErr;
	372	}
	373
	374
	375	//_____________________________________________________________________________
	376	Double_t* GetConditionalEffErr(Double_t* effErr1, Double_t* effErr2, Double_t* effErrBoth, Int_t exclude)
	377	{
	378	Double_t* effErr = new Double_t[2*kNch];
	379	for ( Int_t ich=0; ich<kNch; ich++ ) {
	380	if ( ich == exclude ) {
	381	effErr[ich] = ( effErr1[ich] < 1. ) ? ( effErr2[ich] - effErrBoth[ich] ) / ( 1. - effErr1[ich] ) : 0.;
	382	effErr[ich+kNch] = 0;
	383	if ( effErr1[ich] < 1. ) {
	384	Double_t err2 = effErr2[ich+kNch] / ( 1. - effErr1[ich] );
	385	Double_t errBoth = effErrBoth[ich+kNch] / ( 1. - effErr1[ich] );
	386	Double_t err1 = effErr1[ich+kNch] * effErr[ich] / ( 1. - effErr1[ich] );
	387	effErr[ich+kNch] = TMath::Sqrt(err2err2 + errBotherrBoth + err1*err1);
	388	}
	389	}
	390	else {
8415c156	391	effErr[ich] = ( effErr1[ich] > kZero ) ? effErrBoth[ich]/effErr1[ich] : 0.;
	392	Double_t relErr1 = ( effErr1[ich] > kZero ) ? effErr1[ich+kNch]/effErr1[ich] : 0.;
	393	Double_t relErrBoth = ( effErrBoth[ich] > kZero ) ? effErrBoth[ich+kNch]/effErrBoth[ich] : 0.;
352cd18c	394	effErr[ich+kNch] = effErr[ich] * TMath::Sqrt(relErr1relErr1 + relErrBothrelErrBoth);
	395	}
	396	//printf("%f %f %f -> %f\n", effErr1[ich], effErr2[ich], effErrBoth[ich], effErr[ich]); // REMEMBER TO CUT
	397	} // loop on chambers
	398	return effErr;
	399	}
	400
	401
	402	//_____________________________________________________________________________
	403	TH1* GetHisto(TString histoName, TFile* file, TList* histoList)
	404	{
	405	TH1* histo = 0x0;
	406	if ( histoList )
	407	histo = (TH1*)histoList->FindObject(histoName.Data());
	408	else
	409	histo = (TH1*)file->FindObjectAny(histoName.Data());
	410
	411	return histo;
	412	}