TRD/qaRec/AliTRDpidChecker.cxx

   1 #include "TPDGCode.h"
   2 #include "TF1.h"
   3 #include "TH1F.h"
   4 #include "TH2F.h"
   5 #include "TProfile.h"
   6 #include "TGraph.h"
   7 #include "TGraphErrors.h"
   8
   9 #include <TClonesArray.h>
  10 #include <TObjArray.h>
  11 #include <TList.h>
  12
  13 // #include "AliPID.h"
  14 #include "AliESDEvent.h"
  15 #include "AliESDInputHandler.h"
  16 #include "AliTrackReference.h"
  17
  18 #include "AliAnalysisTask.h"
  19 // #include "AliAnalysisManager.h"
  20
  21 #include "AliTRDtrackerV1.h"
  22 #include "AliTRDtrackV1.h"
  23 #include "AliTRDcluster.h"
  24 #include "AliTRDReconstructor.h"
  25 #include "AliCDBManager.h"
  26 // #include "../Cal/AliTRDCalPID.h"
  27
  28
  29 #include "AliTRDpidChecker.h"
  30 #include "AliTRDtrackInfo/AliTRDtrackInfo.h"
  31
  32 // calculate pion efficiency at 90% electron efficiency for 11 momentum bins
  33 // this task should be used with simulated data only
  34
  35 ClassImp(AliTRDpidChecker)
  36
  37 //________________________________________________________________________
  38 AliTRDpidChecker::AliTRDpidChecker()
  39   :AliTRDrecoTask("PID", "PID Checker")
  40   ,fReconstructor(0x0)
  41 {
  42   //
  43   // Default constructor
  44   //
  45
  46   fReconstructor = new AliTRDReconstructor();
  47   fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
  48 }
  49
  50
  51 //________________________________________________________________________
  52 AliTRDpidChecker::~AliTRDpidChecker()
  53 {
  54   if(fReconstructor) delete fReconstructor;
  55 }
  56
  57
  58 //________________________________________________________________________
  59 void AliTRDpidChecker::CreateOutputObjects()
  60 {
  61   // Create histograms
  62   // Called once
  63
  64   OpenFile(0, "RECREATE");
  65   fContainer = new TObjArray();
  66   fContainer -> Expand(kGraphNNerr + 1);
  67
  68   const Float_t epsilon = 1/(2*(kBins-1));     // get nice histos with bin center at 0 and 1
  69
  70   // histos of the electron probability of all 5 particle species and 11 momenta for the 2-dim LQ method
  71   for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
  72     for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
  73       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);
  74       if(fDebugLevel>=4) Printf("LQ Particle[%d] Momentum[%d] : [%d]", iPart, iMom, kLQlikelihood + iPart*AliTRDCalPID::kNMom + iMom);
  75     }
  76   }
  77
  78   // histos of the electron probability of all 5 particle species and 11 momenta for the neural network method
  79   for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
  80     for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
  81       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);
  82       if(fDebugLevel>=4) Printf("NN Particle[%d] Momentum[%d] : [%d]", iPart, iMom, kNNlikelihood + iPart*AliTRDCalPID::kNMom + iMom);
  83     }
  84   }
  85
  86   // histos of the dE/dx distribution for all 5 particle species and 11 momenta
  87   for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
  88     for(Int_t iMom = 0;  iMom < AliTRDCalPID::kNMom; iMom++){
  89       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);
  90     }
  91   }
  92
  93   // histos of the pulse height distribution for all 5 particle species and 11 momenta
  94   for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
  95     for(Int_t iMom = 0;  iMom < AliTRDCalPID::kNMom; iMom++){
  96       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);
  97     }
  98   }
  99
 100   // momentum distributions - absolute and in momentum bins
 101   fContainer->AddAt(new TH1F("hMom", "momentum distribution", 100, 0., 12.),kMomentum);
 102   fContainer->AddAt(new TH1F("hMomBin", "momentum distribution in momentum bins", AliTRDCalPID::kNMom, 0.5, 11.5),kMomentumBin);
 103
 104
 105   // TGraph of the pion efficiencies
 106
 107   TGraph *gEffisLQ = new TGraph(AliTRDCalPID::kNMom);
 108   gEffisLQ -> SetLineColor(4);
 109   gEffisLQ -> SetMarkerColor(4);
 110   gEffisLQ -> SetMarkerStyle(29);
 111
 112   TGraphErrors *gEffisErrLQ = new TGraphErrors(AliTRDCalPID::kNMom);
 113   gEffisErrLQ -> SetLineColor(4);
 114   gEffisErrLQ -> SetMarkerColor(4);
 115   gEffisErrLQ -> SetMarkerStyle(29);
 116
 117   TGraph *gEffisNN = new TGraph(AliTRDCalPID::kNMom);
 118   gEffisNN -> SetLineColor(2);
 119   gEffisNN -> SetMarkerColor(2);
 120   gEffisNN -> SetMarkerStyle(29);
 121
 122   TGraphErrors *gEffisErrNN = new TGraphErrors(AliTRDCalPID::kNMom);
 123   gEffisErrNN -> SetLineColor(2);
 124   gEffisErrNN -> SetMarkerColor(2);
 125   gEffisErrNN -> SetMarkerStyle(29);
 126
 127
 128   fContainer -> AddAt(gEffisLQ,kGraphLQ);
 129   fContainer -> AddAt(gEffisErrLQ,kGraphLQerr);
 130   fContainer -> AddAt(gEffisNN,kGraphNN);
 131   fContainer -> AddAt(gEffisErrNN,kGraphNNerr);
 132 }
 133
 134 //________________________________________________________________________
 135 void AliTRDpidChecker::Exec(Option_t *)
 136 {
 137   // Main loop
 138   // Called for each event
 139
 140
 141 //   if(!AliTracker::GetFieldMap()){
 142 //     AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 2, 1., 10., AliMagFMaps::k5kG);
 143 //     AliTracker::SetFieldMap(field, kTRUE);
 144 //   }
 145
 146   TH1F *hPIDLQ[AliPID::kSPECIES][AliTRDCalPID::kNMom];
 147   TH1F *hPIDNN[AliPID::kSPECIES][AliTRDCalPID::kNMom];
 148   TH1F *hdEdx[AliPID::kSPECIES][AliTRDCalPID::kNMom];
 149   TProfile *hPH[AliPID::kSPECIES][AliTRDCalPID::kNMom];
 150
 151   for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
 152     for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
 153       hPIDLQ[iPart][iMom] = (TH1F*)fContainer->At(kLQlikelihood + iPart*AliTRDCalPID::kNMom+iMom);
 154       hPIDNN[iPart][iMom] = (TH1F*)fContainer->At(kNNlikelihood + iPart*AliTRDCalPID::kNMom+iMom);
 155       hdEdx[iPart][iMom]  = (TH1F*)fContainer->At(kdEdx + iPart*AliTRDCalPID::kNMom+iMom);
 156       hPH[iPart][iMom]    = (TProfile*)fContainer->At(kPH + iPart*AliTRDCalPID::kNMom+iMom);
 157     }
 158   }
 159
 160   TH1F *hMom    = (TH1F*)fContainer->At(kMomentum);
 161   TH1F *hMomBin = (TH1F*)fContainer->At(kMomentumBin);
 162
 163   Int_t labelsacc[10000];
 164   memset(labelsacc, 0, sizeof(Int_t) * 10000);
 165
 166   Float_t mom;
 167   ULong_t status;
 168   Int_t nTRD = 0;
 169   Float_t *fdEdx;
 170
 171   AliTRDtrackInfo     *track = 0x0;
 172   AliTRDtrackV1    *TRDtrack = 0x0;
 173   AliTrackReference     *ref = 0x0;
 174   AliExternalTrackParam *esd = 0x0;
 175
 176   AliTRDseedV1 *TRDtracklet[AliTRDCalPID::kNPlane];
 177   for(Int_t iChamb = 0; iChamb < AliTRDCalPID::kNPlane; iChamb++)
 178     TRDtracklet[iChamb] = 0x0;
 179
 180   AliTRDcluster *TRDcluster = 0x0;
 181
 182   for(Int_t itrk=0; itrk<fTracks->GetEntriesFast(); itrk++){
 183     track = (AliTRDtrackInfo*)fTracks->UncheckedAt(itrk);
 184     if(!track->HasESDtrack()) continue;
 185     status = track->GetStatus();
 186     if(!(status&AliESDtrack::kTPCout)) continue;
 187
 188     if(!(TRDtrack = track->GetTRDtrack())) continue;
 189     //&&(track->GetNumberOfClustersRefit()
 190
 191     // use only tracks that hit 6 chambers
 192     if(!(TRDtrack->GetNumberOfTracklets() == AliTRDCalPID::kNPlane)) continue;
 193
 194     ref = track->GetTrackRef(0);
 195     esd = track->GetOuterParam();
 196     mom = ref ? ref->P(): esd->P();
 197
 198     labelsacc[nTRD] = track->GetLabel();
 199     nTRD++;
 200
 201
 202     // set the 11 momentum bins
 203     Int_t iMomBin = -1;
 204     if(mom < .7) iMomBin = 0;
 205     else if(mom < .9) iMomBin = 1;
 206     else if(mom < 1.25) iMomBin = 2;
 207     else if(mom < 1.75) iMomBin = 3;
 208     else if(mom < 2.5) iMomBin = 4;
 209     else if(mom < 3.5) iMomBin = 5;
 210     else if(mom < 4.5) iMomBin = 6;
 211     else if(mom < 5.5) iMomBin = 7;
 212     else if(mom < 7.0) iMomBin = 8;
 213     else if(mom < 9.0) iMomBin = 9;
 214     else  iMomBin = 10;
 215
 216     // fill momentum histo to have the momentum distribution
 217     hMom -> Fill(mom);
 218     hMomBin -> Fill(iMomBin);
 219
 220
 221     // set the reconstructor
 222     TRDtrack -> SetReconstructor(fReconstructor);
 223
 224
 225     // if no monte carlo data available -> use TRDpid
 226     if(!HasMCdata()){
 227       fReconstructor -> SetOption("nn");
 228       TRDtrack -> CookPID();
 229       if(TRDtrack -> GetPID(0) > TRDtrack -> GetPID(1) + TRDtrack -> GetPID(2)  + TRDtrack -> GetPID(3) + TRDtrack -> GetPID(4)){
 230         track -> SetPDG(kElectron);
 231       }
 232       else if(TRDtrack -> GetPID(4) > TRDtrack -> GetPID(2)  && TRDtrack -> GetPID(4) > TRDtrack -> GetPID(3)  && TRDtrack -> GetPID(4) > TRDtrack -> GetPID(1)){
 233         track -> SetPDG(kProton);
 234       }
 235       else if(TRDtrack -> GetPID(3) > TRDtrack -> GetPID(1)  && TRDtrack -> GetPID(3) > TRDtrack -> GetPID(2)){
 236         track -> SetPDG(kKPlus);
 237       }
 238       else if(TRDtrack -> GetPID(1) > TRDtrack -> GetPID(2)){
 239         track -> SetPDG(kMuonPlus);
 240       }
 241       else{
 242         track -> SetPDG(kPiPlus);
 243       }
 244     }
 245
 246
 247     // calculate the probabilities for electron probability using 2-dim LQ, the deposited charge per chamber and the pulse height spectra and fill histograms
 248     fReconstructor -> SetOption("!nn");
 249     TRDtrack -> CookPID();
 250
 251     if(fDebugLevel>=4) Printf("PIDmethod[%d] Slices[%d] PDG[%d] LQLike[%f]", fReconstructor->GetPIDMethod(), fReconstructor->GetNdEdxSlices(), track->GetPDG(), TRDtrack -> GetPID(0));
 252
 253
 254     Float_t SumdEdx[AliTRDCalPID::kNPlane];
 255     for(Int_t iChamb = 0; iChamb < AliTRDCalPID::kNPlane; iChamb++){
 256       TRDtracklet[iChamb] = TRDtrack -> GetTracklet(iChamb);
 257       SumdEdx[iChamb] = 0.;
 258       fdEdx = TRDtracklet[iChamb] -> GetdEdx();
 259       SumdEdx[iChamb] += fdEdx[0] + fdEdx[1] + fdEdx[2];
 260     }
 261
 262     switch(track->GetPDG()){
 263     case kElectron:
 264     case kPositron:
 265       hPIDLQ[AliPID::kElectron][iMomBin] -> Fill(TRDtrack -> GetPID(0));
 266       for(Int_t iChamb = 0; iChamb < AliTRDCalPID::kNPlane; iChamb++){
 267         hdEdx[AliPID::kElectron][iMomBin] -> Fill(SumdEdx[iChamb]);
 268         for(Int_t iClus = 0; iClus < AliTRDtrackerV1::GetNTimeBins(); iClus++){
 269           if(!(TRDcluster = (AliTRDcluster*)TRDtracklet[iChamb] -> GetClusters(iClus)))
 270             continue;
 271           hPH[AliPID::kElectron][iMomBin] -> Fill(TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus));
 272         }
 273       }
 274       break;
 275     case kMuonPlus:
 276     case kMuonMinus:
 277       hPIDLQ[AliPID::kMuon][iMomBin] -> Fill(TRDtrack -> GetPID(0));
 278       for(Int_t iChamb = 0; iChamb < AliTRDCalPID::kNPlane; iChamb++){
 279         hdEdx[AliPID::kMuon][iMomBin] -> Fill(SumdEdx[iChamb]);
 280         for(Int_t iClus = 0; iClus < AliTRDtrackerV1::GetNTimeBins(); iClus++){
 281           if(!(TRDcluster = (AliTRDcluster*)TRDtracklet[iChamb] -> GetClusters(iClus)))
 282             continue;
 283           hPH[AliPID::kMuon][iMomBin] -> Fill(TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus));
 284         }
 285       }
 286       break;
 287     case kPiPlus:
 288     case kPiMinus:
 289       hPIDLQ[AliPID::kPion][iMomBin] -> Fill(TRDtrack -> GetPID(0));
 290       for(Int_t iChamb = 0; iChamb < AliTRDCalPID::kNPlane; iChamb++){
 291         hdEdx[AliPID::kPion][iMomBin] -> Fill(SumdEdx[iChamb]);
 292         for(Int_t iClus = 0; iClus < AliTRDtrackerV1::GetNTimeBins(); iClus++){
 293           if(!(TRDcluster = (AliTRDcluster*)TRDtracklet[iChamb] -> GetClusters(iClus)))
 294             continue;
 295           hPH[AliPID::kPion][iMomBin] -> Fill(TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus));
 296         }
 297       }
 298       break;
 299     case kKPlus:
 300     case kKMinus:
 301       hPIDLQ[AliPID::kKaon][iMomBin] -> Fill(TRDtrack -> GetPID(0));
 302       for(Int_t iChamb = 0; iChamb < AliTRDCalPID::kNPlane; iChamb++){
 303         hdEdx[AliPID::kKaon][iMomBin] -> Fill(SumdEdx[iChamb]);
 304         for(Int_t iClus = 0; iClus < AliTRDtrackerV1::GetNTimeBins(); iClus++){
 305           if(!(TRDcluster = (AliTRDcluster*)TRDtracklet[iChamb] -> GetClusters(iClus)))
 306             continue;
 307           hPH[AliPID::kKaon][iMomBin] -> Fill(TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus));
 308         }
 309       }
 310       break;
 311     case kProton:
 312     case kProtonBar:
 313       hPIDLQ[AliPID::kProton][iMomBin] -> Fill(TRDtrack -> GetPID(0));
 314       for(Int_t iChamb = 0; iChamb < AliTRDCalPID::kNPlane; iChamb++){
 315         hdEdx[AliPID::kProton][iMomBin] -> Fill(SumdEdx[iChamb]);
 316         for(Int_t iClus = 0; iClus < AliTRDtrackerV1::GetNTimeBins(); iClus++){
 317           if(!(TRDcluster = (AliTRDcluster*)TRDtracklet[iChamb] -> GetClusters(iClus)))
 318             continue;
 319           hPH[AliPID::kProton][iMomBin] -> Fill(TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus));
 320         }
 321       }
 322       break;
 323     }
 324
 325
 326     // calculate the probabilities and fill histograms for electrons using NN
 327     fReconstructor -> SetOption("nn");
 328     TRDtrack->CookPID();
 329
 330
 331     if(fDebugLevel>=4) Printf("PIDmethod[%d] Slices[%d] PDG[%d] NNLike[%f]", fReconstructor->GetPIDMethod(), fReconstructor->GetNdEdxSlices(), track->GetPDG(), TRDtrack -> GetPID(0));
 332
 333
 334     switch(track->GetPDG()){
 335     case kElectron:
 336     case kPositron:
 337       hPIDNN[AliPID::kElectron][iMomBin] -> Fill(TRDtrack -> GetPID(0));
 338       break;
 339     case kMuonPlus:
 340     case kMuonMinus:
 341       hPIDNN[AliPID::kMuon][iMomBin] -> Fill(TRDtrack -> GetPID(0));
 342       break;
 343     case kPiPlus:
 344     case kPiMinus:
 345       hPIDNN[AliPID::kPion][iMomBin] -> Fill(TRDtrack -> GetPID(0));
 346       break;
 347     case kKPlus:
 348     case kKMinus:
 349       hPIDNN[AliPID::kKaon][iMomBin] -> Fill(TRDtrack -> GetPID(0));
 350       break;
 351     case kProton:
 352     case kProtonBar:
 353       hPIDNN[AliPID::kProton][iMomBin] -> Fill(TRDtrack -> GetPID(0));
 354       break;
 355     }
 356   }
 357
 358   PostData(0, fContainer);
 359 }
 360
 361 //________________________________________________________
 362 void AliTRDpidChecker::GetRefFigure(Int_t ifig, Int_t &first, Int_t &last, Option_t */*opt*/)
 363 {
 364   switch(ifig){
 365   case 0:
 366     first = (Int_t)kGraphStart+1; last = first+1;
 367     break;
 368   case 1:
 369     first = (Int_t)kGraphStart+3; last = first+1;
 370     break;
 371   default:
 372     first = (Int_t)kGraphStart; last = first;
 373     break;
 374   }
 375 }
 376
 377
 378 //________________________________________________________________________
 379 Bool_t AliTRDpidChecker::PostProcess()
 380 {
 381   // Draw result to the screen
 382   // Called once at the end of the query
 383
 384   if (!fContainer) {
 385     Printf("ERROR: list not available");
 386     return kFALSE;
 387   }
 388 //   return kTRUE; // testing protection
 389
 390   // normalize the dE/dx histos
 391   const Int_t kNSpectra = AliPID::kSPECIES*AliTRDCalPID::kNMom;
 392   TH1F *hdEdx[kNSpectra];
 393   for(Int_t iHist = 0; iHist < kNSpectra; iHist++){
 394     hdEdx[iHist] = (TH1F*)fContainer->At(kdEdx + iHist);
 395     if(hdEdx[iHist] -> GetEntries() > 0)
 396       hdEdx[iHist] -> Scale(1./hdEdx[iHist] -> GetEntries());
 397     else continue;
 398   }
 399
 400
 401   // container for the pion efficiencies and the errors
 402   Double_t  PionEffiLQ[AliTRDCalPID::kNMom],
 403             PionEffiNN[AliTRDCalPID::kNMom],
 404             PionEffiErrorLQ[AliTRDCalPID::kNMom],
 405             PionEffiErrorNN[AliTRDCalPID::kNMom];
 406
 407
 408   // calculate the pion efficiencies and the errors for 90% electron efficiency (2-dim LQ)
 409   for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
 410     PionEffiLQ[iMom] = GetPionEfficiency((kLQlikelihood + iMom) + (AliPID::kElectron * AliTRDCalPID::kNMom), (kLQlikelihood + iMom) + (AliPID::kPion * AliTRDCalPID::kNMom));
 411     PionEffiErrorLQ[iMom] = GetError((kLQlikelihood + iMom) + (AliPID::kElectron * AliTRDCalPID::kNMom), (kLQlikelihood + iMom) + (AliPID::kPion * AliTRDCalPID::kNMom));
 412     if(fDebugLevel>=1) Printf("Pion Efficiency for 2-dim LQ is : %f +/- %f\n\n", PionEffiLQ[iMom], PionEffiErrorLQ[iMom]);
 413   }
 414
 415   // calculate the pion efficiencies and the errors for 90% electron efficiency (NN)
 416   for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
 417     PionEffiNN[iMom] = GetPionEfficiency((kNNlikelihood + iMom) + (AliPID::kElectron * AliTRDCalPID::kNMom), (kNNlikelihood + iMom) + (AliPID::kPion * AliTRDCalPID::kNMom));
 418     PionEffiErrorNN[iMom] = GetError((kNNlikelihood + iMom) + (AliPID::kElectron * AliTRDCalPID::kNMom), (kNNlikelihood + iMom) + (AliPID::kPion * AliTRDCalPID::kNMom));
 419     if(fDebugLevel>=1) Printf("Pion Efficiency for NN is : %f +/- %f\n\n", PionEffiNN[iMom], PionEffiErrorNN[iMom]);
 420   }
 421
 422
 423   // create TGraph to plot the pion efficiencies
 424   TGraph *gEffisLQ=0x0, *gEffisNN=0x0;
 425   TGraphErrors *gEffisErrLQ=0x0, *gEffisErrNN=0x0;
 426
 427   gEffisLQ = (TGraph*)fContainer->At(kGraphLQ);
 428   gEffisErrLQ = (TGraphErrors*)fContainer->At(kGraphLQerr);
 429   gEffisNN = (TGraph*)fContainer->At(kGraphNN);
 430   gEffisErrNN = (TGraphErrors*)fContainer->At(kGraphNNerr);
 431
 432   for(Int_t iBin = 0; iBin < AliTRDCalPID::kNMom; iBin++){
 433     Float_t momentum = AliTRDCalPID::GetMomentum(iBin);
 434     gEffisLQ->SetPoint(iBin, momentum, PionEffiLQ[iBin]);
 435     gEffisErrLQ->SetPoint(iBin, momentum, PionEffiLQ[iBin]);
 436     gEffisErrLQ->SetPointError(iBin, 0., PionEffiErrorLQ[iBin]);
 437
 438     gEffisNN->SetPoint(iBin, momentum, PionEffiNN[iBin]);
 439     gEffisErrNN->SetPoint(iBin, momentum, PionEffiNN[iBin]);
 440     gEffisErrNN->SetPointError(iBin, 0., PionEffiErrorNN[iBin]);
 441   }
 442
 443   gEffisLQ -> SetNameTitle("gEffisLQ", "Efficiencies of the 2-dim LQ method");
 444   gEffisErrLQ -> SetNameTitle("gEffisLQErr", "Efficiencies and Errors of the 2-dim LQ method");
 445   gEffisNN -> SetNameTitle("gEffisNN", "Efficiencies of the NN method");
 446   gEffisErrNN -> SetNameTitle("gEffisNNErr", "Efficiencies and Errors of the NN method");
 447
 448   fNRefFigures = 2;
 449   return kTRUE; // testing protection
 450 }
 451
 452
 453 //________________________________________________________________________
 454 void AliTRDpidChecker::Terminate(Option_t *)
 455 {
 456   // Draw result to the screen
 457   // Called once at the end of the query
 458
 459   fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
 460   if (!fContainer) {
 461     Printf("ERROR: list not available");
 462     return;
 463   }
 464 }
 465
 466
 467 //________________________________________________________________________
 468 Double_t  AliTRDpidChecker::GetPionEfficiency(Int_t Index1, Int_t Index2){
 469
 470   const Float_t Multi = 0.9;     // electron efficiency
 471   Int_t abin, bbin;
 472   Double_t SumElecs, SumPions;   // integrated sum of elecs and pions in the histos
 473   Double_t aBinSum, bBinSum;     // content of a single bin in the histos
 474
 475   TH1F *Histo1 = (TH1F*)fContainer->At(Index1);  // electron histo
 476   TH1F *Histo2 = (TH1F*)fContainer->At(Index2);  // pion histo
 477
 478
 479   if(Index1 >= kNNlikelihood)              // print the correct momentum index for neural networks
 480     Index1 = Index1 - kNNlikelihood;
 481   SumElecs = 0.;
 482   if(!(Histo1 -> GetEntries() > 0 && Histo2 -> GetEntries() > 0)){
 483     if(fDebugLevel>=2) Printf("AliTRDpidChecker::GetPionEfficiency : Warning: Histo momentum intervall %d has no Entries!", Index1);
 484     return -1.;
 485   }
 486   Histo1 -> Scale(1./Histo1->GetEntries());
 487   Histo2 -> Scale(1./Histo2->GetEntries());
 488
 489
 490   // calculate threshold for pion efficiency
 491   for (abin = (Histo1 -> GetNbinsX()); abin >= 0; abin--){
 492     aBinSum = 0;
 493     aBinSum = Histo1 -> GetBinContent(abin);
 494     if(!(aBinSum == 0)){
 495       SumElecs = SumElecs + aBinSum;
 496     }
 497     if (SumElecs >= Multi){
 498       break;
 499     }
 500   }
 501
 502
 503   // calculate pion efficiency
 504   SumPions = 0.;
 505   for (bbin = (Histo2 -> GetNbinsX()); bbin >= abin; bbin--){
 506     bBinSum = 0;
 507     bBinSum = Histo2 -> GetBinContent(bbin);
 508     if(!(bBinSum == 0)){
 509       SumPions = SumPions + bBinSum;
 510     }
 511   }
 512
 513
 514   // print the electron efficiency and its cuts
 515   if(fDebugLevel>=1) Printf("Cut for momentum intervall %d and electron efficiency of %f at: %5.3f", Index1, SumElecs, Histo1 -> GetBinCenter(abin));
 516   if(fDebugLevel>=1) Printf("(%.0f electrons and %.0f pions)",Histo1 -> GetEntries(), Histo2 -> GetEntries());
 517
 518
 519   // return the pion efficiency
 520   return SumPions;
 521
 522 }
 523
 524
 525 //________________________________________________________________________
 526 Double_t  AliTRDpidChecker::GetError(Int_t Index1, Int_t Index2){
 527
 528
 529 //   const Int_t kBins = 12000;         // binning of the histos and eficiency calculation
 530   const Float_t Multi = 0.9;                          // electron efficiency
 531   Int_t abinE, bbinE, cbinE;
 532   Double_t SumElecsE[kBins], SumPionsE[kBins];  // array of the integrated sum in each bin
 533   Double_t aBinSumE, bBinSumE;                  // content of a single bin
 534   Double_t EleEffi, PioEffi;                    // electron and pion efficiency
 535   Bool_t bFirst = 1;                            // checks if threshold is crossed for the first time
 536   Double_t fError = 0.;                         // error of the pion efficiency
 537
 538
 539   TH1F *Histo1 = (TH1F*)fContainer->At(Index1);  // electron histo
 540   TH1F *Histo2 = (TH1F*)fContainer->At(Index2);  // pion histo
 541
 542   if(Index1 > kNNlikelihood)              // print the correct momentum index for neural networks
 543     Index1 = Index1 - kNNlikelihood;
 544   if(!(Histo1 -> GetEntries() > 0 && Histo2 -> GetEntries() > 0)){
 545     if(fDebugLevel>=2) Printf("Warning: Histo momentum intervall %d has no Entries!", Index1);
 546     return -1.;
 547   }
 548
 549   for(Int_t iBin = 0; iBin < kBins; iBin++){
 550     SumElecsE[iBin] = 0.;
 551     SumPionsE[iBin] = 0.;
 552   }
 553
 554   EleEffi = 0.;
 555   PioEffi = 0.;
 556   cbinE = -1;
 557
 558
 559   // calculate electron efficiency of each bin
 560   for (abinE = (Histo1 -> GetNbinsX())-2; abinE >= 0; abinE--){
 561     aBinSumE = 0;
 562     aBinSumE = Histo1 -> GetBinContent(abinE);
 563
 564     SumElecsE[abinE] = SumElecsE[abinE+1] + aBinSumE;
 565     if((SumElecsE[abinE] >= Multi) && (bFirst == 1)){
 566       bFirst = 0;
 567       cbinE = abinE;
 568       EleEffi = (SumElecsE[cbinE]);
 569     }
 570   }
 571
 572
 573   // calculate pion efficiency of each bin
 574   for (bbinE = (Histo2 -> GetNbinsX())-2; bbinE >= abinE; bbinE--){
 575     bBinSumE = 0;
 576     bBinSumE = Histo2 -> GetBinContent(bbinE);
 577
 578     SumPionsE[bbinE] = SumPionsE[bbinE+1] + bBinSumE;
 579     if(bbinE == cbinE){
 580       PioEffi = (SumPionsE[cbinE]);
 581     }
 582   }
 583
 584
 585   // pion efficiency vs electron efficiency
 586   TGraph *gEffis = new TGraph(kBins, SumElecsE, SumPionsE);
 587
 588   // use fit function to get derivate of the TGraph for error calculation
 589   TF1 *f1 = new TF1("f1","[0]*x*x+[1]*x+[2]", Multi-.05, Multi+.05);
 590   gEffis -> Fit("f1","Q","",Multi-.05, Multi+.05);
 591   if(fDebugLevel>=1) Printf("Derivative at %4.2f : %f", Multi, f1 -> Derivative(Multi));
 592
 593   // return the error of the pion efficiency
 594   if(((1.-PioEffi) < 0) || ((1.-EleEffi) < 0)){
 595     if(fDebugLevel>=2) Printf("AliTRDpidChecker::GetError : Warning: EleEffi or PioEffi > 1. Error can not be calculated. Please increase statistics or check your simulation!");
 596     return -1;
 597   }
 598   fError = sqrt(((1/Histo2 -> GetEntries())*PioEffi*(1-PioEffi))+((f1 -> Derivative(Multi))*(f1 -> Derivative(Multi))*(1/Histo1 -> GetEntries())*EleEffi*(1-EleEffi)));
 599   return fError;
 600 }
 601