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