PWG1/TRD/AliTRDinfoGen.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: AliTRDinfoGen.cxx 27496 2008-07-22 08:35:45Z cblume $ */
  17
  18 ////////////////////////////////////////////////////////////////////////////
  19 //
  20 //  Tender wagon for TRD performance/calibration train
  21 //
  22 // In this wagon the information from
  23 //   - ESD
  24 //   - Friends [if available]
  25 //   - MC [if available]
  26 // are grouped into AliTRDtrackInfo objects and fed to worker tasks
  27 //
  28 //  Authors:
  29 //    Markus Fasel <M.Fasel@gsi.de>
  30 //    Alexandru Bercuci <A.Bercuci@gsi.de>
  31 //
  32 ////////////////////////////////////////////////////////////////////////////
  33
  34 #include <TClonesArray.h>
  35 #include <TObjArray.h>
  36 #include <TObject.h>
  37 #include <TString.h>
  38 #include <TH1S.h>
  39 #include <TPad.h>
  40 #include <TFile.h>
  41 #include <TTree.h>
  42 #include <TROOT.h>
  43 #include <TChain.h>
  44 #include <TParticle.h>
  45
  46 #include "AliLog.h"
  47 #include "AliAnalysisManager.h"
  48 #include "AliESDEvent.h"
  49 #include "AliMCEvent.h"
  50 #include "AliESDInputHandler.h"
  51 #include "AliMCEventHandler.h"
  52
  53 #include "AliESDfriend.h"
  54 #include "AliESDfriendTrack.h"
  55 #include "AliESDHeader.h"
  56 #include "AliESDtrack.h"
  57 #include "AliESDv0.h"
  58 #include "AliESDtrackCuts.h"
  59 #include "AliMCParticle.h"
  60 #include "AliPID.h"
  61 #include "AliStack.h"
  62 #include "AliTRDtrackV1.h"
  63 #include "AliTrackReference.h"
  64 #include "AliTRDgeometry.h"
  65 #include "AliTRDcluster.h"
  66 #include "AliTRDseedV1.h"
  67 #include "TTreeStream.h"
  68
  69 #include <cstdio>
  70 #include <climits>
  71 #include <cstring>
  72 #include <iostream>
  73
  74 #include "AliTRDinfoGen.h"
  75 #include "info/AliTRDtrackInfo.h"
  76 #include "info/AliTRDeventInfo.h"
  77 #include "info/AliTRDv0Info.h"
  78 #include "info/AliTRDeventCuts.h"
  79 #include "macros/AliTRDperformanceTrain.h"
  80
  81 ClassImp(AliTRDinfoGen)
  82
  83 const Float_t AliTRDinfoGen::fgkITS = 100.; // to be checked
  84 const Float_t AliTRDinfoGen::fgkTPC = 290.;
  85 const Float_t AliTRDinfoGen::fgkTRD = 365.;
  86
  87 const Float_t AliTRDinfoGen::fgkEvVertexZ = 15.;
  88 const Int_t   AliTRDinfoGen::fgkEvVertexN = 1;
  89
  90 const Float_t AliTRDinfoGen::fgkTrkDCAxy  = 3.;
  91 const Float_t AliTRDinfoGen::fgkTrkDCAz   = 10.;
  92 const Int_t   AliTRDinfoGen::fgkNclTPC    = 70;
  93 const Float_t AliTRDinfoGen::fgkPt        = 0.2;
  94 const Float_t AliTRDinfoGen::fgkEta       = 0.9;
  95
  96 //____________________________________________________________________
  97 AliTRDinfoGen::AliTRDinfoGen()
  98   :AliAnalysisTaskSE()
  99   ,fEvTrigger(NULL)
 100   ,fESDev(NULL)
 101   ,fMCev(NULL)
 102   ,fEventCut(NULL)
 103   ,fTrackCut(NULL)
 104   ,fV0Cut(NULL)
 105   ,fTrackInfo(NULL)
 106   ,fEventInfo(NULL)
 107   ,fV0Info(NULL)
 108   ,fTracksBarrel(NULL)
 109   ,fTracksSA(NULL)
 110   ,fTracksKink(NULL)
 111   ,fV0List(NULL)
 112   ,fContainer(NULL)
 113   ,fDebugStream(NULL)
 114 {
 115   //
 116   // Default constructor
 117   //
 118   SetNameTitle("TRDinfoGen", "MC-REC TRD-track list generator");
 119 }
 120
 121 //____________________________________________________________________
 122 AliTRDinfoGen::AliTRDinfoGen(char* name)
 123   :AliAnalysisTaskSE(name)
 124   ,fEvTrigger(NULL)
 125   ,fESDev(NULL)
 126   ,fMCev(NULL)
 127   ,fEventCut(NULL)
 128   ,fTrackCut(NULL)
 129   ,fV0Cut(NULL)
 130   ,fTrackInfo(NULL)
 131   ,fEventInfo(NULL)
 132   ,fV0Info(NULL)
 133   ,fTracksBarrel(NULL)
 134   ,fTracksSA(NULL)
 135   ,fTracksKink(NULL)
 136   ,fV0List(NULL)
 137   ,fContainer(NULL)
 138   ,fDebugStream(NULL)
 139 {
 140   //
 141   // Default constructor
 142   //
 143   SetTitle("MC-REC TRD-track list generator");
 144   DefineOutput(kTracksBarrel, TObjArray::Class());
 145   DefineOutput(kTracksSA, TObjArray::Class());
 146   DefineOutput(kTracksKink, TObjArray::Class());
 147   DefineOutput(kEventInfo, AliTRDeventInfo::Class());
 148   DefineOutput(kV0List, TObjArray::Class());
 149   DefineOutput(kMonitor, TObjArray::Class()); // histogram list
 150 }
 151
 152 //____________________________________________________________________
 153 AliTRDinfoGen::~AliTRDinfoGen()
 154 {
 155 // Destructor
 156
 157   if(fDebugStream) delete fDebugStream;
 158   if(fEvTrigger) delete fEvTrigger;
 159   if(fV0Cut) delete fV0Cut;
 160   if(fTrackCut) delete fTrackCut;
 161   if(fEventCut) delete fEventCut;
 162   if(fTrackInfo) delete fTrackInfo; fTrackInfo = NULL;
 163   if(fEventInfo) delete fEventInfo; fEventInfo = NULL;
 164   if(fV0Info) delete fV0Info; fV0Info = NULL;
 165   if(fTracksBarrel){
 166     fTracksBarrel->Delete(); delete fTracksBarrel;
 167     fTracksBarrel = NULL;
 168   }
 169   if(fTracksSA){
 170     fTracksSA->Delete(); delete fTracksSA;
 171     fTracksSA = NULL;
 172   }
 173   if(fTracksKink){
 174     fTracksKink->Delete(); delete fTracksKink;
 175     fTracksKink = NULL;
 176   }
 177   if(fV0List){
 178     fV0List->Delete();
 179     delete fV0List;
 180     fV0List = NULL;
 181   }
 182   if(fContainer){
 183     fContainer->Delete();
 184     delete fContainer;
 185     fContainer = NULL;
 186   }
 187 }
 188
 189 //____________________________________________________________________
 190 Bool_t AliTRDinfoGen::GetRefFigure(Int_t)
 191 {
 192   if(!gPad){
 193     AliWarning("Please provide a canvas to draw results.");
 194     return kFALSE;
 195   }
 196   fContainer->At(0)->Draw("bar");
 197   return kTRUE;
 198 }
 199
 200 //____________________________________________________________________
 201 void AliTRDinfoGen::UserCreateOutputObjects()
 202 {
 203   //
 204   // Create Output Containers (TObjectArray containing 1D histograms)
 205   //
 206
 207   fTrackInfo = new AliTRDtrackInfo();
 208   fEventInfo = new AliTRDeventInfo();
 209   fV0Info    = new AliTRDv0Info();
 210   fTracksBarrel = new TObjArray(200); fTracksBarrel->SetOwner(kTRUE);
 211   fTracksSA = new TObjArray(20); fTracksSA->SetOwner(kTRUE);
 212   fTracksKink = new TObjArray(20); fTracksKink->SetOwner(kTRUE);
 213   fV0List = new TObjArray(10); fV0List->SetOwner(kTRUE);
 214
 215   // define general monitor
 216   fContainer = new TObjArray(1); fContainer->SetOwner(kTRUE);
 217   TH1 *h=new TH1S("hStat", "Run statistics;Observable;Entries", 12, -0.5, 11.5);
 218   TAxis *ax(h->GetXaxis());
 219   ax->SetBinLabel( 1, "ESD");
 220   ax->SetBinLabel( 2, "MC");
 221   ax->SetBinLabel( 3, "V0");
 222   ax->SetBinLabel( 4, "TPC");
 223   ax->SetBinLabel( 5, "TRDin");
 224   ax->SetBinLabel( 6, "TRDout");
 225   ax->SetBinLabel( 7, "Barrel");
 226   ax->SetBinLabel( 8, "BarrelMC");
 227   ax->SetBinLabel( 9, "SA");
 228   ax->SetBinLabel(10, "SAMC");
 229   ax->SetBinLabel(11, "Kink");
 230   ax->SetBinLabel(12, "KinkMC");
 231   fContainer->AddAt(h, 0);
 232   PostData(kMonitor, fContainer);
 233 }
 234
 235 //____________________________________________________________________
 236 Bool_t AliTRDinfoGen::Load(const Char_t *file, const Char_t *dir, const Char_t *name)
 237 {
 238 // Load data from performance file
 239
 240   if(!TFile::Open(file)){
 241     AliWarning(Form("Couldn't open file %s.", file));
 242     return kFALSE;
 243   }
 244   if(dir){
 245     if(!gFile->cd(dir)){
 246       AliWarning(Form("Couldn't cd to %s in %s.", dir, file));
 247       return kFALSE;
 248     }
 249   }
 250   TObjArray *o(NULL);
 251   const Char_t *tn=(name ? name : GetName());
 252   if(!(o = (TObjArray*)gDirectory->Get(tn))){
 253     AliWarning(Form("Missing histogram container %s.", tn));
 254     return kFALSE;
 255   }
 256   fContainer = (TObjArray*)o->Clone(GetName());
 257   gFile->Close();
 258   return kTRUE;
 259 }
 260
 261 //____________________________________________________________________
 262 void AliTRDinfoGen::UserExec(Option_t *){
 263   //
 264   // Run the Analysis
 265   //
 266
 267   fESDev = dynamic_cast<AliESDEvent*>(InputEvent());
 268   fMCev = MCEvent();
 269
 270   if(!fESDev){
 271     AliError("Failed retrieving ESD event");
 272     return;
 273   }
 274
 275   // event selection : trigger cut
 276   if(UseLocalEvSelection() && fEvTrigger){
 277     Bool_t kTRIGGERED(kFALSE);
 278     const TObjArray *trig = fEvTrigger->Tokenize(" ");
 279     for(Int_t itrig=trig->GetEntriesFast(); itrig--;){
 280       const Char_t *trigClass(((TObjString*)(*trig)[itrig])->GetName());
 281       if(fESDev->IsTriggerClassFired(trigClass)) {
 282         AliDebug(2, Form("Ev[%4d] Trigger[%s]", fESDev->GetEventNumberInFile(), trigClass));
 283         kTRIGGERED = kTRUE;
 284         break;
 285       }
 286     }
 287     if(!kTRIGGERED){
 288       AliDebug(2, Form("Reject Ev[%4d] Trigger", fESDev->GetEventNumberInFile()));
 289       return;
 290     }
 291     // select only physical events
 292     if(fESDev->GetEventType() != 7){
 293       AliDebug(2, Form("Reject Ev[%4d] EvType[%d]", fESDev->GetEventNumberInFile(), fESDev->GetEventType()));
 294       return;
 295     }
 296   }
 297
 298   // if the required trigger is a collision trigger then apply event vertex cut
 299   if(UseLocalEvSelection() && IsCollision()){
 300     const AliESDVertex *vertex = fESDev->GetPrimaryVertex();
 301     if(TMath::Abs(vertex->GetZv())<1.e-10 ||
 302        TMath::Abs(vertex->GetZv())>fgkEvVertexZ ||
 303        vertex->GetNContributors()<fgkEvVertexN) {
 304       AliDebug(2, Form("Reject Ev[%4d] Vertex Zv[%f] Nv[%d]", fESDev->GetEventNumberInFile(), TMath::Abs(vertex->GetZv()), vertex->GetNContributors()));
 305       return;
 306     }
 307   }
 308
 309   if(fEventCut && !fEventCut->IsSelected(fESDev, IsCollision())) return;
 310
 311   if(!fESDfriend){
 312     AliError("Failed retrieving ESD friend event");
 313     return;
 314   }
 315   if(HasMCdata() && !fMCev){
 316     AliError("Failed retrieving MC event");
 317     return;
 318   }
 319
 320   fTracksBarrel->Delete();
 321   fTracksSA->Delete();
 322   fTracksKink->Delete();
 323   fV0List->Delete();
 324   fEventInfo->Delete("");
 325   new(fEventInfo)AliTRDeventInfo(fESDev->GetHeader(), const_cast<AliESDRun *>(fESDev->GetESDRun()));
 326
 327   Bool_t *trackMap(NULL);
 328   AliStack * mStack(NULL);
 329   Int_t nTracksMC = HasMCdata() ? fMCev->GetNumberOfTracks() : 0, nTracksESD = fESDev->GetNumberOfTracks();
 330   if(HasMCdata()){
 331     mStack = fMCev->Stack();
 332     if(!mStack){
 333       AliError("Failed retrieving MC Stack");
 334       return;
 335     }
 336     trackMap = new Bool_t[nTracksMC];
 337     memset(trackMap, 0, sizeof(Bool_t) * nTracksMC);
 338   }
 339
 340   Double32_t dedx[100]; Int_t nSlices(0);
 341   Int_t nTRDout(0), nTRDin(0), nTPC(0)
 342        ,nclsTrklt
 343        ,nBarrel(0), nSA(0), nKink(0)
 344        ,nBarrelMC(0), nSAMC(0), nKinkMC(0);
 345   AliESDtrack *esdTrack = NULL;
 346   AliESDfriendTrack *esdFriendTrack = NULL;
 347   TObject *calObject = NULL;
 348   AliTRDtrackV1 *track = NULL;
 349   AliTRDseedV1 *tracklet = NULL;
 350   AliTRDcluster *cl = NULL;
 351   // LOOP 0 - over ESD tracks
 352   for(Int_t itrk = 0; itrk < nTracksESD; itrk++){
 353     new(fTrackInfo) AliTRDtrackInfo();
 354     esdTrack = fESDev->GetTrack(itrk);
 355     AliDebug(3, Form("\n%3d ITS[%d] TPC[%d] TRD[%d]\n", itrk, esdTrack->GetNcls(0), esdTrack->GetNcls(1), esdTrack->GetNcls(2)));
 356
 357     if(esdTrack->GetStatus()&AliESDtrack::kTPCout) nTPC++;
 358     if(esdTrack->GetStatus()&AliESDtrack::kTRDout) nTRDout++;
 359     if(esdTrack->GetStatus()&AliESDtrack::kTRDin) nTRDin++;
 360
 361     // look at external track param
 362     const AliExternalTrackParam *op = esdTrack->GetOuterParam();
 363     Double_t xyz[3];
 364     if(op){
 365       op->GetXYZ(xyz);
 366       op->Global2LocalPosition(xyz, op->GetAlpha());
 367       AliDebug(3, Form("op @ X[%7.3f]\n", xyz[0]));
 368     }
 369
 370     // read MC info
 371     Int_t fPdg = -1;
 372     Int_t label = -1; UInt_t alab=UINT_MAX;
 373     if(HasMCdata()){
 374       label = esdTrack->GetLabel();
 375       alab = TMath::Abs(label);
 376       // register the track
 377       if(alab < UInt_t(nTracksMC)){
 378         trackMap[alab] = kTRUE;
 379       } else {
 380         AliError(Form("MC label[%d] outside scope for Ev[%d] Trk[%d].", label, (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), itrk));
 381         continue;
 382       }
 383       AliMCParticle *mcParticle = NULL;
 384       if(!(mcParticle = (AliMCParticle*) fMCev->GetTrack(alab))){
 385         AliError(Form("MC particle label[%d] missing for Ev[%d] Trk[%d].", label, (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), itrk));
 386         continue;
 387       }
 388       fPdg = mcParticle->Particle()->GetPdgCode();
 389       Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
 390       Int_t iref = 0; AliTrackReference *ref = NULL;
 391       while(iref<nRefs){
 392         ref = mcParticle->GetTrackReference(iref);
 393         if(ref->LocalX() > fgkTPC) break;
 394         iref++;
 395       }
 396
 397       fTrackInfo->SetMC();
 398       fTrackInfo->SetPDG(fPdg);
 399       fTrackInfo->SetPrimary(mcParticle->Particle()->IsPrimary());
 400       fTrackInfo->SetLabel(label);
 401       Int_t jref = iref;//, kref = 0;
 402       while(jref<nRefs){
 403         ref = mcParticle->GetTrackReference(jref);
 404         if(ref->LocalX() > fgkTRD) break;
 405         AliDebug(4, Form("  trackRef[%2d (%2d)] @ %7.3f OK", jref-iref, jref, ref->LocalX()));
 406         fTrackInfo->AddTrackRef(ref);
 407         jref++;
 408       }
 409       AliDebug(3, Form("NtrackRefs[%d(%d)]", fTrackInfo->GetNTrackRefs(), nRefs));
 410     }
 411
 412     // copy some relevant info to TRD track info
 413     fTrackInfo->SetStatus(esdTrack->GetStatus());
 414     fTrackInfo->SetTrackId(esdTrack->GetID());
 415     Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);
 416     fTrackInfo->SetESDpid(p);
 417     fTrackInfo->SetESDpidQuality(esdTrack->GetTRDntrackletsPID());
 418     if(!nSlices) nSlices = esdTrack->GetNumberOfTRDslices();
 419     memset(dedx, 0, 100*sizeof(Double32_t));
 420     Int_t in(0);
 421     for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++)
 422       for(Int_t is=0; is<nSlices; is++)
 423         dedx[in++]=esdTrack->GetTRDslice(il, is);
 424     for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++) dedx[in++]=esdTrack->GetTRDmomentum(il);
 425     fTrackInfo->SetSlices(in, dedx);
 426     fTrackInfo->SetNumberOfClustersRefit(esdTrack->GetNcls(2));
 427     // some other Informations which we may wish to store in order to find problematic cases
 428     fTrackInfo->SetKinkIndex(esdTrack->GetKinkIndex(0));
 429     fTrackInfo->SetTPCncls(static_cast<UShort_t>(esdTrack->GetNcls(1)));
 430     nclsTrklt = 0;
 431
 432
 433     // read REC info
 434     esdFriendTrack = fESDfriend->GetTrack(itrk);
 435     if(esdFriendTrack){
 436       Int_t icalib = 0;
 437       while((calObject = esdFriendTrack->GetCalibObject(icalib++))){
 438         if(strcmp(calObject->IsA()->GetName(),"AliTRDtrackV1") != 0) continue; // Look for the TRDtrack
 439         if(!(track = dynamic_cast<AliTRDtrackV1*>(calObject))) break;
 440         AliDebug(4, Form("TRD track OK"));
 441         // Set the clusters to unused
 442         for(Int_t ipl = 0; ipl < AliTRDgeometry::kNlayer; ipl++){
 443           if(!(tracklet = track->GetTracklet(ipl))) continue;
 444           tracklet->ResetClusterIter();
 445           while((cl = tracklet->NextCluster())) cl->Use(0);
 446         }
 447         fTrackInfo->SetTrack(track);
 448         break;
 449       }
 450       AliDebug(3, Form("Ntracklets[%d]\n", fTrackInfo->GetNTracklets()));
 451     } else AliDebug(3, "No ESD friends");
 452     if(op) fTrackInfo->SetOuterParam(op);
 453
 454     if(DebugLevel() >= 1){
 455       AliTRDtrackInfo info(*fTrackInfo);
 456       (*DebugStream()) << "trackInfo"
 457       << "TrackInfo.=" << &info
 458       << "\n";
 459       info.Delete("");
 460     }
 461
 462     ULong_t status(esdTrack->GetStatus());
 463     if((status&AliESDtrack::kTPCout)){
 464       if(!esdTrack->GetKinkIndex(0)){ // Barrel  Track Selection
 465         Bool_t selected(kTRUE);
 466         if(UseLocalTrkSelection()){
 467           if(esdTrack->Pt() < fgkPt){
 468             AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] Pt[%5.2f]", fESDev->GetEventNumberInFile(), itrk, esdTrack->Pt()));
 469             selected = kFALSE;
 470           }
 471           if(selected && TMath::Abs(esdTrack->Eta()) > fgkEta){
 472             AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] Eta[%5.2f]", fESDev->GetEventNumberInFile(), itrk, TMath::Abs(esdTrack->Eta())));
 473             selected = kFALSE;
 474           }
 475           if(selected && esdTrack->GetTPCNcls() < fgkNclTPC){
 476             AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] NclTPC[%d]", fESDev->GetEventNumberInFile(), itrk, esdTrack->GetTPCNcls()));
 477             selected = kFALSE;
 478           }
 479           Float_t par[2], cov[3];
 480           esdTrack->GetImpactParameters(par, cov);
 481           if(IsCollision()){ // cuts on DCA
 482             if(selected && TMath::Abs(par[0]) > fgkTrkDCAxy){
 483               AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] DCAxy[%f]", fESDev->GetEventNumberInFile(), itrk, TMath::Abs(par[0])));
 484               selected = kFALSE;
 485             }
 486             if(selected && TMath::Abs(par[1]) > fgkTrkDCAz){
 487               AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] DCAz[%f]", fESDev->GetEventNumberInFile(), itrk, TMath::Abs(par[1])));
 488               selected = kFALSE;
 489             }
 490           } else if(selected && fMCev && !fMCev->IsPhysicalPrimary(alab)){;
 491             AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] Primary", fESDev->GetEventNumberInFile(), itrk));
 492             selected = kFALSE;
 493           }
 494         }
 495         if(fTrackCut && !fTrackCut->IsSelected(esdTrack)) selected = kFALSE;
 496         if(selected){
 497           fTracksBarrel->Add(new AliTRDtrackInfo(*fTrackInfo));
 498           nBarrel++;
 499         }
 500       } else {
 501         fTracksKink->Add(new AliTRDtrackInfo(*fTrackInfo));
 502         nKink++;
 503       }
 504     } else if((status&AliESDtrack::kTRDout) && !(status&AliESDtrack::kTRDin)){
 505       fTracksSA->Add(new AliTRDtrackInfo(*fTrackInfo));
 506       nSA++;
 507     }
 508     fTrackInfo->Delete("");
 509   }
 510
 511   // LOOP 1 - over ESD v0s
 512   Bool_t kFOUND(kFALSE);
 513   Float_t bField(fESDev->GetMagneticField());
 514   AliESDv0 *v0(NULL);
 515   for(Int_t iv0(0); iv0<fESDev->GetNumberOfV0s(); iv0++){
 516     if(!(v0 = fESDev->GetV0(iv0))) continue;
 517
 518     // purge V0 irelevant for TRD
 519     Int_t jb(0); AliTRDtrackInfo *ti(NULL);
 520     for(Int_t ib(0); ib<nBarrel; ib++){
 521       ti =   (AliTRDtrackInfo*)fTracksBarrel->At(ib);
 522       Int_t id(ti->GetTrackId());
 523       if(id!=v0->GetPindex() && id!=v0->GetNindex()) continue;
 524       kFOUND=kTRUE; ti->SetV0(); jb=ib+1;
 525       break;
 526     }
 527     if(!kFOUND) continue;
 528
 529     // register v0
 530     if(fV0Cut) new(fV0Info) AliTRDv0Info(*fV0Cut);
 531     else new(fV0Info) AliTRDv0Info();
 532     fV0Info->SetMagField(bField);
 533     fV0Info->SetV0tracks(fESDev->GetTrack(v0->GetPindex()), fESDev->GetTrack(v0->GetNindex()));
 534     fV0Info->SetV0Info(v0);
 535     // mark the other leg too
 536     for(Int_t ib(jb); ib<nBarrel; ib++){
 537       ti =   (AliTRDtrackInfo*)fTracksBarrel->At(ib);
 538       if(!fV0Info->HasTrack(ti)) continue;
 539       ti->SetV0();
 540       break;
 541     }
 542     //fV0Info->Print("a");
 543     fV0List->Add(new AliTRDv0Info(*fV0Info)); kFOUND=kFALSE;
 544   }
 545
 546   // LOOP 2 - over MC tracks which are passing TRD where the track is not reconstructed
 547   if(HasMCdata()){
 548     AliDebug(10, "Output of the MC track map:");
 549     for(Int_t itk = 0; itk < nTracksMC;  itk++) AliDebug(10, Form("trackMap[%d] = %s", itk, trackMap[itk] == kTRUE ? "TRUE" : "kFALSE"));
 550
 551     for(Int_t itk = 0; itk < nTracksMC; itk++){
 552       if(trackMap[itk]) continue;
 553       AliMCParticle *mcParticle =  (AliMCParticle*) fMCev->GetTrack(TMath::Abs(itk));
 554       Int_t fPdg = mcParticle->Particle()->GetPdgCode();
 555       Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
 556       Int_t iref = 0; AliTrackReference *ref = NULL;
 557       Int_t nRefsTRD = 0;
 558       new(fTrackInfo) AliTRDtrackInfo();
 559       fTrackInfo->SetMC();
 560       fTrackInfo->SetPDG(fPdg);
 561       while(iref<nRefs){ // count TRD TR
 562         Bool_t kIN(kFALSE);
 563         ref = mcParticle->GetTrackReference(iref);
 564         if(ref->LocalX() > fgkTPC && ref->LocalX() < fgkTRD){
 565           fTrackInfo->AddTrackRef(ref);
 566           nRefsTRD++;kIN=kTRUE;
 567         }
 568         AliDebug(4, Form("  trackRef[%2d] @ x[%7.3f] %s", iref, ref->LocalX(), kIN?"IN":"OUT"));
 569         iref++;
 570       }
 571       if(!nRefsTRD){
 572         // In this stage we at least require 1 hit inside TRD. What will be done with this tracks is a task for the
 573         // analysis job
 574         fTrackInfo->Delete("");
 575         continue;
 576       }
 577       fTrackInfo->SetPrimary(mcParticle->Particle()->IsPrimary());
 578       fTrackInfo->SetLabel(itk);
 579       if(DebugLevel() >= 1){
 580         AliTRDtrackInfo info(*fTrackInfo);
 581         (*DebugStream()) << "trackInfo"
 582         << "TrackInfo.=" << &info
 583         << "\n";
 584         info.Delete("");
 585       }
 586       AliDebug(3, Form("Add MC track @ label[%d] nTRDrefs[%d].", itk, nRefsTRD));
 587       // check where the track starts
 588       ref = mcParticle->GetTrackReference(0);
 589       if(ref->LocalX() < fgkITS){
 590         fTracksBarrel->Add(new AliTRDtrackInfo(*fTrackInfo));
 591         nBarrelMC++;
 592       } else if(ref->LocalX() < fgkTPC) {
 593         fTracksKink->Add(new AliTRDtrackInfo(*fTrackInfo));
 594         nKinkMC++;
 595       } else if(nRefsTRD>6){
 596         fTracksSA->Add(new AliTRDtrackInfo(*fTrackInfo));
 597         nSAMC++;
 598       }
 599       fTrackInfo->Delete("");
 600     }
 601     delete[] trackMap;
 602   }
 603   AliDebug(1, Form(
 604     "\nEv[%3d] Tracks: ESD[%d] MC[%d] V0[%d]\n"
 605     "        TPCout[%d] TRDin[%d] TRDout[%d]\n"
 606     "        Barrel[%3d+%3d=%3d] SA[%2d+%2d=%2d] Kink[%2d+%2d=%2d]"
 607     ,(Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nTracksESD, nTracksMC, fV0List->GetEntries()
 608     , nTPC, nTRDin, nTRDout
 609     ,nBarrel, nBarrelMC, fTracksBarrel->GetEntries()
 610     ,nSA, nSAMC, fTracksSA->GetEntries()
 611     ,nKink, nKinkMC, fTracksKink->GetEntries()
 612   ));
 613   // save track statistics
 614   TH1 *h((TH1S*)fContainer->At(0));
 615   h->Fill( 0., nTracksESD);
 616   h->Fill( 1., nTracksMC);
 617   h->Fill( 2., fV0List->GetEntries());
 618   h->Fill( 3., nTPC);
 619   h->Fill( 4., nTRDin);
 620   h->Fill( 5., nTRDout);
 621   h->Fill( 6., nBarrel);
 622   h->Fill( 7., nBarrelMC);
 623   h->Fill( 8., nSA);
 624   h->Fill( 9., nSAMC);
 625   h->Fill(10., nKink);
 626   h->Fill(11., nKinkMC);
 627
 628   PostData(kTracksBarrel, fTracksBarrel);
 629   PostData(kTracksSA, fTracksSA);
 630   PostData(kTracksKink, fTracksKink);
 631   PostData(kEventInfo, fEventInfo);
 632   PostData(kV0List, fV0List);
 633 }
 634
 635 //____________________________________________________________________
 636 void AliTRDinfoGen::SetLocalV0Selection(AliTRDv0Info *v0)
 637 {
 638 // Set V0 cuts from outside
 639
 640   if(!fV0Cut) fV0Cut = new AliTRDv0Info(*v0);
 641   else new(fV0Cut) AliTRDv0Info(*v0);
 642 }
 643
 644 //____________________________________________________________________
 645 void AliTRDinfoGen::SetTrigger(const Char_t *trigger)
 646 {
 647   if(!fEvTrigger) fEvTrigger = new TString(trigger);
 648   else (*fEvTrigger) = trigger;
 649 }
 650
 651 //____________________________________________________________________
 652 TTreeSRedirector* AliTRDinfoGen::DebugStream()
 653 {
 654   if(!fDebugStream){
 655     TDirectory *savedir = gDirectory;
 656     fDebugStream = new TTreeSRedirector("TRD.DebugInfoGen.root");
 657     savedir->cd();
 658   }
 659   return fDebugStream;
 660 }
 661
 662
 663