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("infoGen", "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 }
 233
 234 //____________________________________________________________________
 235 Bool_t AliTRDinfoGen::Load(const Char_t *file, const Char_t *dir, const Char_t *name)
 236 {
 237 // Load data from performance file
 238
 239   if(!TFile::Open(file)){
 240     AliWarning(Form("Couldn't open file %s.", file));
 241     return kFALSE;
 242   }
 243   if(dir){
 244     if(!gFile->cd(dir)){
 245       AliWarning(Form("Couldn't cd to %s in %s.", dir, file));
 246       return kFALSE;
 247     }
 248   }
 249   TObjArray *o(NULL);
 250   const Char_t *tn=(name ? name : GetName());
 251   if(!(o = (TObjArray*)gDirectory->Get(tn))){
 252     AliWarning(Form("Missing histogram container %s.", tn));
 253     return kFALSE;
 254   }
 255   fContainer = (TObjArray*)o->Clone(GetName());
 256   gFile->Close();
 257   return kTRUE;
 258 }
 259
 260 //____________________________________________________________________
 261 void AliTRDinfoGen::UserExec(Option_t *){
 262   //
 263   // Run the Analysis
 264   //
 265
 266   fESDev = dynamic_cast<AliESDEvent*>(InputEvent());
 267   fMCev = MCEvent();
 268
 269   if(!fESDev){
 270     AliError("Failed retrieving ESD event");
 271     return;
 272   }
 273
 274   // event selection : trigger cut
 275   if(UseLocalEvSelection() && fEvTrigger){
 276     Bool_t kTRIGGERED(kFALSE);
 277     const TObjArray *trig = fEvTrigger->Tokenize(" ");
 278     for(Int_t itrig=trig->GetEntriesFast(); itrig--;){
 279       const Char_t *trigClass(((TObjString*)(*trig)[itrig])->GetName());
 280       if(fESDev->IsTriggerClassFired(trigClass)) {
 281         AliDebug(2, Form("Ev[%4d] Trigger[%s]", fESDev->GetEventNumberInFile(), trigClass));
 282         kTRIGGERED = kTRUE;
 283         break;
 284       }
 285     }
 286     if(!kTRIGGERED){
 287       AliDebug(2, Form("Reject Ev[%4d] Trigger", fESDev->GetEventNumberInFile()));
 288       return;
 289     }
 290     // select only physical events
 291     if(fESDev->GetEventType() != 7){
 292       AliDebug(2, Form("Reject Ev[%4d] EvType[%d]", fESDev->GetEventNumberInFile(), fESDev->GetEventType()));
 293       return;
 294     }
 295   }
 296
 297   // if the required trigger is a collision trigger then apply event vertex cut
 298   if(UseLocalEvSelection() && IsCollision()){
 299     const AliESDVertex *vertex = fESDev->GetPrimaryVertex();
 300     if(TMath::Abs(vertex->GetZv())<1.e-10 ||
 301        TMath::Abs(vertex->GetZv())>fgkEvVertexZ ||
 302        vertex->GetNContributors()<fgkEvVertexN) {
 303       AliDebug(2, Form("Reject Ev[%4d] Vertex Zv[%f] Nv[%d]", fESDev->GetEventNumberInFile(), TMath::Abs(vertex->GetZv()), vertex->GetNContributors()));
 304       return;
 305     }
 306   }
 307
 308   if(fEventCut && !fEventCut->IsSelected(fESDev, IsCollision())) return;
 309
 310   if(!fESDfriend){
 311     AliError("Failed retrieving ESD friend event");
 312     return;
 313   }
 314   if(HasMCdata() && !fMCev){
 315     AliError("Failed retrieving MC event");
 316     return;
 317   }
 318
 319   fTracksBarrel->Delete();
 320   fTracksSA->Delete();
 321   fTracksKink->Delete();
 322   fV0List->Delete();
 323   fEventInfo->Delete("");
 324   new(fEventInfo)AliTRDeventInfo(fESDev->GetHeader(), const_cast<AliESDRun *>(fESDev->GetESDRun()));
 325
 326   Bool_t *trackMap(NULL);
 327   AliStack * mStack(NULL);
 328   Int_t nTracksMC = HasMCdata() ? fMCev->GetNumberOfTracks() : 0, nTracksESD = fESDev->GetNumberOfTracks();
 329   if(HasMCdata()){
 330     mStack = fMCev->Stack();
 331     if(!mStack){
 332       AliError("Failed retrieving MC Stack");
 333       return;
 334     }
 335     trackMap = new Bool_t[nTracksMC];
 336     memset(trackMap, 0, sizeof(Bool_t) * nTracksMC);
 337   }
 338
 339   Double32_t dedx[100]; Int_t nSlices(0);
 340   Int_t nTRDout(0), nTRDin(0), nTPC(0)
 341        ,nclsTrklt
 342        ,nBarrel(0), nSA(0), nKink(0)
 343        ,nBarrelMC(0), nSAMC(0), nKinkMC(0);
 344   AliESDtrack *esdTrack = NULL;
 345   AliESDfriendTrack *esdFriendTrack = NULL;
 346   TObject *calObject = NULL;
 347   AliTRDtrackV1 *track = NULL;
 348   AliTRDseedV1 *tracklet = NULL;
 349   AliTRDcluster *cl = NULL;
 350   // LOOP 0 - over ESD tracks
 351   for(Int_t itrk = 0; itrk < nTracksESD; itrk++){
 352     new(fTrackInfo) AliTRDtrackInfo();
 353     esdTrack = fESDev->GetTrack(itrk);
 354     AliDebug(3, Form("\n%3d ITS[%d] TPC[%d] TRD[%d]\n", itrk, esdTrack->GetNcls(0), esdTrack->GetNcls(1), esdTrack->GetNcls(2)));
 355
 356     if(esdTrack->GetStatus()&AliESDtrack::kTPCout) nTPC++;
 357     if(esdTrack->GetStatus()&AliESDtrack::kTRDout) nTRDout++;
 358     if(esdTrack->GetStatus()&AliESDtrack::kTRDin) nTRDin++;
 359
 360     // look at external track param
 361     const AliExternalTrackParam *op = esdTrack->GetOuterParam();
 362     Double_t xyz[3];
 363     if(op){
 364       op->GetXYZ(xyz);
 365       op->Global2LocalPosition(xyz, op->GetAlpha());
 366       AliDebug(3, Form("op @ X[%7.3f]\n", xyz[0]));
 367     }
 368
 369     // read MC info
 370     Int_t fPdg = -1;
 371     Int_t label = -1; UInt_t alab=UINT_MAX;
 372     if(HasMCdata()){
 373       label = esdTrack->GetLabel();
 374       alab = TMath::Abs(label);
 375       // register the track
 376       if(alab < UInt_t(nTracksMC)){
 377         trackMap[alab] = kTRUE;
 378       } else {
 379         AliError(Form("MC label[%d] outside scope for Ev[%d] Trk[%d].", label, (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), itrk));
 380         continue;
 381       }
 382       AliMCParticle *mcParticle = NULL;
 383       if(!(mcParticle = (AliMCParticle*) fMCev->GetTrack(alab))){
 384         AliError(Form("MC particle label[%d] missing for Ev[%d] Trk[%d].", label, (Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), itrk));
 385         continue;
 386       }
 387       fPdg = mcParticle->Particle()->GetPdgCode();
 388       Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
 389       Int_t iref = 0; AliTrackReference *ref = NULL;
 390       while(iref<nRefs){
 391         ref = mcParticle->GetTrackReference(iref);
 392         if(ref->LocalX() > fgkTPC) break;
 393         iref++;
 394       }
 395
 396       fTrackInfo->SetMC();
 397       fTrackInfo->SetPDG(fPdg);
 398       fTrackInfo->SetPrimary(mcParticle->Particle()->IsPrimary());
 399       fTrackInfo->SetLabel(label);
 400       Int_t jref = iref;//, kref = 0;
 401       while(jref<nRefs){
 402         ref = mcParticle->GetTrackReference(jref);
 403         if(ref->LocalX() > fgkTRD) break;
 404         AliDebug(4, Form("  trackRef[%2d (%2d)] @ %7.3f OK", jref-iref, jref, ref->LocalX()));
 405         fTrackInfo->AddTrackRef(ref);
 406         jref++;
 407       }
 408       AliDebug(3, Form("NtrackRefs[%d(%d)]", fTrackInfo->GetNTrackRefs(), nRefs));
 409     }
 410
 411     // copy some relevant info to TRD track info
 412     fTrackInfo->SetStatus(esdTrack->GetStatus());
 413     fTrackInfo->SetTrackId(esdTrack->GetID());
 414     Double_t p[AliPID::kSPECIES]; esdTrack->GetTRDpid(p);
 415     fTrackInfo->SetESDpid(p);
 416     fTrackInfo->SetESDpidQuality(esdTrack->GetTRDntrackletsPID());
 417     if(!nSlices) nSlices = esdTrack->GetNumberOfTRDslices();
 418     memset(dedx, 0, 100*sizeof(Double32_t));
 419     Int_t in(0);
 420     for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++)
 421       for(Int_t is=0; is<nSlices; is++)
 422         dedx[in++]=esdTrack->GetTRDslice(il, is);
 423     for(Int_t il=0; il<AliTRDgeometry::kNlayer; il++) dedx[in++]=esdTrack->GetTRDmomentum(il);
 424     fTrackInfo->SetSlices(in, dedx);
 425     fTrackInfo->SetNumberOfClustersRefit(esdTrack->GetNcls(2));
 426     // some other Informations which we may wish to store in order to find problematic cases
 427     fTrackInfo->SetKinkIndex(esdTrack->GetKinkIndex(0));
 428     fTrackInfo->SetTPCncls(static_cast<UShort_t>(esdTrack->GetNcls(1)));
 429     nclsTrklt = 0;
 430
 431
 432     // read REC info
 433     esdFriendTrack = fESDfriend->GetTrack(itrk);
 434     if(esdFriendTrack){
 435       Int_t icalib = 0;
 436       while((calObject = esdFriendTrack->GetCalibObject(icalib++))){
 437         if(strcmp(calObject->IsA()->GetName(),"AliTRDtrackV1") != 0) continue; // Look for the TRDtrack
 438         if(!(track = dynamic_cast<AliTRDtrackV1*>(calObject))) break;
 439         AliDebug(4, Form("TRD track OK"));
 440         // Set the clusters to unused
 441         for(Int_t ipl = 0; ipl < AliTRDgeometry::kNlayer; ipl++){
 442           if(!(tracklet = track->GetTracklet(ipl))) continue;
 443           tracklet->ResetClusterIter();
 444           while((cl = tracklet->NextCluster())) cl->Use(0);
 445         }
 446         fTrackInfo->SetTrack(track);
 447         break;
 448       }
 449       AliDebug(3, Form("Ntracklets[%d]\n", fTrackInfo->GetNTracklets()));
 450     } else AliDebug(3, "No ESD friends");
 451     if(op) fTrackInfo->SetOuterParam(op);
 452
 453     if(DebugLevel() >= 1){
 454       AliTRDtrackInfo info(*fTrackInfo);
 455       (*DebugStream()) << "trackInfo"
 456       << "TrackInfo.=" << &info
 457       << "\n";
 458       info.Delete("");
 459     }
 460
 461     ULong_t status(esdTrack->GetStatus());
 462     if((status&AliESDtrack::kTPCout)){
 463       if(!esdTrack->GetKinkIndex(0)){ // Barrel  Track Selection
 464         Bool_t selected(kTRUE);
 465         if(UseLocalTrkSelection()){
 466           if(esdTrack->Pt() < fgkPt){
 467             AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] Pt[%5.2f]", fESDev->GetEventNumberInFile(), itrk, esdTrack->Pt()));
 468             selected = kFALSE;
 469           }
 470           if(selected && TMath::Abs(esdTrack->Eta()) > fgkEta){
 471             AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] Eta[%5.2f]", fESDev->GetEventNumberInFile(), itrk, TMath::Abs(esdTrack->Eta())));
 472             selected = kFALSE;
 473           }
 474           if(selected && esdTrack->GetTPCNcls() < fgkNclTPC){
 475             AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] NclTPC[%d]", fESDev->GetEventNumberInFile(), itrk, esdTrack->GetTPCNcls()));
 476             selected = kFALSE;
 477           }
 478           Float_t par[2], cov[3];
 479           esdTrack->GetImpactParameters(par, cov);
 480           if(IsCollision()){ // cuts on DCA
 481             if(selected && TMath::Abs(par[0]) > fgkTrkDCAxy){
 482               AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] DCAxy[%f]", fESDev->GetEventNumberInFile(), itrk, TMath::Abs(par[0])));
 483               selected = kFALSE;
 484             }
 485             if(selected && TMath::Abs(par[1]) > fgkTrkDCAz){
 486               AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] DCAz[%f]", fESDev->GetEventNumberInFile(), itrk, TMath::Abs(par[1])));
 487               selected = kFALSE;
 488             }
 489           } else if(selected && fMCev && !fMCev->IsPhysicalPrimary(alab)){;
 490             AliDebug(2, Form("Reject Ev[%4d] Trk[%3d] Primary", fESDev->GetEventNumberInFile(), itrk));
 491             selected = kFALSE;
 492           }
 493         }
 494         if(fTrackCut && !fTrackCut->IsSelected(esdTrack)) selected = kFALSE;
 495         if(selected){
 496           fTracksBarrel->Add(new AliTRDtrackInfo(*fTrackInfo));
 497           nBarrel++;
 498         }
 499       } else {
 500         fTracksKink->Add(new AliTRDtrackInfo(*fTrackInfo));
 501         nKink++;
 502       }
 503     } else if((status&AliESDtrack::kTRDout) && !(status&AliESDtrack::kTRDin)){
 504       fTracksSA->Add(new AliTRDtrackInfo(*fTrackInfo));
 505       nSA++;
 506     }
 507     fTrackInfo->Delete("");
 508   }
 509
 510   // LOOP 1 - over ESD v0s
 511   Bool_t kFOUND(kFALSE);
 512   Float_t bField(fESDev->GetMagneticField());
 513   AliESDv0 *v0(NULL);
 514   for(Int_t iv0(0); iv0<fESDev->GetNumberOfV0s(); iv0++){
 515     if(!(v0 = fESDev->GetV0(iv0))) continue;
 516
 517     // purge V0 irelevant for TRD
 518     Int_t jb(0); AliTRDtrackInfo *ti(NULL);
 519     for(Int_t ib(0); ib<nBarrel; ib++){
 520       ti =   (AliTRDtrackInfo*)fTracksBarrel->At(ib);
 521       Int_t id(ti->GetTrackId());
 522       if(id!=v0->GetPindex() && id!=v0->GetNindex()) continue;
 523       kFOUND=kTRUE; ti->SetV0(); jb=ib+1;
 524       break;
 525     }
 526     if(!kFOUND) continue;
 527
 528     // register v0
 529     if(fV0Cut) new(fV0Info) AliTRDv0Info(*fV0Cut);
 530     else new(fV0Info) AliTRDv0Info();
 531     fV0Info->SetMagField(bField);
 532     fV0Info->SetV0tracks(fESDev->GetTrack(v0->GetPindex()), fESDev->GetTrack(v0->GetNindex()));
 533     fV0Info->SetV0Info(v0);
 534     // mark the other leg too
 535     for(Int_t ib(jb); ib<nBarrel; ib++){
 536       ti =   (AliTRDtrackInfo*)fTracksBarrel->At(ib);
 537       if(!fV0Info->HasTrack(ti)) continue;
 538       ti->SetV0();
 539       break;
 540     }
 541     //fV0Info->Print("a");
 542     fV0List->Add(new AliTRDv0Info(*fV0Info)); kFOUND=kFALSE;
 543   }
 544
 545   // LOOP 2 - over MC tracks which are passing TRD where the track is not reconstructed
 546   if(HasMCdata()){
 547     AliDebug(10, "Output of the MC track map:");
 548     for(Int_t itk = 0; itk < nTracksMC;  itk++) AliDebug(10, Form("trackMap[%d] = %s", itk, trackMap[itk] == kTRUE ? "TRUE" : "kFALSE"));
 549
 550     for(Int_t itk = 0; itk < nTracksMC; itk++){
 551       if(trackMap[itk]) continue;
 552       AliMCParticle *mcParticle =  (AliMCParticle*) fMCev->GetTrack(TMath::Abs(itk));
 553       Int_t fPdg = mcParticle->Particle()->GetPdgCode();
 554       Int_t nRefs = mcParticle->GetNumberOfTrackReferences();
 555       Int_t iref = 0; AliTrackReference *ref = NULL;
 556       Int_t nRefsTRD = 0;
 557       new(fTrackInfo) AliTRDtrackInfo();
 558       fTrackInfo->SetMC();
 559       fTrackInfo->SetPDG(fPdg);
 560       while(iref<nRefs){ // count TRD TR
 561         Bool_t kIN(kFALSE);
 562         ref = mcParticle->GetTrackReference(iref);
 563         if(ref->LocalX() > fgkTPC && ref->LocalX() < fgkTRD){
 564           fTrackInfo->AddTrackRef(ref);
 565           nRefsTRD++;kIN=kTRUE;
 566         }
 567         AliDebug(4, Form("  trackRef[%2d] @ x[%7.3f] %s", iref, ref->LocalX(), kIN?"IN":"OUT"));
 568         iref++;
 569       }
 570       if(!nRefsTRD){
 571         // In this stage we at least require 1 hit inside TRD. What will be done with this tracks is a task for the
 572         // analysis job
 573         fTrackInfo->Delete("");
 574         continue;
 575       }
 576       fTrackInfo->SetPrimary(mcParticle->Particle()->IsPrimary());
 577       fTrackInfo->SetLabel(itk);
 578       if(DebugLevel() >= 1){
 579         AliTRDtrackInfo info(*fTrackInfo);
 580         (*DebugStream()) << "trackInfo"
 581         << "TrackInfo.=" << &info
 582         << "\n";
 583         info.Delete("");
 584       }
 585       AliDebug(3, Form("Add MC track @ label[%d] nTRDrefs[%d].", itk, nRefsTRD));
 586       // check where the track starts
 587       ref = mcParticle->GetTrackReference(0);
 588       if(ref->LocalX() < fgkITS){
 589         fTracksBarrel->Add(new AliTRDtrackInfo(*fTrackInfo));
 590         nBarrelMC++;
 591       } else if(ref->LocalX() < fgkTPC) {
 592         fTracksKink->Add(new AliTRDtrackInfo(*fTrackInfo));
 593         nKinkMC++;
 594       } else if(nRefsTRD>6){
 595         fTracksSA->Add(new AliTRDtrackInfo(*fTrackInfo));
 596         nSAMC++;
 597       }
 598       fTrackInfo->Delete("");
 599     }
 600     delete[] trackMap;
 601   }
 602   AliDebug(1, Form(
 603     "\nEv[%3d] Tracks: ESD[%d] MC[%d] V0[%d]\n"
 604     "        TPCout[%d] TRDin[%d] TRDout[%d]\n"
 605     "        Barrel[%3d+%3d=%3d] SA[%2d+%2d=%2d] Kink[%2d+%2d=%2d]"
 606     ,(Int_t)AliAnalysisManager::GetAnalysisManager()->GetCurrentEntry(), nTracksESD, nTracksMC, fV0List->GetEntries()
 607     , nTPC, nTRDin, nTRDout
 608     ,nBarrel, nBarrelMC, fTracksBarrel->GetEntries()
 609     ,nSA, nSAMC, fTracksSA->GetEntries()
 610     ,nKink, nKinkMC, fTracksKink->GetEntries()
 611   ));
 612   // save track statistics
 613   TH1 *h((TH1S*)fContainer->At(0));
 614   h->Fill( 0., nTracksESD);
 615   h->Fill( 1., nTracksMC);
 616   h->Fill( 2., fV0List->GetEntries());
 617   h->Fill( 3., nTPC);
 618   h->Fill( 4., nTRDin);
 619   h->Fill( 5., nTRDout);
 620   h->Fill( 6., nBarrel);
 621   h->Fill( 7., nBarrelMC);
 622   h->Fill( 8., nSA);
 623   h->Fill( 9., nSAMC);
 624   h->Fill(10., nKink);
 625   h->Fill(11., nKinkMC);
 626
 627   PostData(kTracksBarrel, fTracksBarrel);
 628   PostData(kTracksSA, fTracksSA);
 629   PostData(kTracksKink, fTracksKink);
 630   PostData(kEventInfo, fEventInfo);
 631   PostData(kV0List, fV0List);
 632   PostData(kMonitor, fContainer);
 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