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