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