PWG1/TRD/AliTRDefficiencyMC.cxx

   1 /**************************************************************************
   2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 /* $Id: AliTRDefficiencyMC.cxx 27496 2008-07-22 08:35:45Z cblume $ */
  17
  18 ////////////////////////////////////////////////////////////////////////////
  19 //                                                                        //
  20 //  Reconstruction QA                                                     //
  21 //                                                                        //
  22 //  Authors:                                                              //
  23 //    Markus Fasel <M.Fasel@gsi.de>                                       //
  24 //                                                                        //
  25 ////////////////////////////////////////////////////////////////////////////
  26
  27 #include <TObjArray.h>
  28 #include <TClonesArray.h>
  29 #include <TPad.h>
  30 #include <TLegend.h>
  31 #include <TProfile.h>
  32 #include <TMath.h>
  33 #include <TDatabasePDG.h>
  34 #include "TTreeStream.h"
  35
  36 #include "AliMagF.h"
  37 #include "AliPID.h"
  38 #include "AliESDtrack.h"
  39 #include "AliMathBase.h"
  40 #include "AliTrackReference.h"
  41 #include "AliAnalysisManager.h"
  42
  43 #include "AliTRDcluster.h"
  44 #include "AliTRDseedV1.h"
  45 #include "AliTRDtrackV1.h"
  46 #include "AliTRDpidUtil.h"
  47 #include "Cal/AliTRDCalPID.h"
  48 #include "info/AliTRDtrackInfo.h"
  49 #include "AliTRDinfoGen.h"
  50 #include "AliTRDefficiencyMC.h"
  51
  52 ClassImp(AliTRDefficiencyMC)
  53 Float_t AliTRDefficiencyMC::fgPCut   = 0.2; //[GeV/c]
  54 Float_t AliTRDefficiencyMC::fgPhiCut = 50.; //[deg]
  55 Float_t AliTRDefficiencyMC::fgThtCut = 50.; //[deg]
  56 //_____________________________________________________________________________
  57 AliTRDefficiencyMC::AliTRDefficiencyMC()
  58   :AliTRDrecoTask()
  59 {
  60   //
  61   // Default constructor
  62   //
  63 }
  64
  65 AliTRDefficiencyMC::AliTRDefficiencyMC(char* name)
  66   :AliTRDrecoTask(name, "Combined Tracking Efficiency")
  67 {
  68   //
  69   // Default constructor
  70   //
  71 }
  72
  73 //_____________________________________________________________________________
  74 void AliTRDefficiencyMC::UserCreateOutputObjects(){
  75   //
  76   // Create output objects
  77   //
  78
  79   fContainer = Histos();
  80   PostData(1, fContainer);
  81 }
  82
  83 //_____________________________________________________________________________
  84 void AliTRDefficiencyMC::UserExec(Option_t *){
  85   //
  86   // Execute the task:
  87   //
  88   // Loop over TrackInfos
  89   // 1st: check if there is a trackTRD
  90   // 2nd: put conditions on the track:
  91   //      - check if we did not register it before
  92   //      - check if we have Track References for the track
  93   // 3rd: Register track:
  94   //      - accepted if both conditions are fulfilled
  95   //      - contamination if at least one condition is not fulfilled
  96   // 4th: check Monte-Carlo Track wheter findable or not if there is no TRD track in this track info
  97   // 5th: register MC track as rejected if findable and not jet registered
  98   // Match the registers accepted and rejected and clean register rejected
  99   // Fill the histograms
 100   //
 101   const Int_t kArraySize = 10000;     // Store indices of track references in an array
 102   Int_t indexAccept[kArraySize],
 103         indexReject[kArraySize],
 104         indexContam[kArraySize];
 105   memset(indexAccept, 0, sizeof(Int_t) * kArraySize);
 106   memset(indexReject, 0, sizeof(Int_t) * kArraySize);
 107   memset(indexContam, 0, sizeof(Int_t) * kArraySize);
 108   Int_t naccept(0),
 109         nreject(0),
 110         nfindnt(0),
 111         nkink(0),
 112         ncontam(0);
 113   Bool_t isContamination = kFALSE;
 114
 115   fTracks = dynamic_cast<TObjArray *>(GetInputData(1));
 116   if(!fTracks) return;
 117   Int_t nTrackInfos = fTracks->GetEntriesFast();
 118   AliDebug(2, Form("   CANDIDATE TRACKS[%d]", nTrackInfos));
 119
 120   //AliTRDtrackV1 *trackTRD(NULL);
 121   AliTRDtrackInfo *trkInf(NULL);
 122   for(Int_t itinf = 0; itinf < nTrackInfos; itinf++){
 123     trkInf = dynamic_cast<AliTRDtrackInfo *>(fTracks->UncheckedAt(itinf));
 124     if(!trkInf) continue;
 125
 126     if(trkInf->GetTrack() || trkInf->GetNumberOfClustersRefit()){
 127       isContamination = (IsRegistered(trkInf,indexAccept,naccept)>=0);
 128       if(!trkInf->GetNTrackRefs()){
 129         // We reject the track since the Monte Carlo Information is missing
 130         AliDebug(2, Form("MC(Track Reference) missing @ label[%d]", trkInf->GetLabel()));
 131         isContamination = kTRUE;
 132         // Debugging
 133         if(trkInf && DebugLevel()>5) FillStreamTrackWOMC(trkInf);
 134       }
 135       if(isContamination){
 136         // reject kink (we count these only once)
 137         if(trkInf->GetKinkIndex()){
 138           AliDebug(4, Form("  track @ idx[%d] MC[%d] is kink.", itinf, trkInf->GetLabel()));
 139           nkink++;
 140           continue;
 141         }
 142         // Register track as contamination
 143         AliDebug(4, Form("  track @ idx[%d] MC[%d] is contamination.", itinf, trkInf->GetLabel()));
 144         indexContam[ncontam++]=itinf;
 145         continue;
 146       }
 147       // Accept track
 148       AliDebug(4, Form("  track @ idx[%d] is ACCEPTED.", itinf));
 149       // Register track as accepted
 150       indexAccept[naccept++] = itinf;
 151     }else{
 152       Int_t code(0);
 153       if((code=IsFindableNot(trkInf))){
 154         AliDebug(4, Form("  track @ idx[%d] MC[%d] not findable [%d].", itinf, trkInf->GetLabel(), code));
 155         nfindnt++;
 156       } else {
 157         // register track as rejected if not already registered there
 158         // Attention:
 159         // These track infos are not!!! registered as contamination
 160         if(IsRegistered(trkInf, indexReject, nreject)<0){
 161           AliDebug(4, Form("  track @ idx[%d] MC[%d] is missed.", itinf, trkInf->GetLabel()));
 162           indexReject[nreject++] = itinf;
 163         }
 164       }
 165     }
 166   }
 167   AliDebug(2, Form("TRACKS STATISTICS naccept[%d] ncontam[%d] nkink[%d] nmissed[%d] nfindnt[%d] ALL[%d] LOST[%d]", naccept, ncontam, nkink, nreject, nfindnt, naccept+ncontam+nkink, nTrackInfos-(naccept+nreject+ncontam+nkink+nfindnt)));
 168
 169   // we have to check if the rejected tracks are registered as found
 170   // a possible source for this:
 171   // first the track is not found by the barrel tracking but it is later found
 172   // by the stand alone tracking, then two track info objects with the same
 173   // label would be created
 174   // Attention:
 175   // these tracks are not! registered as contamination
 176   Int_t tmprejected[kArraySize]; Int_t nrej = nreject;
 177   memcpy(tmprejected, indexReject, sizeof(Int_t) * nreject);
 178   nreject = 0;
 179   for(Int_t irej = 0; irej < nrej; irej++){
 180     if(!(trkInf = dynamic_cast<AliTRDtrackInfo *>(fTracks->At(tmprejected[irej])))) continue;
 181     Int_t idx(-1);
 182     if((idx = IsRegistered(trkInf,indexAccept,naccept))<0){
 183       indexReject[nreject++] = tmprejected[irej];
 184     }else{
 185       //printf("tracks @ accept[%d] missed[%d] are the same.\n", indexAccept[idx], tmprejected[irej]);
 186     }
 187   }
 188
 189   // Fill Histograms
 190   FillHistograms(naccept, &indexAccept[0], kAccept);
 191   FillHistograms(nreject, &indexReject[0], kMiss);
 192   FillHistograms(ncontam, &indexContam[0], kFake);
 193
 194   Int_t nall(naccept + nreject);
 195   AliInfo(Form("%3d Tracks: MC[%3d] TRD[%3d | %5.2f%%] Fake[%3d | %5.2f%%]",
 196     (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(),
 197     nall,
 198     naccept,
 199     (nall ? 1.E2*Float_t(naccept)/Float_t(nall) : 0.),
 200     ncontam,
 201     (nall ? 1.E2*Float_t(ncontam)/Float_t(nall) : 0.)));
 202 }
 203
 204
 205 //_____________________________________________________________________________
 206 Bool_t AliTRDefficiencyMC::PostProcess()
 207 {
 208   //
 209   // Post Process
 210   //
 211   // Change the histogram style
 212   // For species histograms apply the colors connected with the given particle species
 213   fNRefFigures = 8;
 214   return kTRUE;
 215 }
 216
 217 //_____________________________________________________________________________
 218 Bool_t AliTRDefficiencyMC::GetRefFigure(Int_t ifig){
 219   //
 220   // Plot the histograms
 221   //
 222   if(ifig >= fNRefFigures) return kFALSE;
 223   if(!gPad) return kFALSE;
 224   gPad->SetLogx(kTRUE);
 225   TH1 *h(NULL);
 226   if(ifig < 2){
 227     if(!(h = (dynamic_cast<TH1 *>(fContainer->At(ifig))))) return kFALSE;
 228     h->Draw("e1");
 229     return kTRUE;
 230   }
 231   TLegend *leg=new TLegend(.65, .12, .85, .3);
 232   leg->SetHeader("Charge");
 233   leg->SetBorderSize(1);leg->SetFillColor(kWhite);
 234   switch(ifig){
 235   case 2:
 236     if(!(h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram)))) return kFALSE;;
 237     h->Draw("e1"); leg->AddEntry(h, "  -", "pl");
 238     if(!(h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+1)))) return kFALSE;;
 239     h->Draw("e1same"); leg->AddEntry(h, "  +", "pl");
 240     leg->Draw();
 241     break;
 242   case 3:
 243     if(!(h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+2)))) return kFALSE;;
 244     h->Draw("e1"); leg->AddEntry(h, "  -", "pl");
 245     if(!(h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+3)))) return kFALSE;;
 246     h->Draw("e1same"); leg->AddEntry(h, "  +", "pl");
 247     leg->Draw();
 248     break;
 249   case 4:
 250     if(!(h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+4)))) return kFALSE;;
 251     h->Draw("e1"); leg->AddEntry(h, "  -", "pl");
 252     if(!(h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+5)))) return kFALSE;;
 253     h->Draw("e1same"); leg->AddEntry(h, "  +", "pl");
 254     leg->Draw();
 255     break;
 256   case 5:
 257     if(!(h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+6)))) return kFALSE;;
 258     h->Draw("e1"); leg->AddEntry(h, "  -", "pl");
 259     if(!(h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+7)))) return kFALSE;;
 260     h->Draw("e1same"); leg->AddEntry(h, "  +", "pl");
 261     leg->Draw();
 262     break;
 263   case 6:
 264     if(!(h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+8)))) return kFALSE;;
 265     h->Draw("e1"); leg->AddEntry(h, "  -", "pl");
 266     if(!(h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+9)))) return kFALSE;;
 267     h->Draw("e1same"); leg->AddEntry(h, "  +", "pl");
 268     leg->Draw();
 269     break;
 270   case 7:
 271     if(!(h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+10)))) return kFALSE;;
 272     h->Draw("e1"); leg->AddEntry(h, "  -", "pl");
 273     if(!(h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+11)))) return kFALSE;;
 274     h->Draw("e1same"); leg->AddEntry(h, "  +", "pl");
 275     leg->Draw();
 276     break;
 277   case 8:
 278     if(!(h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+12)))) return kFALSE;;
 279     h->Draw("e1"); leg->AddEntry(h, "  -", "pl");
 280     if(!(h=dynamic_cast<TH1 *>(fContainer->At(kEfficiencySpeciesHistogram+13)))) return kFALSE;;
 281     h->Draw("e1same"); leg->AddEntry(h, "  +", "pl");
 282     leg->Draw();
 283     break;
 284   }
 285   return kTRUE;
 286 }
 287
 288 //_____________________________________________________________________________
 289 TObjArray *AliTRDefficiencyMC::Histos(){
 290   //
 291   // Create the histograms
 292   //
 293
 294   if(fContainer) return fContainer;
 295   const Int_t nbins = AliTRDCalPID::kNMom;
 296   Float_t xbins[nbins+1] = {fgPCut, .7, .9, 1.3, 1.7, 2.4, 3.5, 4.5, 5.5, 7., 9., 11.};
 297   const Int_t marker[2][AliPID::kSPECIES+1] = {
 298     {20, 21, 22, 23, 29, 2},
 299     {24, 25, 26, 27, 30, 5}
 300   };
 301
 302   fContainer = new TObjArray();fContainer->Expand(14);
 303
 304   TH1 *h(NULL);
 305   fContainer->AddAt(h=new TProfile("hEff", "Tracking Efficiency ALL", nbins, xbins), kEfficiencyHistogram);
 306   h->SetMarkerStyle(22);
 307   h->SetMarkerColor(kBlue);
 308   h->GetXaxis()->SetTitle("p [GeV/c]");
 309   h->GetXaxis()->SetMoreLogLabels();
 310   h->GetYaxis()->SetTitle("Efficiency");
 311   h->GetYaxis()->SetRangeUser(0.2, 1.1);
 312   fContainer->AddAt(h=new TProfile("hFake", "Fake Tracks", nbins, xbins), kContaminationHistogram);
 313   h->SetMarkerStyle(22);
 314   h->SetMarkerColor(kBlue);
 315   h->GetXaxis()->SetTitle("p [GeV/c]");
 316   h->GetXaxis()->SetMoreLogLabels();
 317   h->GetYaxis()->SetTitle("Contamination");
 318
 319   Char_t sign[]={'+', '-'};
 320   for(Int_t isign = 0; isign < 2; isign++){
 321     for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
 322       fContainer->AddAt(h=new TProfile(
 323         Form("hEff_%s%c", AliPID::ParticleShortName(ispec), sign[isign]),
 324         Form("Tracking Efficiency for %s", AliPID::ParticleName(ispec)), nbins, xbins),
 325         kEfficiencySpeciesHistogram+ispec*2+isign);
 326       h->SetMarkerStyle(marker[isign][ispec]);
 327       h->SetLineColor(AliTRDCalPID::GetPartColor(ispec));
 328       h->SetMarkerColor(kBlack);
 329       h->GetXaxis()->SetTitle("p [GeV/c]");
 330       h->GetXaxis()->SetMoreLogLabels();
 331       h->GetYaxis()->SetTitle("Efficiency");
 332       h->GetYaxis()->SetRangeUser(0.2, 1.1);
 333     }
 334
 335     fContainer->AddAt(h=new TProfile(Form("hEff_PID%c", sign[isign]), "Tracking Efficiency no PID", nbins, xbins), kEfficiencySpeciesHistogram+AliPID::kSPECIES*2+isign);
 336     h->SetMarkerStyle(marker[isign][AliPID::kSPECIES]);
 337     h->SetMarkerColor(kBlack);h->SetLineColor(kBlack);
 338     h->GetXaxis()->SetTitle("p [GeV/c]");
 339     h->GetXaxis()->SetMoreLogLabels();
 340     h->GetYaxis()->SetTitle("Efficiency");
 341     h->GetYaxis()->SetRangeUser(0.2, 1.1);
 342   }
 343   return fContainer;
 344 }
 345
 346 //_____________________________________________________________________________
 347 Int_t AliTRDefficiencyMC::IsFindableNot(AliTRDtrackInfo * const trkInf){
 348   //
 349   // Apply Cuts on the Monte Carlo track references
 350   // return whether track is findable or not
 351   //
 352
 353
 354   const Float_t chmbHght = AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght();
 355   const Float_t eps(1.E-3);
 356   Int_t ntr(trkInf->GetNTrackRefs());
 357
 358   AliDebug(10, Form("  CANDIDATE TrackRefs[%d]", ntr));
 359   // Check if track is findable
 360   Double_t mom(0.), phi(0.), tht(0.);
 361   Float_t xmin = 10000.0, xmax = 0.0;
 362   Float_t ymin = 0.0, ymax = 0.0;
 363   Float_t zmin = 0.0, zmax = 0.0;
 364   Float_t lastx = 0.0, x = 0.0;
 365   Int_t nLayers(0), ntrTRD(0);
 366   Int_t sector[20];
 367   AliTrackReference *trackRef(NULL);
 368   for(Int_t itr(0); itr<ntr; itr++){
 369     if(!(trackRef = trkInf->GetTrackRef(itr))) continue;
 370     x = trackRef->LocalX();
 371     // Be Sure that we are inside TRD
 372     if(x < AliTRDinfoGen::GetEndTPC() || x > AliTRDinfoGen::GetEndTRD()) continue;
 373     sector[ntrTRD] = Int_t(trackRef->Alpha()/AliTRDgeometry::GetAlpha());
 374     AliDebug(10, Form("    [%2d] x[%7.2f] y[%7.2f] z[%7.2f] Sec[%2d]", itr, trackRef->LocalX(), trackRef->LocalY(), trackRef->Z(), sector[ntrTRD]));
 375     if(x < xmin){
 376       xmin = trackRef->LocalX();
 377       ymin = trackRef->LocalY();
 378       zmin = trackRef->Z();
 379       mom  = trackRef->P();
 380     } else if(x > xmax){
 381       xmax = trackRef->LocalX();
 382       ymax = trackRef->LocalY();
 383       zmax = trackRef->Z();
 384     }
 385     if(itr > 0){
 386       Float_t dist = TMath::Abs(x - lastx);
 387       if(TMath::Abs(dist - chmbHght) < eps && sector[ntrTRD]==sector[0]){
 388         AliDebug(10, Form("    dx = %7.2f", dist));
 389         nLayers++;
 390       }
 391     }
 392     lastx = x;
 393     ntrTRD++; if(ntrTRD>=20) break;
 394   }
 395   Double_t dx(xmax - xmin);
 396   if(TMath::Abs(dx)<eps) return kNoChmb;
 397
 398   phi = (ymax -ymin)/dx;
 399   tht = (zmax -zmin)/dx;
 400   phi=TMath::ATan(phi)*TMath::RadToDeg();
 401   tht=TMath::ATan(tht)*TMath::RadToDeg();
 402   Bool_t primary = trkInf->IsPrimary();
 403   const AliTRDtrackInfo::AliESDinfo *esd(trkInf->GetESDinfo());
 404   AliDebug(10, Form("    p=%6.3f[GeV/c] phi=%6.2f[deg] theta=%6.2f[deg] nLy[%d]",
 405       mom, phi, tht, nLayers));
 406   if(DebugLevel()){
 407     (*DebugStream()) << "IsFindable"
 408       << "P="       << mom
 409       << "Phi="     << phi
 410       << "Tht="     << tht
 411       << "Ntr="     << ntrTRD
 412       << "NLy="     << nLayers
 413       << "Primary=" << primary
 414       << "\n";
 415   }
 416
 417   // Apply cuts
 418   if(!nLayers) return kNoChmb;
 419   if(xmax < xmin) return kCurved;
 420   if(mom < fgPCut) return kPCut;
 421
 422
 423   if(TMath::Abs(phi) > fgPhiCut) return kPhiCut;
 424   if(TMath::Abs(tht) > fgThtCut) return kThtCut;
 425
 426   if(nLayers < 4){
 427     if(!esd)return kLayer;
 428     if(!(esd->GetStatus() & AliESDtrack::kTPCout)) return kLayer;
 429   }
 430
 431   //if(!trkInf->IsPrimary()) {failCode=kPrimary; return kFALSE;}
 432
 433   return kFindable;
 434 }
 435
 436 //_____________________________________________________________________________
 437 void AliTRDefficiencyMC::FillHistograms(Int_t nTracks, Int_t *indices, ETRDefficiencyMCstatus mode){
 438   //
 439   // Fill Histograms in three different modes:
 440   // 1st tracks which are found and accepted
 441   // 2nd tracks which are not found and not already counted
 442   // 3rd contaminating tracks: either double counts (kinks!) or tracks with no MC hit inside TRD
 443   //
 444
 445   TDatabasePDG *dbPDG(TDatabasePDG::Instance());
 446   Double_t trkmom(0.);   // the track momentum
 447   Int_t trkpdg(-1);      // particle PDG code
 448   AliTRDtrackInfo *trkInf(NULL);
 449   for(Int_t itk = 0; itk < nTracks; itk++){
 450     if(!(trkInf = dynamic_cast<AliTRDtrackInfo *>(fTracks->At(indices[itk])))) continue;
 451     if(trkInf->GetNTrackRefs()){
 452       // use Monte-Carlo Information for Momentum and PID
 453       trkmom = trkInf->GetTrackRef(0)->P();
 454       trkpdg = trkInf->GetPDG();
 455     }else{
 456       // Use TPC Momentum
 457       trkmom = trkInf->GetTrack()->P();
 458     }
 459
 460     TProfile *hp(NULL);
 461     const Char_t *cmode(NULL);
 462     switch(mode){
 463       case kAccept:
 464         if(!(hp = (dynamic_cast<TProfile *>(fContainer->At(kEfficiencyHistogram))))) continue;
 465         hp->Fill(trkmom, 1);
 466         if(!(hp = ((dynamic_cast<TProfile *>(fContainer->At(kContaminationHistogram)))))) continue;
 467         hp->Fill(trkmom, 0);
 468         cmode="ACCEPT";
 469         break;
 470       case kMiss:
 471         if(!(hp = ((dynamic_cast<TProfile *>(fContainer->At(kEfficiencyHistogram)))))) continue;
 472         hp->Fill(trkmom, 0);
 473         if(!(hp = ((dynamic_cast<TProfile *>(fContainer->At(kContaminationHistogram)))))) continue;
 474         hp->Fill(trkmom, 0);
 475         cmode="MISS";
 476         break;
 477       case kFake:
 478         if(!(hp = ((dynamic_cast<TProfile *>(fContainer->At(kContaminationHistogram)))))) continue;
 479         hp->Fill(trkmom, 1);
 480         cmode="FAKE";
 481         break;
 482     }
 483     AliDebug(3, Form(" track[%d] MC[%d] Mode[%s]", indices[itk], trkInf->GetLabel(), cmode));
 484
 485     // Fill species histogram
 486     Int_t idxSpec = AliTRDpidUtil::Pdg2Pid(TMath::Abs(trkpdg));
 487     Int_t sign = dbPDG->GetParticle(trkpdg)->Charge() > 0. ? 1 : 0;
 488     //printf("[%d]%s pdg[%d] sign[%d]\n", idxSpec, AliPID::ParticleName(idxSpec), trkpdg, sign);
 489     if(idxSpec < 0) idxSpec = AliPID::kSPECIES;
 490     if(!(hp = (dynamic_cast<TProfile *>(fContainer->At(kEfficiencySpeciesHistogram + idxSpec*2+sign))))) continue;
 491     hp->Fill(trkmom, mode==kAccept?1:0);
 492   }
 493 }
 494
 495 //_____________________________________________________________________________
 496 void AliTRDefficiencyMC::FillStreamTrackWOMC(AliTRDtrackInfo * const trkInf){
 497   // fill debug stream
 498   // we want to know:
 499   //  1. The event number
 500   //  2. The track label
 501   //  3. The TRD track label
 502   //  4. The frequency of the TRD Label
 503   //  5. Momentum from TPC (NO STAND ALONE TRACK)
 504   //  6. TPC Phi angle
 505   //  7. the TRD track
 506   //  8. Monte Carlo PID
 507   //  9. We check the Labels of the TRD track according to them we search the maching Monte-Carlo track.
 508   //     From the matching Monte-Carlo track we store trackRefs, phi and momentum
 509   // 10. We may also want to keep the kink index
 510   Double_t mom = trkInf->GetESDinfo()->GetOuterParam()->P();
 511   Int_t event = (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry();
 512   Int_t label = trkInf->GetLabel();
 513   Int_t kinkIndex = trkInf->GetKinkIndex();
 514   Int_t pdg = trkInf->GetPDG();
 515   Double_t phiTPC = trkInf->GetESDinfo()->GetOuterParam()->Phi();
 516   Int_t labelsTRD[180]; // Container for the cluster labels
 517   Int_t sortlabels[360];        // Cluster Labels sorted according their occurancy
 518   AliTRDseedV1 *tracklet(NULL);
 519   AliTRDcluster *c(NULL);
 520   Int_t nclusters(0);
 521   AliTRDtrackV1 *trackTRD = trkInf->GetTrack();
 522   for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
 523     tracklet = trackTRD->GetTracklet(il);
 524     if(!tracklet) continue;
 525     tracklet->ResetClusterIter();
 526     c = NULL;
 527     while((c = tracklet->NextCluster())) labelsTRD[nclusters++] = c->GetLabel(0);
 528   }
 529   // Determine Label and Frequency
 530   AliMathBase::Freq(nclusters, const_cast<const Int_t *>(&labelsTRD[0]), &sortlabels[0], kTRUE);
 531   Int_t labelTRD = sortlabels[0];
 532   Int_t freqTRD = sortlabels[1];
 533   // find the track info object matching to the TRD track
 534   AliTRDtrackInfo *realtrack = 0;
 535   TObjArrayIter rtiter(fTracks);
 536   while((realtrack = (AliTRDtrackInfo *)rtiter())){
 537     if(realtrack->GetLabel() != labelTRD) continue;
 538     break;
 539   }
 540   TClonesArray trackRefs("AliTrackReference");
 541   Int_t realPdg = -1;
 542   Double_t realP = 0.;
 543   Double_t realPhi = 0.;
 544   if(realtrack){
 545     // pack the track references into the trackRefsContainer
 546     for(Int_t iref = 0; iref < realtrack->GetNTrackRefs(); iref++){
 547     new(trackRefs[iref])AliTrackReference(*(realtrack->GetTrackRef(iref)));
 548     }
 549     realPdg = realtrack->GetPDG();
 550     if(realtrack->GetNTrackRefs()){
 551       realP = realtrack->GetTrackRef(0)->P();
 552       realPhi = realtrack->GetTrackRef(0)->Phi();
 553     }
 554   }
 555   (*DebugStream()) << "EffMCfake"
 556     << "Event=" << event
 557     << "Label=" << label
 558     << "labelTRD=" << labelTRD
 559     << "FreqTRDlabel=" << freqTRD
 560     << "TPCp="  << mom
 561     << "phiTPC=" << phiTPC
 562     << "trackTRD=" << trackTRD
 563     << "PDG="   << pdg
 564     << "TrackRefs=" << &trackRefs
 565     << "RealPDG=" << realPdg
 566     << "RealP=" << realP
 567     << "RealPhi" << realPhi
 568     << "KinkIndex=" << kinkIndex
 569     << "\n";
 570 }
 571
 572 //_____________________________________________________________________________
 573 Int_t AliTRDefficiencyMC::IsRegistered(const AliTRDtrackInfo * const trkInf, Int_t *indices, Int_t nTracks){
 574   //
 575   // Checks if track is registered in a given mode
 576   //
 577
 578   AliTRDtrackInfo *ti(NULL);
 579   Int_t label(trkInf->GetLabel());
 580   for(Int_t il(nTracks); il--;){
 581     if(!(ti = dynamic_cast<AliTRDtrackInfo *>(fTracks->At(indices[il])))) continue;
 582     if(ti->GetLabel() == label) return il;
 583   }
 584   return -1;
 585 }
 586