TRD/qaRec/AliTRDpidChecker.cxx

   1 #include "TPDGCode.h"
   2 #include "TF1.h"
   3 #include "TH1F.h"
   4 #include "TH1D.h"
   5 #include "TH2F.h"
   6 #include "TProfile.h"
   7 #include "TProfile2D.h"
   8 #include "TGraph.h"
   9 #include "TGraphErrors.h"
  10
  11 #include <TClonesArray.h>
  12 #include <TObjArray.h>
  13 #include <TList.h>
  14
  15 // #include "AliPID.h"
  16 #include "AliESDEvent.h"
  17 #include "AliESDInputHandler.h"
  18 #include "AliTrackReference.h"
  19
  20 #include "AliAnalysisTask.h"
  21 // #include "AliAnalysisManager.h"
  22
  23 #include "AliTRDtrackerV1.h"
  24 #include "AliTRDtrackV1.h"
  25 #include "AliTRDcluster.h"
  26 #include "AliTRDReconstructor.h"
  27 #include "AliCDBManager.h"
  28 // #include "../Cal/AliTRDCalPID.h"
  29 #include "AliTRDpidUtil.h"
  30
  31
  32 #include "AliTRDpidChecker.h"
  33 #include "AliTRDtrackInfo/AliTRDtrackInfo.h"
  34
  35 // calculate pion efficiency at 90% electron efficiency for 11 momentum bins
  36 // this task should be used with simulated data only
  37
  38 ClassImp(AliTRDpidChecker)
  39
  40 //________________________________________________________________________
  41 AliTRDpidChecker::AliTRDpidChecker()
  42   :AliTRDrecoTask("PID", "PID Checker")
  43   ,fReconstructor(0x0)
  44 {
  45   //
  46   // Default constructor
  47   //
  48
  49   fReconstructor = new AliTRDReconstructor();
  50   fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
  51 }
  52
  53
  54 //________________________________________________________________________
  55 AliTRDpidChecker::~AliTRDpidChecker()
  56 {
  57   if(fReconstructor) delete fReconstructor;
  58 }
  59
  60
  61 //________________________________________________________________________
  62 void AliTRDpidChecker::CreateOutputObjects()
  63 {
  64   // Create histograms
  65   // Called once
  66
  67   OpenFile(0, "RECREATE");
  68   Int_t xBins = AliPID::kSPECIES*AliTRDCalPID::kNMom;
  69   fContainer = new TObjArray();
  70   fContainer -> Expand(kGraphNN + 1);
  71
  72   const Float_t epsilon = 1./(2*(AliTRDpidUtil::kBins-1));     // get nice histos with bin center at 0 and 1
  73
  74   // histos of the electron probability of all 5 particle species and 11 momenta for the 2-dim LQ method
  75   fContainer->AddAt(
  76     new TH2F("PID_LQ", "",
  77       xBins, -0.5, xBins - 0.5,
  78       AliTRDpidUtil::kBins, 0.-epsilon, 1.+epsilon)
  79   ,kLQlikelihood);
  80
  81
  82   // histos of the electron probability of all 5 particle species and 11 momenta for the neural network method
  83   fContainer->AddAt(
  84     new TH2F("PID_NN", "",
  85       xBins, -0.5, xBins - 0.5,
  86       AliTRDpidUtil::kBins, 0.-epsilon, 1.+epsilon)
  87   ,kNNlikelihood);
  88
  89   // histos of the dE/dx distribution for all 5 particle species and 11 momenta
  90   fContainer->AddAt(
  91     new TH2F("dEdx", "",
  92       xBins, -0.5, xBins - 0.5,
  93       200, 0, 10000)
  94     ,kdEdx);
  95
  96   // histos of the pulse height distribution for all 5 particle species and 11 momenta
  97   fContainer->AddAt(
  98     new TProfile2D("PH", "",
  99       xBins, -0.5, xBins - 0.5,
 100       AliTRDtrackerV1::GetNTimeBins(), -0.5, AliTRDtrackerV1::GetNTimeBins() - 0.5)
 101     ,kPH);
 102
 103
 104   // momentum distributions - absolute and in momentum bins
 105   fContainer->AddAt(new TH1F("hMom", "momentum distribution", 100, 0., 12.),kMomentum);
 106   fContainer->AddAt(new TH1F("hMomBin", "momentum distribution in momentum bins", AliTRDCalPID::kNMom, 0.5, 11.5),kMomentumBin);
 107
 108
 109   // TGraph of the pion efficiencies
 110
 111   TGraphErrors *gEffisLQ = 0x0;
 112   TGraphErrors *gEffisNN = 0x0;
 113
 114   fContainer->AddAt(gEffisLQ = new TGraphErrors(), kGraphLQ);
 115   gEffisLQ->SetLineColor(kBlue);
 116   gEffisLQ->SetMarkerColor(kBlue);
 117   gEffisLQ->SetMarkerStyle(29);
 118
 119   fContainer -> AddAt(gEffisNN = new TGraphErrors(),kGraphNN);
 120   gEffisNN->SetLineColor(kRed);
 121   gEffisNN->SetMarkerColor(kRed);
 122   gEffisNN->SetMarkerStyle(29);
 123
 124   gEffisLQ -> SetNameTitle("gEffisLQErr", "Efficiencies and Errors of the 2-dim LQ method");
 125   gEffisNN -> SetNameTitle("gEffisNNErr", "Efficiencies and Errors of the NN method");
 126
 127 }
 128
 129 //________________________________________________________________________
 130 void AliTRDpidChecker::Exec(Option_t *)
 131 {
 132   // Main loop
 133   // Called for each event
 134
 135
 136 //   if(!AliTracker::GetFieldMap()){
 137 //     AliMagFMaps* field = new AliMagFMaps("Maps","Maps", 2, 1., 10., AliMagFMaps::k5kG);
 138 //     AliTracker::SetFieldMap(field, kTRUE);
 139 //   }
 140
 141   TH2F *hPIDLQ;
 142   TH2F *hPIDNN;
 143   TH2F *hdEdx;
 144   TProfile2D *hPH;
 145
 146   hPIDLQ = (TH2F*)fContainer->At(kLQlikelihood);
 147   hPIDNN = (TH2F*)fContainer->At(kNNlikelihood);
 148   hdEdx  = (TH2F*)fContainer->At(kdEdx);
 149   hPH    = (TProfile2D*)fContainer->At(kPH);
 150
 151   TH1F *hMom    = (TH1F*)fContainer->At(kMomentum);
 152   TH1F *hMomBin = (TH1F*)fContainer->At(kMomentumBin);
 153
 154   Int_t labelsacc[10000];
 155   memset(labelsacc, 0, sizeof(Int_t) * 10000);
 156
 157   Float_t mom;
 158   ULong_t status;
 159   Int_t nTRD = 0;
 160   Float_t *fdEdx;
 161
 162   AliTRDtrackInfo     *track = 0x0;
 163   AliTRDtrackV1    *TRDtrack = 0x0;
 164   AliTrackReference     *ref = 0x0;
 165   AliExternalTrackParam *esd = 0x0;
 166
 167   AliTRDseedV1 *TRDtracklet[AliTRDgeometry::kNlayer];
 168   for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++)
 169     TRDtracklet[iChamb] = 0x0;
 170
 171   AliTRDcluster *TRDcluster = 0x0;
 172
 173   AliTRDpidUtil *util = new AliTRDpidUtil();
 174   for(Int_t itrk=0; itrk<fTracks->GetEntriesFast(); itrk++){
 175     track = (AliTRDtrackInfo*)fTracks->UncheckedAt(itrk);
 176     if(!track->HasESDtrack()) continue;
 177     status = track->GetStatus();
 178     if(!(status&AliESDtrack::kTPCout)) continue;
 179
 180     if(!(TRDtrack = track->GetTRDtrack())) continue;
 181     //&&(track->GetNumberOfClustersRefit()
 182
 183     // use only tracks that hit 6 chambers
 184     if(!(TRDtrack->GetNumberOfTracklets() == AliTRDgeometry::kNlayer)) continue;
 185
 186     ref = track->GetTrackRef(0);
 187     esd = track->GetOuterParam();
 188     mom = ref ? ref->P(): esd->P();
 189
 190     labelsacc[nTRD] = track->GetLabel();
 191     nTRD++;
 192
 193
 194     // set the 11 momentum bins
 195     Int_t iMomBin = -1;
 196     iMomBin = util -> GetMomentumBin(mom);
 197     if(fDebugLevel>=4) Printf("MomBin[%d] MomTot[%f]", iMomBin, mom);
 198
 199
 200     // fill momentum histo to have the momentum distribution
 201     hMom -> Fill(mom);
 202     hMomBin -> Fill(iMomBin);
 203
 204
 205     // set the reconstructor
 206     TRDtrack -> SetReconstructor(fReconstructor);
 207
 208
 209     // if no monte carlo data available -> use TRDpid
 210     if(!HasMCdata()){
 211       fReconstructor -> SetOption("nn");
 212       TRDtrack -> CookPID();
 213       if(TRDtrack -> GetPID(0) > TRDtrack -> GetPID(1) + TRDtrack -> GetPID(2)  + TRDtrack -> GetPID(3) + TRDtrack -> GetPID(4)){
 214         track -> SetPDG(kElectron);
 215       }
 216       else if(TRDtrack -> GetPID(4) > TRDtrack -> GetPID(2)  && TRDtrack -> GetPID(4) > TRDtrack -> GetPID(3)  && TRDtrack -> GetPID(4) > TRDtrack -> GetPID(1)){
 217         track -> SetPDG(kProton);
 218       }
 219       else if(TRDtrack -> GetPID(3) > TRDtrack -> GetPID(1)  && TRDtrack -> GetPID(3) > TRDtrack -> GetPID(2)){
 220         track -> SetPDG(kKPlus);
 221       }
 222       else if(TRDtrack -> GetPID(1) > TRDtrack -> GetPID(2)){
 223         track -> SetPDG(kMuonPlus);
 224       }
 225       else{
 226         track -> SetPDG(kPiPlus);
 227       }
 228     }
 229
 230
 231     // calculate the probabilities for electron probability using 2-dim LQ, the deposited charge per chamber and the pulse height spectra and fill histograms
 232     fReconstructor -> SetOption("!nn");
 233     TRDtrack -> CookPID();
 234
 235     if(fDebugLevel>=4) Printf("PIDmethod[%d] Slices[%d] PDG[%d] LQLike[%f]", fReconstructor->GetPIDMethod(), fReconstructor->GetNdEdxSlices(), track->GetPDG(), TRDtrack -> GetPID(0));
 236
 237
 238     Float_t SumdEdx[AliTRDgeometry::kNlayer];
 239     for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
 240       TRDtracklet[iChamb] = TRDtrack -> GetTracklet(iChamb);
 241       SumdEdx[iChamb] = 0.;
 242       fdEdx = TRDtracklet[iChamb] -> GetdEdx();
 243       SumdEdx[iChamb] += fdEdx[0] + fdEdx[1] + fdEdx[2];
 244     }
 245
 246     switch(track->GetPDG()){
 247     case kElectron:
 248     case kPositron:
 249       hPIDLQ -> Fill(AliPID::kElectron * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
 250       for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
 251         hdEdx -> Fill(AliPID::kElectron * AliTRDCalPID::kNMom + iMomBin, SumdEdx[iChamb]);
 252         for(Int_t iClus = 0; iClus < AliTRDtrackerV1::GetNTimeBins(); iClus++){
 253           if(!(TRDcluster = (AliTRDcluster*)TRDtracklet[iChamb] -> GetClusters(iClus)))
 254             continue;
 255           hPH -> Fill(AliPID::kElectron * AliTRDCalPID::kNMom + iMomBin, TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus));
 256         }
 257       }
 258       break;
 259     case kMuonPlus:
 260     case kMuonMinus:
 261       hPIDLQ -> Fill(AliPID::kMuon * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
 262       for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
 263         hdEdx -> Fill(AliPID::kMuon * AliTRDCalPID::kNMom + iMomBin, SumdEdx[iChamb]);
 264         for(Int_t iClus = 0; iClus < AliTRDtrackerV1::GetNTimeBins(); iClus++){
 265           if(!(TRDcluster = (AliTRDcluster*)TRDtracklet[iChamb] -> GetClusters(iClus)))
 266             continue;
 267           hPH -> Fill(AliPID::kMuon * AliTRDCalPID::kNMom + iMomBin, TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus));
 268         }
 269       }
 270       break;
 271     case kPiPlus:
 272     case kPiMinus:
 273       hPIDLQ -> Fill(AliPID::kPion * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
 274       for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
 275         hdEdx -> Fill(AliPID::kPion * AliTRDCalPID::kNMom + iMomBin, SumdEdx[iChamb]);
 276         for(Int_t iClus = 0; iClus < AliTRDtrackerV1::GetNTimeBins(); iClus++){
 277           if(!(TRDcluster = (AliTRDcluster*)TRDtracklet[iChamb] -> GetClusters(iClus)))
 278             continue;
 279           hPH -> Fill(AliPID::kPion * AliTRDCalPID::kNMom + iMomBin, TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus));
 280         }
 281       }
 282       break;
 283     case kKPlus:
 284     case kKMinus:
 285       hPIDLQ -> Fill(AliPID::kKaon * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
 286       for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
 287         hdEdx -> Fill(AliPID::kKaon * AliTRDCalPID::kNMom + iMomBin, SumdEdx[iChamb]);
 288         for(Int_t iClus = 0; iClus < AliTRDtrackerV1::GetNTimeBins(); iClus++){
 289           if(!(TRDcluster = (AliTRDcluster*)TRDtracklet[iChamb] -> GetClusters(iClus)))
 290             continue;
 291           hPH -> Fill(AliPID::kKaon * AliTRDCalPID::kNMom + iMomBin, TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus));
 292         }
 293       }
 294       break;
 295     case kProton:
 296     case kProtonBar:
 297       hPIDLQ -> Fill(AliPID::kProton * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
 298       for(Int_t iChamb = 0; iChamb < AliTRDgeometry::kNlayer; iChamb++){
 299         hdEdx -> Fill(AliPID::kProton * AliTRDCalPID::kNMom + iMomBin, SumdEdx[iChamb]);
 300         for(Int_t iClus = 0; iClus < AliTRDtrackerV1::GetNTimeBins(); iClus++){
 301           if(!(TRDcluster = (AliTRDcluster*)TRDtracklet[iChamb] -> GetClusters(iClus)))
 302             continue;
 303           hPH -> Fill(AliPID::kProton * AliTRDCalPID::kNMom + iMomBin, TRDcluster -> GetLocalTimeBin(), TRDtracklet[iChamb] -> GetdQdl(iClus));
 304         }
 305       }
 306       break;
 307     }
 308
 309
 310     // calculate the probabilities and fill histograms for electrons using NN
 311     fReconstructor -> SetOption("nn");
 312     TRDtrack->CookPID();
 313
 314
 315     if(fDebugLevel>=4) Printf("PIDmethod[%d] Slices[%d] PDG[%d] NNLike[%f]", fReconstructor->GetPIDMethod(), fReconstructor->GetNdEdxSlices(), track->GetPDG(), TRDtrack -> GetPID(0));
 316
 317
 318     switch(track->GetPDG()){
 319     case kElectron:
 320     case kPositron:
 321       hPIDNN -> Fill(AliPID::kElectron * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
 322       break;
 323     case kMuonPlus:
 324     case kMuonMinus:
 325       hPIDNN -> Fill(AliPID::kMuon * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
 326       break;
 327     case kPiPlus:
 328     case kPiMinus:
 329       hPIDNN -> Fill(AliPID::kPion * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
 330       break;
 331     case kKPlus:
 332     case kKMinus:
 333       hPIDNN -> Fill(AliPID::kKaon * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
 334       break;
 335     case kProton:
 336     case kProtonBar:
 337       hPIDNN -> Fill(AliPID::kProton * AliTRDCalPID::kNMom + iMomBin, TRDtrack -> GetPID(0));
 338       break;
 339     }
 340
 341   }
 342
 343   util -> Delete();
 344   PostData(0, fContainer);
 345 }
 346
 347 //________________________________________________________
 348 void AliTRDpidChecker::GetRefFigure(Int_t ifig, Int_t &first, Int_t &last, Option_t *opt)
 349 {
 350   opt = "pl";
 351   switch(ifig){
 352   case 0:
 353     first = (Int_t)kGraphStart; last = first+2;
 354     break;
 355   default:
 356     first = (Int_t)kGraphStart; last = first+2;
 357     break;
 358   }
 359 }
 360
 361
 362 //________________________________________________________________________
 363 Bool_t AliTRDpidChecker::PostProcess()
 364 {
 365   // Draw result to the screen
 366   // Called once at the end of the query
 367
 368   if (!fContainer) {
 369     Printf("ERROR: list not available");
 370     return kFALSE;
 371   }
 372 //   return kTRUE; // testing protection
 373
 374
 375   // container for the pion efficiencies and the errors
 376   Double_t  PionEffiLQ[AliTRDCalPID::kNMom],
 377             PionEffiErrorLQ[AliTRDCalPID::kNMom],
 378             EleEffiLQ[AliTRDCalPID::kNMom],
 379             ThresholdLQ[AliTRDCalPID::kNMom];
 380
 381   Double_t  PionEffiNN[AliTRDCalPID::kNMom],
 382             PionEffiErrorNN[AliTRDCalPID::kNMom],
 383             EleEffiNN[AliTRDCalPID::kNMom],
 384             ThresholdNN[AliTRDCalPID::kNMom];
 385
 386   Float_t mom = 0.;
 387
 388   TH1D *Histo1=0x0, *Histo2=0x0;
 389
 390   TH2F *hPIDLQ=0x0, *hPIDNN=0x0;
 391   hPIDLQ = (TH2F*)fContainer->At(kLQlikelihood);
 392   hPIDNN = (TH2F*)fContainer->At(kNNlikelihood);
 393
 394   // calculate the pion efficiencies and the errors for 90% electron efficiency (2-dim LQ)
 395   for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
 396
 397     AliTRDpidUtil *util = new AliTRDpidUtil();
 398     mom = AliTRDCalPID::GetMomentum(iMom);
 399
 400     Histo1 = hPIDLQ -> ProjectionY("LQ_ele",AliTRDCalPID::kNMom*AliPID::kElectron+iMom+1,AliTRDCalPID::kNMom*AliPID::kElectron+iMom+1);
 401     Histo2 = hPIDLQ -> ProjectionY("LQ_pio",AliTRDCalPID::kNMom*AliPID::kPion+iMom+1,AliTRDCalPID::kNMom*AliPID::kPion+iMom+1);
 402
 403     util -> CalculatePionEffi(Histo1, Histo2);
 404
 405     PionEffiLQ[iMom] = util -> GetPionEfficiency();
 406     PionEffiErrorLQ[iMom] = util -> GetError();
 407     EleEffiLQ[iMom] = util -> GetCalcElectronEfficiency();
 408     ThresholdLQ[iMom] = util -> GetThreshold();
 409
 410     if(fDebugLevel>=1) Printf("Pion Efficiency for 2-dim LQ is : %f +/- %f\n\n", PionEffiLQ[iMom], PionEffiErrorLQ[iMom]);
 411     util -> Delete();
 412   }
 413
 414
 415
 416   // calculate the pion efficiencies and the errors for 90% electron efficiency (NN)
 417   for(Int_t iMom = 0; iMom < AliTRDCalPID::kNMom; iMom++){
 418
 419     AliTRDpidUtil *util = new AliTRDpidUtil();
 420     mom = AliTRDCalPID::GetMomentum(iMom);
 421
 422     Histo1 = hPIDNN -> ProjectionY("NN_ele",AliTRDCalPID::kNMom*AliPID::kElectron+iMom+1,AliTRDCalPID::kNMom*AliPID::kElectron+iMom+1);
 423     Histo2 = hPIDNN -> ProjectionY("NN_pio",AliTRDCalPID::kNMom*AliPID::kPion+iMom+1,AliTRDCalPID::kNMom*AliPID::kPion+iMom+1);
 424
 425     util -> CalculatePionEffi(Histo1, Histo2);
 426
 427     PionEffiNN[iMom] = util -> GetPionEfficiency();
 428     PionEffiErrorNN[iMom] = util -> GetError();
 429     EleEffiNN[iMom] = util -> GetCalcElectronEfficiency();
 430     ThresholdNN[iMom] = util -> GetThreshold();
 431
 432     if(fDebugLevel>=1) Printf("Pion Efficiency for NN is : %f +/- %f\n\n", PionEffiNN[iMom], PionEffiErrorNN[iMom]);
 433
 434     util -> Delete();
 435   }
 436
 437
 438   // create TGraph to plot the pion efficiencies
 439   TGraphErrors *gEffisLQ=0x0, *gEffisNN=0x0;
 440   gEffisLQ = (TGraphErrors*)fContainer->At(kGraphLQ);
 441   gEffisNN = (TGraphErrors*)fContainer->At(kGraphNN);
 442
 443
 444   for(Int_t iBin = 0; iBin < AliTRDCalPID::kNMom; iBin++){
 445
 446     Float_t momentum = AliTRDCalPID::GetMomentum(iBin);
 447     gEffisLQ->SetPoint(iBin, momentum, PionEffiLQ[iBin]);
 448     gEffisLQ->SetPointError(iBin, 0., PionEffiErrorLQ[iBin]);
 449
 450     gEffisNN->SetPoint(iBin, momentum, PionEffiNN[iBin]);
 451     gEffisNN->SetPointError(iBin, 0., PionEffiErrorNN[iBin]);
 452   }
 453
 454
 455   fNRefFigures = 1;
 456   return kTRUE; // testing protection
 457 }
 458
 459
 460 //________________________________________________________________________
 461 void AliTRDpidChecker::Terminate(Option_t *)
 462 {
 463   // Draw result to the screen
 464   // Called once at the end of the query
 465
 466   fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
 467   if (!fContainer) {
 468     Printf("ERROR: list not available");
 469     return;
 470   }
 471 }
 472
 473