TRD/qaRec/AliTRDpidChecker.cxx

   1 #include "TPDGCode.h"
   2 #include "TF1.h"
   3 #include "TH1F.h"
   4 #include "TH2F.h"
   5 #include "TGraph.h"
   6 #include "TGraphErrors.h"
   7
   8 #include <TClonesArray.h>
   9 #include <TObjArray.h>
  10 #include <TList.h>
  11
  12 #include "AliPID.h"
  13 #include "AliESDEvent.h"
  14 #include "AliESDInputHandler.h"
  15 #include "AliTrackReference.h"
  16
  17 #include "AliAnalysisTask.h"
  18 #include "AliAnalysisManager.h"
  19
  20 #include "AliTRDtrackV1.h"
  21 #include "AliTRDReconstructor.h"
  22 #include "AliCDBManager.h"
  23 #include "../Cal/AliTRDCalPID.h"
  24
  25
  26 #include "AliTRDpidChecker.h"
  27 #include "AliTRDtrackInfo/AliTRDtrackInfo.h"
  28
  29 // calculate pion efficiency at 90% electron efficiency for 11 momentum bins
  30 // this task should be used with simulated data only
  31
  32 ClassImp(AliTRDpidChecker)
  33
  34 //________________________________________________________________________
  35 AliTRDpidChecker::AliTRDpidChecker(const char *name)
  36   :AliAnalysisTask(name, "")
  37   ,fObjectContainer(0x0)
  38   ,fTracks(0x0)
  39   ,fReconstructor(0x0)
  40 //   ,fDebugLevel(1)
  41 //   ,fDebugStream(0x0)
  42 {
  43   //
  44   // Default constructor
  45   //
  46
  47   fReconstructor = new AliTRDReconstructor();
  48   fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
  49
  50   DefineInput(0, TObjArray::Class());
  51   DefineOutput(0, TObjArray::Class());
  52 }
  53
  54
  55 //________________________________________________________________________
  56 AliTRDpidChecker::~AliTRDpidChecker()
  57 {
  58   if(fObjectContainer){
  59     //fObjectContainer->Delete();
  60     //delete fObjectContainer;
  61   }
  62 }
  63
  64
  65
  66 //________________________________________________________________________
  67 void AliTRDpidChecker::ConnectInputData(Option_t *)
  68 {
  69   //
  70   // Connect input data
  71   //
  72
  73   fTracks = dynamic_cast<TObjArray*>(GetInputData(0));
  74 }
  75
  76
  77 //________________________________________________________________________
  78 void AliTRDpidChecker::CreateOutputObjects()
  79 {
  80   // Create histograms
  81   // Called once
  82
  83   OpenFile(0, "RECREATE");
  84   fObjectContainer = new TObjArray();
  85   fObjectContainer->Add(new TH1F("hMom", "momentum distribution", AliTRDCalPID::kNMom, 0.5, 11.5));
  86
  87
  88   // histos of the electron probability of all 5 particle species and 11 momenta for the 2-dim LQ method
  89   const Int_t kBins = 12000;         // binning of the histos and eficiency calculation
  90   for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
  91     for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
  92       fObjectContainer->Add(new TH1F(Form("PID%d_%d_LQ",iPart,iMom), Form("PID distribution for %d_%d", iPart, iMom), kBins, -0.1, 1.1));
  93     }
  94   }
  95
  96   // histos of the electron probability of all 5 particle species and 11 momenta for the neural network method
  97   for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
  98     for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
  99       fObjectContainer->Add(new TH1F(Form("PID%d_%d_NN",iPart,iMom), Form("PID distribution for %d_%d", iPart, iMom), kBins, -0.1, 1.1));
 100     }
 101   }
 102
 103   // frame and TGraph of the pion efficiencies
 104   fObjectContainer -> Add(new TH2F("hFrame", "", 10, 0.4, 12., 10, 0.0005, 0.1));
 105   fObjectContainer -> Add(new TGraph(AliTRDCalPID::kNMom));
 106   fObjectContainer -> Add(new TGraphErrors(AliTRDCalPID::kNMom));
 107   fObjectContainer -> Add(new TGraph(AliTRDCalPID::kNMom));
 108   fObjectContainer -> Add(new TGraphErrors(AliTRDCalPID::kNMom));
 109 }
 110
 111 //________________________________________________________________________
 112 void AliTRDpidChecker::Exec(Option_t *)
 113 {
 114   // Main loop
 115   // Called for each event
 116
 117
 118 //   if(!AliTracker::GetFieldMap()){
 119 //     AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 2, 1., 10., AliMagFMaps::k5kG);
 120 //     AliTracker::SetFieldMap(field, kTRUE);
 121 //   }
 122
 123   TH1F *hMom = (TH1F*)fObjectContainer->UncheckedAt(0);
 124   TH1F *hPIDLQ[AliPID::kSPECIES][AliTRDCalPID::kNMom];
 125   TH1F *hPIDNN[AliPID::kSPECIES][AliTRDCalPID::kNMom];
 126
 127   for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
 128     for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
 129       hPIDLQ[iPart][iMom] = (TH1F*)fObjectContainer->At(iPart*AliTRDCalPID::kNMom+iMom+1);
 130       hPIDNN[iPart][iMom] = (TH1F*)fObjectContainer->At(iPart*AliTRDCalPID::kNMom+iMom+1+AliPID::kSPECIES*AliTRDCalPID::kNMom);
 131     }
 132   }
 133
 134   Int_t labelsacc[10000];
 135   memset(labelsacc, 0, sizeof(Int_t) * 10000);
 136
 137   Float_t mom;
 138   ULong_t status;
 139   Int_t nTRD = 0;
 140
 141
 142   AliTRDtrackInfo     *track = 0x0;
 143   AliTRDtrackV1 *TRDtrack = 0x0;
 144   AliTrackReference     *ref = 0x0;
 145   AliExternalTrackParam *esd = 0x0;
 146   for(Int_t itrk=0; itrk<fTracks->GetEntriesFast(); itrk++){
 147     track = (AliTRDtrackInfo*)fTracks->UncheckedAt(itrk);
 148     if(!track->HasESDtrack()) continue;
 149     status = track->GetStatus();
 150     if(!(status&AliESDtrack::kTPCout)) continue;
 151
 152     if(!(TRDtrack = track->GetTRDtrack())) continue;
 153     //&&(track->GetNumberOfClustersRefit()
 154       // use only tracks taht hit 6 chambers
 155     if(!(TRDtrack->GetNumberOfTracklets() == AliTRDCalPID::kNPlane)) continue;
 156
 157     ref = track->GetTrackRef(0);
 158     esd = track->GetOuterParam();
 159     mom = ref ? ref->P(): esd->P();
 160
 161     labelsacc[nTRD] = track->GetLabel();
 162     nTRD++;
 163
 164     // fill momentum histo to have the momentum distribution
 165     hMom -> Fill(mom);
 166
 167
 168     // set the 11 momentum bins
 169     Int_t iMomBin = -1;
 170     if(mom < .7) iMomBin = 0;
 171     else if(mom < .9) iMomBin = 1;
 172     else if(mom < 1.25) iMomBin = 2;
 173     else if(mom < 1.75) iMomBin = 3;
 174     else if(mom < 2.5) iMomBin = 4;
 175     else if(mom < 3.5) iMomBin = 5;
 176     else if(mom < 4.5) iMomBin = 6;
 177     else if(mom < 5.5) iMomBin = 7;
 178     else if(mom < 7.0) iMomBin = 8;
 179     else if(mom < 9.0) iMomBin = 9;
 180     else  iMomBin = 10;
 181
 182     // set the reconstructor
 183     TRDtrack -> SetReconstructor(fReconstructor);
 184
 185
 186     // calculate the probabilities and fill histograms for electrons using 2-dim LQ
 187     fReconstructor -> SetOption("!nn");
 188     TRDtrack -> CookPID();
 189
 190     switch(track->GetPDG()){
 191     case kElectron:
 192     case kPositron:
 193       hPIDLQ[AliPID::kElectron][iMomBin] -> Fill(TRDtrack -> GetPID(0));
 194       break;
 195     case kMuonPlus:
 196     case kMuonMinus:
 197       hPIDLQ[AliPID::kMuon][iMomBin] -> Fill(TRDtrack -> GetPID(0));
 198       break;
 199     case kPiPlus:
 200     case kPiMinus:
 201       hPIDLQ[AliPID::kPion][iMomBin] -> Fill(TRDtrack -> GetPID(0));
 202       break;
 203     case kKPlus:
 204     case kKMinus:
 205       hPIDLQ[AliPID::kKaon][iMomBin] -> Fill(TRDtrack -> GetPID(0));
 206       break;
 207     case kProton:
 208     case kProtonBar:
 209       hPIDLQ[AliPID::kMuon][iMomBin] -> Fill(TRDtrack -> GetPID(0));
 210       break;
 211     }
 212
 213
 214     // calculate the probabilities and fill histograms for electrons using NN
 215     fReconstructor -> SetOption("nn");
 216     TRDtrack->CookPID();
 217     switch(track->GetPDG()){
 218     case kElectron:
 219     case kPositron:
 220       hPIDNN[AliPID::kElectron][iMomBin] -> Fill(TRDtrack -> GetPID(0));
 221       break;
 222     case kMuonPlus:
 223     case kMuonMinus:
 224       hPIDNN[AliPID::kMuon][iMomBin] -> Fill(TRDtrack -> GetPID(0));
 225       break;
 226     case kPiPlus:
 227     case kPiMinus:
 228       hPIDNN[AliPID::kPion][iMomBin] -> Fill(TRDtrack -> GetPID(0));
 229       break;
 230     case kKPlus:
 231     case kKMinus:
 232       hPIDNN[AliPID::kKaon][iMomBin] -> Fill(TRDtrack -> GetPID(0));
 233       break;
 234     case kProton:
 235     case kProtonBar:
 236       hPIDNN[AliPID::kMuon][iMomBin] -> Fill(TRDtrack -> GetPID(0));
 237       break;
 238     }
 239   }
 240
 241   PostData(0, fObjectContainer);
 242 }
 243
 244 //________________________________________________________________________
 245 void AliTRDpidChecker::Terminate(Option_t *)
 246 {
 247   // Draw result to the screen
 248   // Called once at the end of the query
 249
 250   fObjectContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
 251   if (!fObjectContainer) {
 252     Printf("ERROR: list not available");
 253     return;
 254   }
 255
 256
 257   // container for the pion efficiencies and the errors
 258   Double_t  PionEffiLQ[AliTRDCalPID::kNMom],
 259             PionEffiNN[AliTRDCalPID::kNMom],
 260             PionEffiErrorLQ[AliTRDCalPID::kNMom],
 261             PionEffiErrorNN[AliTRDCalPID::kNMom];
 262
 263
 264   // calculate the pion efficiencies and the errors for 90% electron efficiency (2-dim LQ)
 265   for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
 266     PionEffiLQ[iMom] = GetPionEfficiency(iMom+1,iMom+23);
 267     PionEffiErrorLQ[iMom] = GetError(iMom+1,iMom+23);
 268     Printf("Pion Efficiency for 2-dim LQ is : %f +/- %f\n\n", PionEffiLQ[iMom], PionEffiErrorLQ[iMom]);
 269   }
 270
 271   // calculate the pion efficiencies and the errors for 90% electron efficiency (NN)
 272   for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
 273     PionEffiNN[iMom] = GetPionEfficiency(iMom+56,iMom+78);
 274     PionEffiErrorNN[iMom] = GetError(iMom+56,iMom+78);
 275     Printf("Pion Efficiency for NN is : %f +/- %f\n\n", PionEffiNN[iMom], PionEffiErrorNN[iMom]);
 276   }
 277
 278
 279   // create TGraph to plot the pion efficiencies
 280   TGraph *gEffisLQ=0x0, *gEffisNN=0x0;
 281   TGraphErrors *gEffisErrLQ=0x0, *gEffisErrNN=0x0;
 282
 283   gEffisLQ = (TGraph*)fObjectContainer->At(112);
 284   gEffisErrLQ = (TGraphErrors*)fObjectContainer->At(113);
 285   gEffisNN = (TGraph*)fObjectContainer->At(114);
 286   gEffisErrNN = (TGraphErrors*)fObjectContainer->At(115);
 287
 288   for(Int_t iBin = 0; iBin < AliTRDCalPID::kNMom; iBin++){
 289     Float_t momentum = AliTRDCalPID::GetMomentum(iBin);
 290     gEffisLQ->SetPoint(iBin, momentum, PionEffiLQ[iBin]);
 291     gEffisErrLQ->SetPoint(iBin, momentum, PionEffiLQ[iBin]);
 292     gEffisErrLQ->SetPointError(iBin, 0., PionEffiErrorLQ[iBin]);
 293
 294     gEffisNN->SetPoint(iBin, momentum, PionEffiNN[iBin]);
 295     gEffisErrNN->SetPoint(iBin, momentum, PionEffiNN[iBin]);
 296     gEffisErrNN->SetPointError(iBin, 0., PionEffiErrorNN[iBin]);
 297   }
 298
 299   gEffisLQ -> SetNameTitle("gEffisLQ", "Efficiencies of the 2-dim LQ method");
 300   gEffisErrLQ -> SetNameTitle("gEffisLQErr", "Efficiencies and Errors of the 2-dim LQ method");
 301   gEffisNN -> SetNameTitle("gEffisNN", "Efficiencies of the NN method");
 302   gEffisErrNN -> SetNameTitle("gEffisNNErr", "Efficiencies and Errors of the NN method");
 303 }
 304
 305
 306 //________________________________________________________________________
 307 Double_t  AliTRDpidChecker::GetPionEfficiency(Int_t Index1, Int_t Index2){
 308
 309   Float_t Multi = 0.9;           // electron efficiency
 310   Int_t abin, bbin;
 311   Double_t SumElecs, SumPions;   // integrated sum of elecs and pions in the histos
 312   Double_t aBinSum, bBinSum;     // content of a single bin in the histos
 313
 314   TH1F *Histo1 = (TH1F*)fObjectContainer->At(Index1);  // electron histo
 315   TH1F *Histo2 = (TH1F*)fObjectContainer->At(Index2);  // pion histo
 316
 317
 318   SumElecs = 0.;
 319   if(!(Histo1 -> GetEntries() > 0 && Histo2 -> GetEntries() > 0)){
 320     Printf("Warning: Histo momentum intervall %d has no Entries!", Index1-1);
 321     return -1.;
 322   }
 323   Histo1 -> Scale(1./Histo1->GetEntries());
 324   Histo2 -> Scale(1./Histo2->GetEntries());
 325
 326
 327   // calculate threshold for pion efficiency
 328   for (abin = (Histo1 -> GetNbinsX()); abin >= 0; abin--){
 329     aBinSum = 0;
 330     aBinSum = Histo1 -> GetBinContent(abin);
 331     if(!(aBinSum == 0)){
 332       SumElecs = SumElecs + aBinSum;
 333     }
 334     if (SumElecs >= Multi){
 335       break;
 336     }
 337   }
 338
 339
 340   // calculate pion efficiency
 341   SumPions = 0.;
 342   for (bbin = (Histo2 -> GetNbinsX()); bbin >= abin; bbin--){
 343     bBinSum = 0;
 344     bBinSum = Histo2 -> GetBinContent(bbin);
 345     if(!(bBinSum == 0)){
 346       SumPions = SumPions + bBinSum;
 347     }
 348   }
 349
 350
 351   // print the electron efficiency and its cuts
 352   Printf("Cut for momentum intervall %d and electron efficiency of %f for: 0.%d", Index1-1, SumElecs, abin-1000);
 353   Printf("(%.0f electrons and %.0f pions)",Histo1 -> GetEntries(), Histo2 -> GetEntries());
 354
 355
 356   // return the pion efficiency
 357   return SumPions;
 358
 359 }
 360
 361
 362 //________________________________________________________________________
 363 Double_t  AliTRDpidChecker::GetError(Int_t Index1, Int_t Index2){
 364
 365
 366   const Int_t kBins = 12000;         // binning of the histos and eficiency calculation
 367   const Float_t Multi = 0.9;                          // electron efficiency
 368   Int_t abinE, bbinE, cbinE;
 369   Double_t SumElecsE[kBins], SumPionsE[kBins];  // array of the integrated sum in each bin
 370   Double_t aBinSumE, bBinSumE;                  // content of a single bin
 371   Double_t EleEffi, PioEffi;                    // electron and pion efficiency
 372   Bool_t bFirst = 1;                            // checks if threshold is crossed for the first time
 373   Double_t fError = 0.;                         // error of the pion efficiency
 374
 375
 376   TH1F *Histo1 = (TH1F*)fObjectContainer->At(Index1);  // electron histo
 377   TH1F *Histo2 = (TH1F*)fObjectContainer->At(Index2);  // pion histo
 378
 379   if(!(Histo1 -> GetEntries() > 0 && Histo2 -> GetEntries() > 0)){
 380     Printf("Warning: Histo momentum intervall %d has no Entries!", Index1-1);
 381     return -1.;
 382   }
 383
 384   for(Int_t iBin = 0; iBin < kBins; iBin++){
 385     SumElecsE[iBin] = 0.;
 386     SumPionsE[iBin] = 0.;
 387   }
 388
 389   EleEffi = 0.;
 390   PioEffi = 0.;
 391   cbinE = -1;
 392
 393
 394   // calculate electron efficiency of each bin
 395   for (abinE = (Histo1 -> GetNbinsX())-2; abinE >= 0; abinE--){
 396     aBinSumE = 0;
 397     aBinSumE = Histo1 -> GetBinContent(abinE);
 398
 399     SumElecsE[abinE] = SumElecsE[abinE+1] + aBinSumE;
 400     if((SumElecsE[abinE] >= Multi) && (bFirst == 1)){
 401       bFirst = 0;
 402       cbinE = abinE;
 403       EleEffi = (SumElecsE[cbinE]);
 404     }
 405   }
 406
 407
 408   // calculate pion efficiency of each bin
 409   for (bbinE = (Histo2 -> GetNbinsX())-2; bbinE >= abinE; bbinE--){
 410     bBinSumE = 0;
 411     bBinSumE = Histo2 -> GetBinContent(bbinE);
 412
 413     SumPionsE[bbinE] = SumPionsE[bbinE+1] + bBinSumE;
 414     if(bbinE == cbinE){
 415       PioEffi = (SumPionsE[cbinE]);
 416     }
 417   }
 418
 419
 420   // pion efficiency vs electron efficiency
 421   TGraph *gEffis = new TGraph(kBins, SumElecsE, SumPionsE);
 422
 423   // use fit function to get derivate of the TGraph for error calculation
 424   TF1 *f1 = new TF1("f1","[0]*x*x+[1]*x+[2]", Multi-.05, Multi+.05);
 425   gEffis -> Fit("f1","Q","",Multi-.05, Multi+.05);
 426   Printf("Derivative at %4.2f : %f", Multi, f1 -> Derivative(Multi));
 427
 428   // return the error of the pion efficiency
 429   fError = sqrt(((1/Histo2 -> GetEntries())*PioEffi*(1-PioEffi))+((f1 -> Derivative(Multi))*(f1 -> Derivative(Multi))*(1/Histo1 -> GetEntries())*EleEffi*(1-EleEffi)));
 430   return fError;
 431 }