decrease output size for PID task (Alex W)
authorabercuci <abercuci@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 10 Oct 2008 08:06:14 +0000 (08:06 +0000)
committerabercuci <abercuci@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 10 Oct 2008 08:06:14 +0000 (08:06 +0000)
TRD/AliTRDpidUtil.cxx
TRD/AliTRDpidUtil.h
TRD/qaRec/AliTRDpidChecker.cxx
TRD/qaRec/AliTRDpidChecker.h

index 9eae493..15f6103 100644 (file)
@@ -30,7 +30,7 @@ AliTRDpidUtil::AliTRDpidUtil()
 
 
 //________________________________________________________________________
-void  AliTRDpidUtil::CalculatePionEffi(TH1F* histo1, TH1F* histo2)
+void  AliTRDpidUtil::CalculatePionEffi(TH1* histo1, TH1* histo2)
 // Double_t  AliTRDpidUtil::GetError()
 {
   //
@@ -50,17 +50,18 @@ void  AliTRDpidUtil::CalculatePionEffi(TH1F* histo1, TH1F* histo2)
   Int_t abinE, bbinE, cbinE = -1;                    
   Double_t aBinSumE, bBinSumE;                  // content of a single bin
   Bool_t bFirst = 1;                            // checks if threshold is crossed for the first time
-  Double_t SumElecsE[kBins], SumPionsE[kBins];  // array of the integrated sum in each bin
-  memset(SumElecsE, 0, kBins*sizeof(Double_t));
-  memset(SumPionsE, 0, kBins*sizeof(Double_t));
+  Double_t SumElecsE[kBins+2], SumPionsE[kBins+2];  // array of the integrated sum in each bin
+  memset(SumElecsE, 0, (kBins+2)*sizeof(Double_t));
+  memset(SumPionsE, 0, (kBins+2)*sizeof(Double_t));
 
 
   // calculate electron efficiency of each bin
-  for (abinE = (histo1 -> GetNbinsX())-2; abinE >= 0; abinE--){  
-    aBinSumE = 0;
+  for (abinE = (histo1 -> GetNbinsX()); abinE >= 0; abinE--){  
+   aBinSumE = 0;
     aBinSumE = histo1 -> GetBinContent(abinE);
     
     SumElecsE[abinE] = SumElecsE[abinE+1] + aBinSumE;
+
     if((SumElecsE[abinE] >= fEleEffi) && (bFirst == 1)){
       bFirst = 0;
       cbinE = abinE;
@@ -71,7 +72,7 @@ void  AliTRDpidUtil::CalculatePionEffi(TH1F* histo1, TH1F* histo2)
   fThreshold = histo1 -> GetBinCenter(cbinE);
 
   // calculate pion efficiency of each bin
-  for (bbinE = (histo2 -> GetNbinsX())-2; bbinE >= abinE; bbinE--){    
+  for (bbinE = (histo2 -> GetNbinsX()); bbinE >= abinE; bbinE--){      
     bBinSumE = 0;
     bBinSumE = histo2 -> GetBinContent(bbinE);
 
@@ -88,14 +89,17 @@ void  AliTRDpidUtil::CalculatePionEffi(TH1F* histo1, TH1F* histo2)
   // use fit function to get derivate of the TGraph for error calculation
   TF1 *f1 = new TF1("f1","[0]*x*x+[1]*x+[2]", fEleEffi-.05, fEleEffi+.05);
   gEffis -> Fit("f1","Q","",fEleEffi-.05, fEleEffi+.05);
-  AliInfo(Form("Derivative at %4.2f : %f", fEleEffi, f1 -> Derivative(fEleEffi)));
-
+  
   // return the error of the pion efficiency
   if(((1.-fPionEffi) < 0) || ((1.-fCalcEleEffi) < 0)){
     AliWarning(" EleEffi or PioEffi > 1. Error can not be calculated. Please increase statistics or check your simulation!");
     return;
   }
   fError = sqrt(((1/histo2 -> GetEntries())*fPionEffi*(1-fPionEffi))+((f1 -> Derivative(fEleEffi))*(f1 -> Derivative(fEleEffi))*(1/histo1 -> GetEntries())*fCalcEleEffi*(1-fCalcEleEffi)));
+
+  AliInfo(Form("Pion Effi at [%f] : [%f +/- %f], Threshold[%f]", fCalcEleEffi, fPionEffi, fError, fThreshold));
+  AliInfo(Form("Derivative at %4.2f : %f\n", fEleEffi, f1 -> Derivative(fEleEffi)));
+
 }
 
 
index 34a6fa4..7abab04 100644 (file)
@@ -9,7 +9,7 @@
 //
 ///////////////////////////////////////////////////////
 
-class TH1F;
+class TH1;
 class AliTRDpidUtil : public TObject {
 public:
 
@@ -20,7 +20,7 @@ public:
   AliTRDpidUtil();
   virtual ~AliTRDpidUtil(){;}
 
-  void     CalculatePionEffi(TH1F* histo1, TH1F* histo2);
+  void     CalculatePionEffi(TH1* histo1, TH1* histo2);
 
   static Float_t ElectronEfficiency() { return fEleEffi;};
 
index a29c607..ac4e2a5 100644 (file)
@@ -1,8 +1,10 @@
 #include "TPDGCode.h"
 #include "TF1.h"
 #include "TH1F.h"
+#include "TH1D.h"
 #include "TH2F.h"
 #include "TProfile.h"
+#include "TProfile2D.h"
 #include "TGraph.h"
 #include "TGraphErrors.h"
 
@@ -24,6 +26,7 @@
 #include "AliTRDReconstructor.h"
 #include "AliCDBManager.h"
 // #include "../Cal/AliTRDCalPID.h"
+#include "AliTRDpidUtil.h"
 
 
 #include "AliTRDpidChecker.h"
@@ -62,40 +65,41 @@ void AliTRDpidChecker::CreateOutputObjects()
   // Called once
 
   OpenFile(0, "RECREATE");
+  Int_t xBins = AliPID::kSPECIES*AliTRDCalPID::kNMom; 
   fContainer = new TObjArray();
-  fContainer -> Expand(kGraphNNerr + 1);
+  fContainer -> Expand(kGraphNN + 1);
 
-  const Float_t epsilon = 1/(2*(kBins-1));     // get nice histos with bin center at 0 and 1
+  const Float_t epsilon = 1./(2*(AliTRDpidUtil::kBins-1));     // get nice histos with bin center at 0 and 1
 
   // histos of the electron probability of all 5 particle species and 11 momenta for the 2-dim LQ method 
-  for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
-    for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
-      fContainer->AddAt(new TH1F(Form("PID%d_%d_LQ",iPart,iMom), Form("PID distribution for %d_%d", iPart, iMom), kBins, 0.-epsilon, 1.+epsilon), kLQlikelihood + iPart*AliTRDCalPID::kNMom + iMom);
-      if(fDebugLevel>=4) Printf("LQ Particle[%d] Momentum[%d] : [%d]", iPart, iMom, kLQlikelihood + iPart*AliTRDCalPID::kNMom + iMom);
-    }
-  }
+  fContainer->AddAt(
+    new TH2F("PID_LQ", "", 
+      xBins, -0.5, xBins - 0.5,
+      AliTRDpidUtil::kBins, 0.-epsilon, 1.+epsilon)
+  ,kLQlikelihood);
+
 
   // histos of the electron probability of all 5 particle species and 11 momenta for the neural network method
-  for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
-    for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
-      fContainer->AddAt(new TH1F(Form("PID%d_%d_NN",iPart,iMom), Form("PID distribution for %d_%d", iPart, iMom), kBins, 0.-epsilon, 1.+epsilon), kNNlikelihood + iPart*AliTRDCalPID::kNMom + iMom);
-      if(fDebugLevel>=4) Printf("NN Particle[%d] Momentum[%d] : [%d]", iPart, iMom, kNNlikelihood + iPart*AliTRDCalPID::kNMom + iMom);
-    }
-  }
+  fContainer->AddAt(
+    new TH2F("PID_NN", "", 
+      xBins, -0.5, xBins - 0.5,
+      AliTRDpidUtil::kBins, 0.-epsilon, 1.+epsilon)
+  ,kNNlikelihood);
 
   // histos of the dE/dx distribution for all 5 particle species and 11 momenta 
-  for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
-    for(Int_t iMom = 0;  iMom < AliTRDCalPID::kNMom; iMom++){
-      fContainer->AddAt(new TH1F(Form("dEdx%d_%d",iPart,iMom), Form("dEdx distribution for %d_%d", iPart, iMom), 200, 0, 10000), kdEdx + iPart*AliTRDCalPID::kNMom + iMom);
-    }
-  }
+  fContainer->AddAt(
+    new TH2F("dEdx", "", 
+      xBins, -0.5, xBins - 0.5,
+      200, 0, 10000)
+    ,kdEdx);
 
   // histos of the pulse height distribution for all 5 particle species and 11 momenta 
-  for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
-    for(Int_t iMom = 0;  iMom < AliTRDCalPID::kNMom; iMom++){
-      fContainer->AddAt(new TProfile(Form("PH%d_%d",iPart,iMom), Form("PH distribution for %d_%d", iPart, iMom), AliTRDtrackerV1::GetNTimeBins(), -0.5, AliTRDtrackerV1::GetNTimeBins() - 0.5), kPH + iPart*AliTRDCalPID::kNMom + iMom);
-    }
-  }
+  fContainer->AddAt(
+    new TProfile2D("PH", "", 
+      xBins, -0.5, xBins - 0.5,
+      AliTRDtrackerV1::GetNTimeBins(), -0.5, AliTRDtrackerV1::GetNTimeBins() - 0.5)
+    ,kPH);
+
 
   // momentum distributions - absolute and in momentum bins
   fContainer->AddAt(new TH1F("hMom", "momentum distribution", 100, 0., 12.),kMomentum);
@@ -104,31 +108,22 @@ void AliTRDpidChecker::CreateOutputObjects()
 
   // TGraph of the pion efficiencies
 
-  TGraph *gEffisLQ = new TGraph(AliTRDCalPID::kNMom);
-  gEffisLQ -> SetLineColor(4);
-  gEffisLQ -> SetMarkerColor(4);
-  gEffisLQ -> SetMarkerStyle(29);
-
-  TGraphErrors *gEffisErrLQ = new TGraphErrors(AliTRDCalPID::kNMom);
-  gEffisErrLQ -> SetLineColor(4);
-  gEffisErrLQ -> SetMarkerColor(4);
-  gEffisErrLQ -> SetMarkerStyle(29);
+  TGraphErrors *gEffisLQ = 0x0;
+  TGraphErrors *gEffisNN = 0x0;
 
-  TGraph *gEffisNN = new TGraph(AliTRDCalPID::kNMom);
-  gEffisNN -> SetLineColor(2);
-  gEffisNN -> SetMarkerColor(2);
-  gEffisNN -> SetMarkerStyle(29);
+  fContainer->AddAt(gEffisLQ = new TGraphErrors(), kGraphLQ);
+  gEffisLQ->SetLineColor(kBlue);
+  gEffisLQ->SetMarkerColor(kBlue);
+  gEffisLQ->SetMarkerStyle(29);
 
-  TGraphErrors *gEffisErrNN = new TGraphErrors(AliTRDCalPID::kNMom);
-  gEffisErrNN -> SetLineColor(2);
-  gEffisErrNN -> SetMarkerColor(2);
-  gEffisErrNN -> SetMarkerStyle(29);
+  fContainer -> AddAt(gEffisNN = new TGraphErrors(),kGraphNN);
+  gEffisNN->SetLineColor(kRed);
+  gEffisNN->SetMarkerColor(kRed);
+  gEffisNN->SetMarkerStyle(29);
 
+  gEffisLQ -> SetNameTitle("gEffisLQErr", "Efficiencies and Errors of the 2-dim LQ method");
+  gEffisNN -> SetNameTitle("gEffisNNErr", "Efficiencies and Errors of the NN method");
 
-  fContainer -> AddAt(gEffisLQ,kGraphLQ);
-  fContainer -> AddAt(gEffisErrLQ,kGraphLQerr);
-  fContainer -> AddAt(gEffisNN,kGraphNN);
-  fContainer -> AddAt(gEffisErrNN,kGraphNNerr);  
 }
 
 //________________________________________________________________________
@@ -143,20 +138,16 @@ void AliTRDpidChecker::Exec(Option_t *)
 //     AliTracker::SetFieldMap(field, kTRUE);
 //   }
 
-  TH1F *hPIDLQ[AliPID::kSPECIES][AliTRDCalPID::kNMom];
-  TH1F *hPIDNN[AliPID::kSPECIES][AliTRDCalPID::kNMom];
-  TH1F *hdEdx[AliPID::kSPECIES][AliTRDCalPID::kNMom];
-  TProfile *hPH[AliPID::kSPECIES][AliTRDCalPID::kNMom];
-
-  for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
-    for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
-      hPIDLQ[iPart][iMom] = (TH1F*)fContainer->At(kLQlikelihood + iPart*AliTRDCalPID::kNMom+iMom);
-      hPIDNN[iPart][iMom] = (TH1F*)fContainer->At(kNNlikelihood + iPart*AliTRDCalPID::kNMom+iMom);
-      hdEdx[iPart][iMom]  = (TH1F*)fContainer->At(kdEdx + iPart*AliTRDCalPID::kNMom+iMom);
-      hPH[iPart][iMom]    = (TProfile*)fContainer->At(kPH + iPart*AliTRDCalPID::kNMom+iMom);
-    }
-  }
-       
+  TH2F *hPIDLQ;
+  TH2F *hPIDNN;
+  TH2F *hdEdx;
+  TProfile2D *hPH;
+
+  hPIDLQ = (TH2F*)fContainer->At(kLQlikelihood);
+  hPIDNN = (TH2F*)fContainer->At(kNNlikelihood);
+  hdEdx  = (TH2F*)fContainer->At(kdEdx);
+  hPH    = (TProfile2D*)fContainer->At(kPH);
+  
   TH1F *hMom    = (TH1F*)fContainer->At(kMomentum);    
   TH1F *hMomBin = (TH1F*)fContainer->At(kMomentumBin); 
   
@@ -179,6 +170,7 @@ void AliTRDpidChecker::Exec(Option_t *)
 
   AliTRDcluster *TRDcluster = 0x0;
 
+  AliTRDpidUtil *util = new AliTRDpidUtil();
   for(Int_t itrk=0; itrk<fTracks->GetEntriesFast(); itrk++){
     track = (AliTRDtrackInfo*)fTracks->UncheckedAt(itrk);
     if(!track->HasESDtrack()) continue;
@@ -201,17 +193,9 @@ void AliTRDpidChecker::Exec(Option_t *)
 
     // set the 11 momentum bins
     Int_t iMomBin = -1;
-    if(mom < .7) iMomBin = 0;
-    else if(mom < .9) iMomBin = 1;
-    else if(mom < 1.25) iMomBin = 2;
-    else if(mom < 1.75) iMomBin = 3;
-    else if(mom < 2.5) iMomBin = 4;
-    else if(mom < 3.5) iMomBin = 5;
-    else if(mom < 4.5) iMomBin = 6;
-    else if(mom < 5.5) iMomBin = 7;
-    else if(mom < 7.0) iMomBin = 8;
-    else if(mom < 9.0) iMomBin = 9;
-    else  iMomBin = 10;
+    iMomBin = util -> GetMomentumBin(mom);
+    if(fDebugLevel>=4) Printf("MomBin[%d] MomTot[%f]", iMomBin, mom);
+
 
     // fill momentum histo to have the momentum distribution
     hMom -> Fill(mom);
@@ -262,61 +246,61 @@ void AliTRDpidChecker::Exec(Option_t *)
     switch(track->GetPDG()){
     case kElectron:
     case kPositron:
-      hPIDLQ[AliPID::kElectron][iMomBin] -> Fill(TRDtrack -> GetPID(0));
+      hPIDLQ -> Fill(AliPID::kElectron * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
       for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
-        hdEdx[AliPID::kElectron][iMomBin] -> Fill(SumdEdx[iChamb]);
+        hdEdx -> Fill(AliPID::kElectron * AliTRDCalPID::kNMom + iMomBin, SumdEdx[iChamb]);
         for(Int_t iClus = 0; iClus < AliTRDtrackerV1::GetNTimeBins(); iClus++){
           if(!(TRDcluster = (AliTRDcluster*)TRDtracklet[iChamb] -> GetClusters(iClus)))
             continue;
-          hPH[AliPID::kElectron][iMomBin] -> Fill(TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus));
+          hPH -> Fill(AliPID::kElectron * AliTRDCalPID::kNMom + iMomBin, TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus));
         }
       }
       break;
     case kMuonPlus:
     case kMuonMinus:
-      hPIDLQ[AliPID::kMuon][iMomBin] -> Fill(TRDtrack -> GetPID(0));
+      hPIDLQ -> Fill(AliPID::kMuon * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
       for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
-        hdEdx[AliPID::kMuon][iMomBin] -> Fill(SumdEdx[iChamb]);
+        hdEdx -> Fill(AliPID::kMuon * AliTRDCalPID::kNMom + iMomBin, SumdEdx[iChamb]);
         for(Int_t iClus = 0; iClus < AliTRDtrackerV1::GetNTimeBins(); iClus++){
           if(!(TRDcluster = (AliTRDcluster*)TRDtracklet[iChamb] -> GetClusters(iClus)))
             continue;
-          hPH[AliPID::kMuon][iMomBin] -> Fill(TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus));
+          hPH -> Fill(AliPID::kMuon * AliTRDCalPID::kNMom + iMomBin, TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus));
         }
       }
       break;
     case kPiPlus:
     case kPiMinus:
-      hPIDLQ[AliPID::kPion][iMomBin] -> Fill(TRDtrack -> GetPID(0));
+      hPIDLQ -> Fill(AliPID::kPion * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
       for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
-        hdEdx[AliPID::kPion][iMomBin] -> Fill(SumdEdx[iChamb]);
+        hdEdx -> Fill(AliPID::kPion * AliTRDCalPID::kNMom + iMomBin, SumdEdx[iChamb]);
         for(Int_t iClus = 0; iClus < AliTRDtrackerV1::GetNTimeBins(); iClus++){
           if(!(TRDcluster = (AliTRDcluster*)TRDtracklet[iChamb] -> GetClusters(iClus)))
             continue;
-          hPH[AliPID::kPion][iMomBin] -> Fill(TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus));
+          hPH -> Fill(AliPID::kPion * AliTRDCalPID::kNMom + iMomBin, TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus));
         }
       }
       break;
     case kKPlus:
     case kKMinus:
-      hPIDLQ[AliPID::kKaon][iMomBin] -> Fill(TRDtrack -> GetPID(0));
+      hPIDLQ -> Fill(AliPID::kKaon * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
       for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
-        hdEdx[AliPID::kKaon][iMomBin] -> Fill(SumdEdx[iChamb]);
+        hdEdx -> Fill(AliPID::kKaon * AliTRDCalPID::kNMom + iMomBin, SumdEdx[iChamb]);
         for(Int_t iClus = 0; iClus < AliTRDtrackerV1::GetNTimeBins(); iClus++){
           if(!(TRDcluster = (AliTRDcluster*)TRDtracklet[iChamb] -> GetClusters(iClus)))
             continue;
-          hPH[AliPID::kKaon][iMomBin] -> Fill(TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus));
+          hPH -> Fill(AliPID::kKaon * AliTRDCalPID::kNMom + iMomBin, TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus));
         }
       }
       break;
     case kProton:
     case kProtonBar:
-      hPIDLQ[AliPID::kProton][iMomBin] -> Fill(TRDtrack -> GetPID(0));
+      hPIDLQ -> Fill(AliPID::kProton * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
       for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
-        hdEdx[AliPID::kProton][iMomBin] -> Fill(SumdEdx[iChamb]);
+        hdEdx -> Fill(AliPID::kProton * AliTRDCalPID::kNMom + iMomBin, SumdEdx[iChamb]);
         for(Int_t iClus = 0; iClus < AliTRDtrackerV1::GetNTimeBins(); iClus++){
           if(!(TRDcluster = (AliTRDcluster*)TRDtracklet[iChamb] -> GetClusters(iClus)))
             continue;
-          hPH[AliPID::kProton][iMomBin] -> Fill(TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus));
+          hPH -> Fill(AliPID::kProton * AliTRDCalPID::kNMom + iMomBin, TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus));
         }
       }
       break;
@@ -334,42 +318,42 @@ void AliTRDpidChecker::Exec(Option_t *)
     switch(track->GetPDG()){
     case kElectron:
     case kPositron:
-      hPIDNN[AliPID::kElectron][iMomBin] -> Fill(TRDtrack -> GetPID(0));
+      hPIDNN -> Fill(AliPID::kElectron * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
       break;
     case kMuonPlus:
     case kMuonMinus:
-      hPIDNN[AliPID::kMuon][iMomBin] -> Fill(TRDtrack -> GetPID(0));
+      hPIDNN -> Fill(AliPID::kMuon * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
       break;
     case kPiPlus:
     case kPiMinus:
-      hPIDNN[AliPID::kPion][iMomBin] -> Fill(TRDtrack -> GetPID(0));
+      hPIDNN -> Fill(AliPID::kPion * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
       break;
     case kKPlus:
     case kKMinus:
-      hPIDNN[AliPID::kKaon][iMomBin] -> Fill(TRDtrack -> GetPID(0));
+      hPIDNN -> Fill(AliPID::kKaon * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
       break;
     case kProton:
     case kProtonBar:
-      hPIDNN[AliPID::kProton][iMomBin] -> Fill(TRDtrack -> GetPID(0));
+      hPIDNN -> Fill(AliPID::kProton * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
       break;
     }
+
   }
 
+  util -> Delete();
   PostData(0, fContainer);
 }
 
 //________________________________________________________
-void AliTRDpidChecker::GetRefFigure(Int_t ifig, Int_t &first, Int_t &last, Option_t */*opt*/)
+void AliTRDpidChecker::GetRefFigure(Int_t ifig, Int_t &first, Int_t &last, Option_t *opt)
 {
+  opt = "pl";
   switch(ifig){
   case 0:
-    first = (Int_t)kGraphStart+1; last = first+1;
-    break;
-  case 1:
-    first = (Int_t)kGraphStart+3; last = first+1;
+    first = (Int_t)kGraphStart; last = first+2;
     break;
   default:
-    first = (Int_t)kGraphStart; last = first;
+    first = (Int_t)kGraphStart; last = first+2;
     break;
   }
 }
@@ -387,65 +371,88 @@ Bool_t AliTRDpidChecker::PostProcess()
   }
 //   return kTRUE; // testing protection
 
-  // normalize the dE/dx histos
-  const Int_t kNSpectra = AliPID::kSPECIES*AliTRDCalPID::kNMom; 
-  TH1F *hdEdx[kNSpectra];
-  for(Int_t iHist = 0; iHist < kNSpectra; iHist++){
-    hdEdx[iHist] = (TH1F*)fContainer->At(kdEdx + iHist);
-    if(hdEdx[iHist] -> GetEntries() > 0)
-      hdEdx[iHist] -> Scale(1./hdEdx[iHist] -> GetEntries());
-    else continue;
-  }
-
 
   // container for the pion efficiencies and the errors
   Double_t  PionEffiLQ[AliTRDCalPID::kNMom], 
-            PionEffiNN[AliTRDCalPID::kNMom],
             PionEffiErrorLQ[AliTRDCalPID::kNMom], 
-            PionEffiErrorNN[AliTRDCalPID::kNMom];
+            EleEffiLQ[AliTRDCalPID::kNMom],
+            ThresholdLQ[AliTRDCalPID::kNMom];
+
+  Double_t  PionEffiNN[AliTRDCalPID::kNMom],
+            PionEffiErrorNN[AliTRDCalPID::kNMom],
+            EleEffiNN[AliTRDCalPID::kNMom],
+            ThresholdNN[AliTRDCalPID::kNMom];
 
+  Float_t mom = 0.;
+
+  TH1D *Histo1=0x0, *Histo2=0x0;
+
+  TH2F *hPIDLQ=0x0, *hPIDNN=0x0;
+  hPIDLQ = (TH2F*)fContainer->At(kLQlikelihood);
+  hPIDNN = (TH2F*)fContainer->At(kNNlikelihood);
 
   // calculate the pion efficiencies and the errors for 90% electron efficiency (2-dim LQ)
   for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
-    PionEffiLQ[iMom] = GetPionEfficiency((kLQlikelihood + iMom) + (AliPID::kElectron * AliTRDCalPID::kNMom), (kLQlikelihood + iMom) + (AliPID::kPion * AliTRDCalPID::kNMom));
-    PionEffiErrorLQ[iMom] = GetError((kLQlikelihood + iMom) + (AliPID::kElectron * AliTRDCalPID::kNMom), (kLQlikelihood + iMom) + (AliPID::kPion * AliTRDCalPID::kNMom));
+
+    AliTRDpidUtil *util = new AliTRDpidUtil();
+    mom = AliTRDCalPID::GetMomentum(iMom);
+
+    Histo1 = hPIDLQ -> ProjectionY("LQ_ele",AliTRDCalPID::kNMom*AliPID::kElectron+iMom+1,AliTRDCalPID::kNMom*AliPID::kElectron+iMom+1);
+    Histo2 = hPIDLQ -> ProjectionY("LQ_pio",AliTRDCalPID::kNMom*AliPID::kPion+iMom+1,AliTRDCalPID::kNMom*AliPID::kPion+iMom+1);
+
+    util -> CalculatePionEffi(Histo1, Histo2);
+
+    PionEffiLQ[iMom] = util -> GetPionEfficiency();
+    PionEffiErrorLQ[iMom] = util -> GetError();
+    EleEffiLQ[iMom] = util -> GetCalcElectronEfficiency();
+    ThresholdLQ[iMom] = util -> GetThreshold();
+
     if(fDebugLevel>=1) Printf("Pion Efficiency for 2-dim LQ is : %f +/- %f\n\n", PionEffiLQ[iMom], PionEffiErrorLQ[iMom]);
+    util -> Delete();
   }
   
+
+
   // calculate the pion efficiencies and the errors for 90% electron efficiency (NN)
   for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
-    PionEffiNN[iMom] = GetPionEfficiency((kNNlikelihood + iMom) + (AliPID::kElectron * AliTRDCalPID::kNMom), (kNNlikelihood + iMom) + (AliPID::kPion * AliTRDCalPID::kNMom));
-    PionEffiErrorNN[iMom] = GetError((kNNlikelihood + iMom) + (AliPID::kElectron * AliTRDCalPID::kNMom), (kNNlikelihood + iMom) + (AliPID::kPion * AliTRDCalPID::kNMom));
+
+    AliTRDpidUtil *util = new AliTRDpidUtil();
+    mom = AliTRDCalPID::GetMomentum(iMom);
+
+    Histo1 = hPIDNN -> ProjectionY("NN_ele",AliTRDCalPID::kNMom*AliPID::kElectron+iMom+1,AliTRDCalPID::kNMom*AliPID::kElectron+iMom+1);
+    Histo2 = hPIDNN -> ProjectionY("NN_pio",AliTRDCalPID::kNMom*AliPID::kPion+iMom+1,AliTRDCalPID::kNMom*AliPID::kPion+iMom+1);
+
+    util -> CalculatePionEffi(Histo1, Histo2);
+
+    PionEffiNN[iMom] = util -> GetPionEfficiency();
+    PionEffiErrorNN[iMom] = util -> GetError();
+    EleEffiNN[iMom] = util -> GetCalcElectronEfficiency();
+    ThresholdNN[iMom] = util -> GetThreshold();
+
     if(fDebugLevel>=1) Printf("Pion Efficiency for NN is : %f +/- %f\n\n", PionEffiNN[iMom], PionEffiErrorNN[iMom]);
+
+    util -> Delete();
   }
   
 
   // create TGraph to plot the pion efficiencies
-  TGraph *gEffisLQ=0x0, *gEffisNN=0x0;
-  TGraphErrors *gEffisErrLQ=0x0, *gEffisErrNN=0x0;
+  TGraphErrors *gEffisLQ=0x0, *gEffisNN=0x0;
+  gEffisLQ = (TGraphErrors*)fContainer->At(kGraphLQ);
+  gEffisNN = (TGraphErrors*)fContainer->At(kGraphNN);
 
-  gEffisLQ = (TGraph*)fContainer->At(kGraphLQ);
-  gEffisErrLQ = (TGraphErrors*)fContainer->At(kGraphLQerr);
-  gEffisNN = (TGraph*)fContainer->At(kGraphNN);
-  gEffisErrNN = (TGraphErrors*)fContainer->At(kGraphNNerr);
 
   for(Int_t iBin = 0; iBin < AliTRDCalPID::kNMom; iBin++){
+
     Float_t momentum = AliTRDCalPID::GetMomentum(iBin);
     gEffisLQ->SetPoint(iBin, momentum, PionEffiLQ[iBin]);
-    gEffisErrLQ->SetPoint(iBin, momentum, PionEffiLQ[iBin]);
-    gEffisErrLQ->SetPointError(iBin, 0., PionEffiErrorLQ[iBin]);
+    gEffisLQ->SetPointError(iBin, 0., PionEffiErrorLQ[iBin]);
 
     gEffisNN->SetPoint(iBin, momentum, PionEffiNN[iBin]);
-    gEffisErrNN->SetPoint(iBin, momentum, PionEffiNN[iBin]);
-    gEffisErrNN->SetPointError(iBin, 0., PionEffiErrorNN[iBin]);
+    gEffisNN->SetPointError(iBin, 0., PionEffiErrorNN[iBin]);
   }
 
-  gEffisLQ -> SetNameTitle("gEffisLQ", "Efficiencies of the 2-dim LQ method");
-  gEffisErrLQ -> SetNameTitle("gEffisLQErr", "Efficiencies and Errors of the 2-dim LQ method");
-  gEffisNN -> SetNameTitle("gEffisNN", "Efficiencies of the NN method");
-  gEffisErrNN -> SetNameTitle("gEffisNNErr", "Efficiencies and Errors of the NN method");
 
-  fNRefFigures = 2;
+  fNRefFigures = 1;
   return kTRUE; // testing protection
 }
 
@@ -464,138 +471,3 @@ void AliTRDpidChecker::Terminate(Option_t *)
 }
 
 
-//________________________________________________________________________
-Double_t  AliTRDpidChecker::GetPionEfficiency(Int_t Index1, Int_t Index2){
-
-  const Float_t Multi = 0.9;     // electron efficiency
-  Int_t abin, bbin;              
-  Double_t SumElecs, SumPions;   // integrated sum of elecs and pions in the histos
-  Double_t aBinSum, bBinSum;     // content of a single bin in the histos
-  
-  TH1F *Histo1 = (TH1F*)fContainer->At(Index1);  // electron histo
-  TH1F *Histo2 = (TH1F*)fContainer->At(Index2);  // pion histo
-
-
-  if(Index1 >= kNNlikelihood)              // print the correct momentum index for neural networks
-    Index1 = Index1 - kNNlikelihood;
-  SumElecs = 0.;
-  if(!(Histo1 -> GetEntries() > 0 && Histo2 -> GetEntries() > 0)){
-    if(fDebugLevel>=2) Printf("AliTRDpidChecker::GetPionEfficiency : Warning: Histo momentum intervall %d has no Entries!", Index1);
-    return -1.;
-  }
-  Histo1 -> Scale(1./Histo1->GetEntries());
-  Histo2 -> Scale(1./Histo2->GetEntries());
-
-
-  // calculate threshold for pion efficiency
-  for (abin = (Histo1 -> GetNbinsX()); abin >= 0; abin--){  
-    aBinSum = 0;
-    aBinSum = Histo1 -> GetBinContent(abin);
-    if(!(aBinSum == 0)){
-      SumElecs = SumElecs + aBinSum;
-    }
-    if (SumElecs >= Multi){
-      break;
-    }
-  }
-
-    
-  // calculate pion efficiency
-  SumPions = 0.;
-  for (bbin = (Histo2 -> GetNbinsX()); bbin >= abin; bbin--){  
-    bBinSum = 0;
-    bBinSum = Histo2 -> GetBinContent(bbin);
-    if(!(bBinSum == 0)){
-      SumPions = SumPions + bBinSum;
-    }
-  }
-  
-  
-  // print the electron efficiency and its cuts
-  if(fDebugLevel>=1) Printf("Cut for momentum intervall %d and electron efficiency of %f at: %5.3f", Index1, SumElecs, Histo1 -> GetBinCenter(abin));
-  if(fDebugLevel>=1) Printf("(%.0f electrons and %.0f pions)",Histo1 -> GetEntries(), Histo2 -> GetEntries());
-
-
-  // return the pion efficiency
-  return SumPions;
-
-} 
-  
-
-//________________________________________________________________________
-Double_t  AliTRDpidChecker::GetError(Int_t Index1, Int_t Index2){
-
-  
-//   const Int_t kBins = 12000;         // binning of the histos and eficiency calculation
-  const Float_t Multi = 0.9;                          // electron efficiency
-  Int_t abinE, bbinE, cbinE;                    
-  Double_t SumElecsE[kBins], SumPionsE[kBins];  // array of the integrated sum in each bin
-  Double_t aBinSumE, bBinSumE;                  // content of a single bin
-  Double_t EleEffi, PioEffi;                    // electron and pion efficiency
-  Bool_t bFirst = 1;                            // checks if threshold is crossed for the first time
-  Double_t fError = 0.;                         // error of the pion efficiency
-
-
-  TH1F *Histo1 = (TH1F*)fContainer->At(Index1);  // electron histo
-  TH1F *Histo2 = (TH1F*)fContainer->At(Index2);  // pion histo
-
-  if(Index1 > kNNlikelihood)              // print the correct momentum index for neural networks
-    Index1 = Index1 - kNNlikelihood;
-  if(!(Histo1 -> GetEntries() > 0 && Histo2 -> GetEntries() > 0)){
-    if(fDebugLevel>=2) Printf("Warning: Histo momentum intervall %d has no Entries!", Index1);
-    return -1.;
-  }
-
-  for(Int_t iBin = 0; iBin < kBins; iBin++){
-    SumElecsE[iBin] = 0.;
-    SumPionsE[iBin] = 0.;
-  }
-
-  EleEffi = 0.;
-  PioEffi = 0.;
-  cbinE = -1;
-
-
-  // calculate electron efficiency of each bin
-  for (abinE = (Histo1 -> GetNbinsX())-2; abinE >= 0; abinE--){  
-    aBinSumE = 0;
-    aBinSumE = Histo1 -> GetBinContent(abinE);
-    
-    SumElecsE[abinE] = SumElecsE[abinE+1] + aBinSumE;
-    if((SumElecsE[abinE] >= Multi) && (bFirst == 1)){
-      bFirst = 0;
-      cbinE = abinE;
-      EleEffi = (SumElecsE[cbinE]); 
-    }
-  }
-  
-
-  // calculate pion efficiency of each bin
-  for (bbinE = (Histo2 -> GetNbinsX())-2; bbinE >= abinE; bbinE--){    
-    bBinSumE = 0;
-    bBinSumE = Histo2 -> GetBinContent(bbinE);
-
-    SumPionsE[bbinE] = SumPionsE[bbinE+1] + bBinSumE;
-    if(bbinE == cbinE){
-      PioEffi = (SumPionsE[cbinE]);
-    }
-  }
-  
-
-  // pion efficiency vs electron efficiency
-  TGraph *gEffis = new TGraph(kBins, SumElecsE, SumPionsE);
-
-  // use fit function to get derivate of the TGraph for error calculation
-  TF1 *f1 = new TF1("f1","[0]*x*x+[1]*x+[2]", Multi-.05, Multi+.05);
-  gEffis -> Fit("f1","Q","",Multi-.05, Multi+.05);
-  if(fDebugLevel>=1) Printf("Derivative at %4.2f : %f", Multi, f1 -> Derivative(Multi));
-
-  // return the error of the pion efficiency
-  if(((1.-PioEffi) < 0) || ((1.-EleEffi) < 0)){
-    if(fDebugLevel>=2) Printf("AliTRDpidChecker::GetError : Warning: EleEffi or PioEffi > 1. Error can not be calculated. Please increase statistics or check your simulation!");
-    return -1;
-  }
-  fError = sqrt(((1/Histo2 -> GetEntries())*PioEffi*(1-PioEffi))+((f1 -> Derivative(Multi))*(f1 -> Derivative(Multi))*(1/Histo1 -> GetEntries())*EleEffi*(1-EleEffi)));
-  return fError;
-}
-
index d57290d..57826b0 100644 (file)
@@ -25,16 +25,14 @@ class AliTRDpidChecker : public AliTRDrecoTask
 {
 
   enum{
-    kLQlikelihood    = 0                                           // place for 2-dim LQ electron likelihood distributions
-    ,kNNlikelihood = 1 * AliTRDCalPID::kNMom * AliPID::kSPECIES  // place for NN electron likelihood distributions
-    ,kdEdx         = 2 * AliTRDCalPID::kNMom * AliPID::kSPECIES  // place for the dE/dx spectra
-    ,kPH           = 3 * AliTRDCalPID::kNMom * AliPID::kSPECIES  // place for pulse height spectra
-    ,kMomentum     = 4 * AliTRDCalPID::kNMom * AliPID::kSPECIES  // place for the momentum distribution
-    ,kMomentumBin  = kMomentum +1                                // place for the momentum distribution
-    ,kGraphLQ      = kMomentumBin +1                             // place for the 2-dim LQ pion efficiencies
-    ,kGraphLQerr   = kGraphLQ +1                                 // place for the 2-dim LQ pion efficiency errors
-    ,kGraphNN      = kGraphLQerr +1                              // place for the NN pion efficiencies
-    ,kGraphNNerr   = kGraphNN +1                                 // place for the NN pion efficiency errors
+    kLQlikelihood    = 0     // place for 2-dim LQ electron likelihood distributions
+    ,kNNlikelihood   = 1     // place for NN electron likelihood distributions
+    ,kdEdx           = 2     // place for the dE/dx spectra
+    ,kPH             = 3     // place for pulse height spectra
+    ,kMomentum       = 4     // place for the momentum distribution
+    ,kMomentumBin    = 5     // place for the momentum distribution
+    ,kGraphLQ        = 6     // place for the 2-dim LQ pion efficiencies
+    ,kGraphNN        = 7     // place for the NN pion efficiencies
   };
 
   enum{
@@ -56,16 +54,8 @@ private:
   AliTRDpidChecker(const AliTRDpidChecker&);               // not implemented
   AliTRDpidChecker& operator=(const AliTRDpidChecker&);    // not implemented
 
-  Double_t GetPionEfficiency(Int_t Index1, Int_t Index2);  // calculates the pion efficiency
-  Double_t GetError(Int_t Index1, Int_t Index2);           // calculates the error
-  
-
   AliTRDReconstructor *fReconstructor;     //! reconstructor needed for recalculation the PID
 
-  enum{
-    kBins = 12001                // binning of the likelihood histograms
-  };
-
   ClassDef(AliTRDpidChecker, 1); // TRD PID checker
 };