]> git.uio.no Git - u/mrichter/AliRoot.git/blob - ANALYSIS/AliAnalysisTaskESDfilter.cxx
Use of TOF of LG digits restored
[u/mrichter/AliRoot.git] / ANALYSIS / AliAnalysisTaskESDfilter.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$ */
17  
18 #include <TChain.h>
19 #include <TTree.h>
20 #include <TList.h>
21 #include <TArrayI.h>
22 #include <TParameter.h>
23 #include <TRandom.h>
24 #include <TParticle.h>
25 #include <TFile.h>
26
27 #include "AliAnalysisTaskESDfilter.h"
28 #include "AliAnalysisManager.h"
29 #include "AliESDEvent.h"
30 #include "AliESDRun.h"
31 #include "AliStack.h"
32 #include "AliAODEvent.h"
33 #include "AliMCEvent.h"
34 #include "AliMCEventHandler.h"
35 #include "AliESDInputHandler.h"
36 #include "AliAODHandler.h"
37 #include "AliAODMCParticle.h"
38 #include "AliAnalysisFilter.h"
39 #include "AliESDMuonTrack.h"
40 #include "AliESDVertex.h"
41 #include "AliCentrality.h"
42 #include "AliEventplane.h"
43 #include "AliESDv0.h"
44 #include "AliESDkink.h"
45 #include "AliESDcascade.h"
46 #include "AliESDPmdTrack.h"
47 #include "AliESDCaloCluster.h"
48 #include "AliESDCaloCells.h"
49 #include "AliMultiplicity.h"
50 #include "AliLog.h"
51 #include "AliCodeTimer.h"
52 #include "AliESDtrackCuts.h"
53 #include "AliESDpid.h"
54 #include "AliAODHMPIDrings.h"
55 #include "AliV0vertexer.h"
56 #include "AliCascadeVertexer.h"
57 #include "Riostream.h"
58 #include "AliExternalTrackParam.h"
59 #include "AliTrackerBase.h"
60 #include "TVector3.h"
61 #include "AliTPCdEdxInfo.h"
62
63 using std::cout;
64 using std::endl;
65 ClassImp(AliAnalysisTaskESDfilter)
66
67 ////////////////////////////////////////////////////////////////////////
68
69 AliAnalysisTaskESDfilter::AliAnalysisTaskESDfilter():
70   AliAnalysisTaskSE(),
71   fTrackFilter(0x0),
72   fKinkFilter(0x0),
73   fV0Filter(0x0),
74   fCascadeFilter(0x0),
75   fHighPthreshold(0),
76   fPtshape(0x0),
77   fEnableFillAOD(kTRUE),
78   fUsedTrack(0x0),
79   fUsedKink(0x0),
80   fUsedV0(0x0),
81   fAODTrackRefs(0x0),
82   fAODV0VtxRefs(0x0),
83   fAODV0Refs(0x0),
84   fMChandler(0x0),
85   fNumberOfTracks(0),
86   fNumberOfPositiveTracks(0),
87   fNumberOfV0s(0),
88   fNumberOfVertices(0),
89   fNumberOfCascades(0),
90   fNumberOfKinks(0),
91   fOldESDformat(kFALSE),
92   fPrimaryVertex(0x0),
93   fTPCConstrainedFilterMask(0),
94   fHybridFilterMaskTPCCG(0),
95   fWriteHybridTPCCOnly(kFALSE),
96   fGlobalConstrainedFilterMask(0),
97   fHybridFilterMaskGCG(0),
98   fWriteHybridGCOnly(kFALSE),
99   fIsVZEROEnabled(kTRUE),
100   fIsTZEROEnabled(kTRUE),
101   fIsZDCEnabled(kTRUE),
102   fIsHMPIDEnabled(kTRUE), 
103   fIsV0CascadeRecoEnabled(kFALSE),
104   fAreCascadesEnabled(kTRUE),
105   fAreV0sEnabled(kTRUE),
106   fAreKinksEnabled(kTRUE),
107   fAreTracksEnabled(kTRUE),
108   fArePmdClustersEnabled(kTRUE),
109   fAreCaloClustersEnabled(kTRUE),
110   fAreEMCALCellsEnabled(kTRUE),
111   fArePHOSCellsEnabled(kTRUE),
112   fAreEMCALTriggerEnabled(kTRUE),
113   fArePHOSTriggerEnabled(kTRUE),
114   fAreTrackletsEnabled(kTRUE),
115   fESDpid(0x0),
116   fIsPidOwner(kFALSE),
117   fTPCaloneTrackCuts(0),
118   fDoPropagateTrackToEMCal(kTRUE),
119   fEMCalSurfaceDistance(440)
120 {
121   // Default constructor
122     fV0Cuts[0] =  33.   ;   // max allowed chi2
123     fV0Cuts[1] =   0.1  ;   // min allowed impact parameter for the 1st daughter
124     fV0Cuts[2] =   0.1  ;   // min allowed impact parameter for the 2nd daughter
125     fV0Cuts[3] =   1.   ;   // max allowed DCA between the daughter tracks
126     fV0Cuts[4] =    .998;   // min allowed cosine of V0's pointing angle
127     fV0Cuts[5] =   0.9  ;   // min radius of the fiducial volume
128     fV0Cuts[6] = 100.   ;   // max radius of the fiducial volume
129
130     fCascadeCuts[0] =  33.   ; // max allowed chi2 (same as PDC07)
131     fCascadeCuts[1] =   0.05 ; // min allowed V0 impact parameter
132     fCascadeCuts[2] =   0.008; // "window" around the Lambda mass
133     fCascadeCuts[3] =   0.03 ; // min allowed bachelor's impact parameter
134     fCascadeCuts[4] =   0.3  ; // max allowed DCA between the V0 and the bachelor
135     fCascadeCuts[5] =   0.999; // min allowed cosine of the cascade pointing angle
136     fCascadeCuts[6] =   0.9  ; // min radius of the fiducial volume
137     fCascadeCuts[7] = 100.   ; // max radius of the fiducial volume
138 }
139
140 //______________________________________________________________________________
141 AliAnalysisTaskESDfilter::AliAnalysisTaskESDfilter(const char* name):
142     AliAnalysisTaskSE(name),
143     fTrackFilter(0x0),
144     fKinkFilter(0x0),
145     fV0Filter(0x0),
146     fCascadeFilter(0x0),
147     fHighPthreshold(0),
148     fPtshape(0x0),
149     fEnableFillAOD(kTRUE),
150     fUsedTrack(0x0),
151     fUsedKink(0x0),
152     fUsedV0(0x0),
153     fAODTrackRefs(0x0),
154     fAODV0VtxRefs(0x0),
155     fAODV0Refs(0x0),
156     fMChandler(0x0),
157     fNumberOfTracks(0),
158     fNumberOfPositiveTracks(0),
159     fNumberOfV0s(0),
160     fNumberOfVertices(0),
161     fNumberOfCascades(0),
162     fNumberOfKinks(0),
163     fOldESDformat(kFALSE),
164     fPrimaryVertex(0x0),
165   fTPCConstrainedFilterMask(0),
166   fHybridFilterMaskTPCCG(0),
167   fWriteHybridTPCCOnly(kFALSE),
168   fGlobalConstrainedFilterMask(0),
169   fHybridFilterMaskGCG(0),
170   fWriteHybridGCOnly(kFALSE),
171     fIsVZEROEnabled(kTRUE),
172     fIsTZEROEnabled(kTRUE),
173     fIsZDCEnabled(kTRUE),
174     fIsHMPIDEnabled(kTRUE), 
175     fIsV0CascadeRecoEnabled(kFALSE),
176     fAreCascadesEnabled(kTRUE),
177     fAreV0sEnabled(kTRUE),
178     fAreKinksEnabled(kTRUE),
179     fAreTracksEnabled(kTRUE),
180     fArePmdClustersEnabled(kTRUE),
181     fAreCaloClustersEnabled(kTRUE),
182     fAreEMCALCellsEnabled(kTRUE),
183     fArePHOSCellsEnabled(kTRUE),
184                 fAreEMCALTriggerEnabled(kTRUE),
185                 fArePHOSTriggerEnabled(kTRUE),
186                 fAreTrackletsEnabled(kTRUE),
187     fESDpid(0x0),
188     fIsPidOwner(kFALSE),
189     fTPCaloneTrackCuts(0),
190     fDoPropagateTrackToEMCal(kTRUE),
191     fEMCalSurfaceDistance(440)
192 {
193   // Constructor
194
195     fV0Cuts[0] =  33.   ;   // max allowed chi2
196     fV0Cuts[1] =   0.1  ;   // min allowed impact parameter for the 1st daughter
197     fV0Cuts[2] =   0.1  ;   // min allowed impact parameter for the 2nd daughter
198     fV0Cuts[3] =   1.   ;   // max allowed DCA between the daughter tracks
199     fV0Cuts[4] =    .998;   // min allowed cosine of V0's pointing angle
200     fV0Cuts[5] =   0.9  ;   // min radius of the fiducial volume
201     fV0Cuts[6] = 100.   ;   // max radius of the fiducial volume
202
203     fCascadeCuts[0] =  33.   ; // max allowed chi2 (same as PDC07)
204     fCascadeCuts[1] =   0.05 ; // min allowed V0 impact parameter
205     fCascadeCuts[2] =   0.008; // "window" around the Lambda mass
206     fCascadeCuts[3] =   0.03 ; // min allowed bachelor's impact parameter
207     fCascadeCuts[4] =   0.3  ; // max allowed DCA between the V0 and the bachelor
208     fCascadeCuts[5] =   0.999; // min allowed cosine of the cascade pointing angle
209     fCascadeCuts[6] =   0.9  ; // min radius of the fiducial volume
210     fCascadeCuts[7] = 100.   ; // max radius of the fiducial volume
211
212
213
214 }
215 AliAnalysisTaskESDfilter::~AliAnalysisTaskESDfilter(){
216     if(fIsPidOwner) delete fESDpid;
217 }
218 //______________________________________________________________________________
219 void AliAnalysisTaskESDfilter::UserCreateOutputObjects()
220 {
221   //
222   // Create Output Objects conenct filter to outputtree
223   // 
224   if(OutputTree())
225   {
226     OutputTree()->GetUserInfo()->Add(fTrackFilter);
227   }
228   else
229   {
230     AliError("No OutputTree() for adding the track filter");
231   }
232   fTPCaloneTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
233 }
234
235 //______________________________________________________________________________
236 void AliAnalysisTaskESDfilter::Init()
237 {
238   // Initialization
239   if (fDebug > 1) AliInfo("Init() \n");
240   // Call configuration file
241 }
242
243 //______________________________________________________________________________
244 Bool_t AliAnalysisTaskESDfilter::Notify()
245 {
246 // Notify method.
247    AddMetadataToUserInfo();
248    return kTRUE;
249 }   
250
251 //______________________________________________________________________________
252 Bool_t AliAnalysisTaskESDfilter::AddMetadataToUserInfo()
253 {
254 // Copy metadata to AOD user info.
255   static Bool_t copyFirst = kFALSE;
256   if (!copyFirst) {
257     AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
258     if (!mgr) {
259        AliError("AliAnalysisTaskESDfilter::AddMetadataToUserInfo() : No analysis manager !");
260        return kFALSE;
261     }   
262     TTree *esdTree = mgr->GetTree()->GetTree();
263     if (!esdTree) return kFALSE;
264     TNamed *alirootVersion = (TNamed*)esdTree->GetUserInfo()->FindObject("alirootVersion");
265     if (!alirootVersion) return kFALSE;    
266     AliAODHandler *aodHandler = dynamic_cast<AliAODHandler*>(mgr->GetOutputEventHandler());
267          if (!aodHandler) return kFALSE;
268     TTree *aodTree = aodHandler->GetTree();
269     if (!aodTree) return kFALSE;
270     aodTree->GetUserInfo()->Add(new TNamed(*alirootVersion));
271     copyFirst = kTRUE;
272   }
273   return kTRUE;
274 }
275
276 //______________________________________________________________________________
277 void AliAnalysisTaskESDfilter::PrintTask(Option_t *option, Int_t indent) const
278 {
279 // Print selection task information
280   AliInfo("");
281   
282   AliAnalysisTaskSE::PrintTask(option,indent);
283   
284   TString spaces(' ',indent+3);
285   
286         cout << spaces.Data() << Form("Cascades       are %s",fAreCascadesEnabled ? "ENABLED":"DISABLED") << endl;
287         cout << spaces.Data() << Form("V0s            are %s",fAreV0sEnabled ? "ENABLED":"DISABLED") << endl;
288         cout << spaces.Data() << Form("Kinks          are %s",fAreKinksEnabled ? "ENABLED":"DISABLED") << endl;
289         cout << spaces.Data() << Form("Tracks         are %s",fAreTracksEnabled ? "ENABLED":"DISABLED") << endl;
290         cout << spaces.Data() << Form("PmdClusters    are %s",fArePmdClustersEnabled ? "ENABLED":"DISABLED") << endl;
291         cout << spaces.Data() << Form("CaloClusters   are %s",fAreCaloClustersEnabled ? "ENABLED":"DISABLED") << endl;
292   cout << spaces.Data() << Form("EMCAL cells    are %s",fAreEMCALCellsEnabled ? "ENABLED":"DISABLED") << endl;
293         cout << spaces.Data() << Form("EMCAL triggers are %s",fAreEMCALTriggerEnabled ? "ENABLED":"DISABLED") << endl;
294         cout << spaces.Data() << Form("PHOS triggers  are %s",fArePHOSTriggerEnabled ? "ENABLED":"DISABLED") << endl;
295         cout << spaces.Data() << Form("Tracklets      are %s",fAreTrackletsEnabled ? "ENABLED":"DISABLED") << endl;  
296         cout << spaces.Data() << Form("PropagateTrackToEMCal  is %s", fDoPropagateTrackToEMCal ? "ENABLED":"DISABLED") << endl; 
297 }
298
299 //______________________________________________________________________________
300 void AliAnalysisTaskESDfilter::UserExec(Option_t */*option*/)
301 {
302 // Execute analysis for current event
303 //
304                                             
305   Long64_t ientry = Entry();
306   
307   if (fDebug > 0) {
308       printf("Filter: Analysing event # %5d\n", (Int_t) ientry);
309       if (fHighPthreshold == 0) AliInfo("detector PID signals are stored in each track");
310       if (!fPtshape) AliInfo("detector PID signals are not stored below the pt threshold");
311   }
312   // Filters must explicitely enable AOD filling in their UserExec (AG)
313   if (!AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()) AliFatal("Cannot run ESD filter without an output event handler");
314   if(fEnableFillAOD) {
315      AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
316      AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillExtension(kTRUE);
317   }   
318   ConvertESDtoAOD();
319 }
320
321 //______________________________________________________________________________
322 TClonesArray& AliAnalysisTaskESDfilter::Cascades()
323 {
324   return *(AODEvent()->GetCascades());
325 }
326
327 //______________________________________________________________________________
328 TClonesArray& AliAnalysisTaskESDfilter::Tracks()
329 {
330   return *(AODEvent()->GetTracks());
331 }
332
333 //______________________________________________________________________________
334 TClonesArray& AliAnalysisTaskESDfilter::V0s()
335 {
336   return *(AODEvent()->GetV0s());
337 }
338
339 //______________________________________________________________________________
340 TClonesArray& AliAnalysisTaskESDfilter::Vertices()
341 {
342   return *(AODEvent()->GetVertices());
343 }
344
345 //______________________________________________________________________________
346 AliAODHeader* AliAnalysisTaskESDfilter::ConvertHeader(const AliESDEvent& esd)
347 {
348 // Convert header information
349
350   AliCodeTimerAuto("",0);
351   
352   AliAODHeader* header = AODEvent()->GetHeader();
353   
354   header->SetRunNumber(esd.GetRunNumber());
355   header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
356
357   TTree* tree = fInputHandler->GetTree();
358   if (tree) {
359     TFile* file = tree->GetCurrentFile();
360     if (file) header->SetESDFileName(file->GetName());
361   }
362   
363   if (fOldESDformat) {
364     header->SetBunchCrossNumber(0);
365     header->SetOrbitNumber(0);
366     header->SetPeriodNumber(0);
367     header->SetEventType(0);
368     header->SetMuonMagFieldScale(-999.);
369     header->SetCentrality(0);       
370     header->SetEventplane(0);
371   } else {
372     header->SetBunchCrossNumber(esd.GetBunchCrossNumber());
373     header->SetOrbitNumber(esd.GetOrbitNumber());
374     header->SetPeriodNumber(esd.GetPeriodNumber());
375     header->SetEventType(esd.GetEventType());
376     
377     header->SetEventNumberESDFile(esd.GetHeader()->GetEventNumberInFile());
378     if(const_cast<AliESDEvent&>(esd).GetCentrality()){
379       header->SetCentrality(const_cast<AliESDEvent&>(esd).GetCentrality());
380     }
381     else{
382       header->SetCentrality(0);
383     }
384     if(const_cast<AliESDEvent&>(esd).GetEventplane()){
385       header->SetEventplane(const_cast<AliESDEvent&>(esd).GetEventplane());
386     }
387     else{
388       header->SetEventplane(0);
389     }
390   }
391   
392   // Trigger
393   header->SetFiredTriggerClasses(esd.GetFiredTriggerClasses());
394   header->SetTriggerMask(esd.GetTriggerMask()); 
395   header->SetTriggerCluster(esd.GetTriggerCluster());
396   header->SetL0TriggerInputs(esd.GetHeader()->GetL0TriggerInputs());    
397   header->SetL1TriggerInputs(esd.GetHeader()->GetL1TriggerInputs());    
398   header->SetL2TriggerInputs(esd.GetHeader()->GetL2TriggerInputs());    
399   
400   header->SetMagneticField(esd.GetMagneticField());
401   header->SetMuonMagFieldScale(esd.GetCurrentDip()/6000.);
402   header->SetZDCN1Energy(esd.GetZDCN1Energy());
403   header->SetZDCP1Energy(esd.GetZDCP1Energy());
404   header->SetZDCN2Energy(esd.GetZDCN2Energy());
405   header->SetZDCP2Energy(esd.GetZDCP2Energy());
406   header->SetZDCEMEnergy(esd.GetZDCEMEnergy(0),esd.GetZDCEMEnergy(1));
407
408   header->SetIRInt2InteractionMap(esd.GetHeader()->GetIRInt2InteractionMap());
409   header->SetIRInt1InteractionMap(esd.GetHeader()->GetIRInt1InteractionMap());
410   
411   // ITS Cluster Multiplicty
412   const AliMultiplicity *mult = esd.GetMultiplicity();
413   for (Int_t ilay = 0; ilay < 6; ilay++) header->SetITSClusters(ilay, mult->GetNumberOfITSClusters(ilay));
414   
415   // TPC only Reference Multiplicty
416   Int_t refMult  = fTPCaloneTrackCuts ? (Short_t)fTPCaloneTrackCuts->GetReferenceMultiplicity(&esd, kTRUE) : -1;
417   header->SetTPConlyRefMultiplicity(refMult);
418   //
419   AliESDtrackCuts::MultEstTrackType estType = esd.GetPrimaryVertexTracks()->GetStatus() ? AliESDtrackCuts::kTrackletsITSTPC : AliESDtrackCuts::kTracklets;
420   header->SetRefMultiplicityComb05(AliESDtrackCuts::GetReferenceMultiplicity(&esd,estType,0.5));
421   header->SetRefMultiplicityComb08(AliESDtrackCuts::GetReferenceMultiplicity(&esd,estType,0.8));
422   //
423   Float_t diamxy[2]={esd.GetDiamondX(),esd.GetDiamondY()};
424   Float_t diamcov[3]; 
425   esd.GetDiamondCovXY(diamcov);
426   header->SetDiamond(diamxy,diamcov);
427   header->SetDiamondZ(esd.GetDiamondZ(),esd.GetSigma2DiamondZ());
428   
429   // VZERO channel equalization factors for event-plane reconstruction   
430   header->SetVZEROEqFactors(esd.GetVZEROEqFactors());
431
432   // T0 Resolution information                                                                                                                                          
433   const AliESDRun* esdRun = esd.GetESDRun();
434   for (Int_t i=0;i<AliESDRun::kT0spreadSize;i++) header->SetT0spread(i,esdRun->GetT0spread(i));
435
436   return header;
437 }
438
439 //______________________________________________________________________________
440 void AliAnalysisTaskESDfilter::ConvertCascades(const AliESDEvent& esd) 
441 {
442
443   // Convert the cascades part of the ESD.
444   // Return the number of cascades
445  
446   AliCodeTimerAuto("",0);
447   
448   // Create vertices starting from the most complex objects
449   Double_t chi2 = 0.;
450   
451   const AliESDVertex* vtx = esd.GetPrimaryVertex();
452   Double_t pos[3] = { 0. };
453   Double_t covVtx[6] = { 0. };
454   Double_t momBach[3]={0.};
455   Double_t covTr[21]={0.};
456   Double_t pid[10]={0.};
457   AliAODPid* detpid(0x0);
458   AliAODVertex* vV0FromCascade(0x0);
459   AliAODv0* aodV0(0x0);
460   AliAODcascade* aodCascade(0x0);
461   AliAODTrack* aodTrack(0x0);
462   Double_t momPos[3]={0.};
463   Double_t momNeg[3] = { 0. };
464   Double_t momPosAtV0vtx[3]={0.};
465   Double_t momNegAtV0vtx[3]={0.};
466
467   TClonesArray& verticesArray = Vertices();
468   TClonesArray& tracksArray = Tracks();
469   TClonesArray& cascadesArray = Cascades();
470   
471   // Cascades (Modified by A.Maire - February 2009)
472   for (Int_t nCascade = 0; nCascade < esd.GetNumberOfCascades(); ++nCascade) {
473     
474     // 0- Preparation
475     //
476     AliESDcascade *esdCascade = esd.GetCascade(nCascade);
477                 Int_t  idxPosFromV0Dghter  = esdCascade->GetPindex();
478                 Int_t  idxNegFromV0Dghter  = esdCascade->GetNindex();
479                 Int_t  idxBachFromCascade  = esdCascade->GetBindex();
480     
481     AliESDtrack  *esdCascadePos  = esd.GetTrack( idxPosFromV0Dghter);
482     AliESDtrack  *esdCascadeNeg  = esd.GetTrack( idxNegFromV0Dghter);
483     AliESDtrack  *esdCascadeBach = esd.GetTrack( idxBachFromCascade);
484     
485     // Identification of the V0 within the esdCascade (via both daughter track indices)
486     AliESDv0 * currentV0   = 0x0;
487     Int_t      idxV0FromCascade = -1;
488     
489     for (Int_t iV0=0; iV0<esd.GetNumberOfV0s(); ++iV0) {
490       
491       currentV0 = esd.GetV0(iV0);
492       Int_t posCurrentV0 = currentV0->GetPindex();
493       Int_t negCurrentV0 = currentV0->GetNindex();
494       
495       if (posCurrentV0==idxPosFromV0Dghter && negCurrentV0==idxNegFromV0Dghter) {
496         idxV0FromCascade = iV0;
497         break;
498       }
499     }
500     
501     if(idxV0FromCascade < 0){
502       printf("Cascade - no matching for the V0 (index V0 = -1) ! Skip ... \n");
503       continue;
504     }// a priori, useless check, but safer ... in case of pb with tracks "out of bounds"
505     
506     AliESDv0 *esdV0FromCascade   = esd.GetV0(idxV0FromCascade);
507         
508     // 1 - Cascade selection 
509     
510     //  AliESDVertex *esdPrimVtx = new AliESDVertex(*(esd.GetPrimaryVertex()));
511     //  TList cascadeObjects;
512     //  cascadeObjects.AddAt(esdV0FromCascade, 0);
513     //  cascadeObjects.AddAt(esdCascadePos,    1);
514     //  cascadeObjects.AddAt(esdCascadeNeg,    2);
515     //  cascadeObjects.AddAt(esdCascade,       3);
516     //  cascadeObjects.AddAt(esdCascadeBach,   4);
517     //  cascadeObjects.AddAt(esdPrimVtx,       5);
518     // 
519     //  UInt_t selectCascade = 0;
520     //  if (fCascadeFilter) {
521     //    // selectCascade = fCascadeFilter->IsSelected(&cascadeObjects); 
522     //          // FIXME AliESDCascadeCuts to be implemented ...
523     // 
524     //          // Here we may encounter a moot point at the V0 level 
525     //          // between the cascade selections and the V0 ones :
526     //          // the V0 selected along with the cascade (secondary V0) may 
527     //          // usually be removed from the dedicated V0 selections (prim V0) ...
528     //          // -> To be discussed !
529     // 
530     //    // this is a little awkward but otherwise the 
531     //    // list wants to access the pointer (delete it) 
532     //    // again when going out of scope
533     //    delete cascadeObjects.RemoveAt(5); // esdPrimVtx created via copy construct
534     //    esdPrimVtx = 0;
535     //    if (!selectCascade) 
536     //      continue;
537     //  }
538     //  else{
539     //    delete cascadeObjects.RemoveAt(5); // esdPrimVtx created via copy construct
540     //    esdPrimVtx = 0;
541     //  }
542     
543     // 2 - Add the cascade vertex
544     
545     esdCascade->GetXYZcascade(pos[0], pos[1], pos[2]);
546     esdCascade->GetPosCovXi(covVtx);
547     chi2 = esdCascade->GetChi2Xi(); 
548     
549     AliAODVertex *vCascade = new(verticesArray[fNumberOfVertices++]) AliAODVertex( pos,
550                                                                      covVtx,
551                                                                      chi2, // FIXME = Chi2/NDF will be needed
552                                                                      fPrimaryVertex,
553                                                                      nCascade, // id
554                                                                      AliAODVertex::kCascade);
555     fPrimaryVertex->AddDaughter(vCascade);
556     
557 //    if (fDebug > 2) {
558 //      printf("---- Cascade / Cascade Vertex (AOD) : \n");
559 //      vCascade->Print();
560 //    }
561     
562 //    if(esd.GetTOFHeader() && fIsPidOwner) fESDpid->SetTOFResponse(const_cast<AliESDEvent*>(&esd), (AliESDpid::EStartTimeType_t)fTimeZeroType); //in case of AOD production starting form LHC10e without Tender. 
563
564
565     // 3 - Add the bachelor track from the cascade
566     
567     if (!fUsedTrack[idxBachFromCascade]) {
568       
569       esdCascadeBach->GetPxPyPz(momBach);
570       esdCascadeBach->GetXYZ(pos);
571       esdCascadeBach->GetCovarianceXYZPxPyPz(covTr);
572       esdCascadeBach->GetESDpid(pid);
573       
574             fUsedTrack[idxBachFromCascade] = kTRUE;
575             UInt_t selectInfo = 0;
576             if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdCascadeBach);
577             if (fMChandler) fMChandler->SelectParticle(esdCascadeBach->GetLabel());
578             aodTrack = new(tracksArray[fNumberOfTracks++]) AliAODTrack(esdCascadeBach->GetID(),
579                                                                            esdCascadeBach->GetLabel(), 
580                                                                            momBach, 
581                                                                            kTRUE,
582                                                                            pos,
583                                                                            kFALSE, // Why kFALSE for "isDCA" ? FIXME
584                                                                            covTr, 
585                                                                            (Short_t)esdCascadeBach->GetSign(),
586                                                                            esdCascadeBach->GetITSClusterMap(), 
587                                                                            pid,
588                                                                            vCascade,
589                                                                            kTRUE,  // usedForVtxFit = kFALSE ? FIXME
590                                                                            vtx->UsesTrack(esdCascadeBach->GetID()),
591                                                                            AliAODTrack::kSecondary,
592                                                                            selectInfo);
593             aodTrack->SetTPCFitMap(esdCascadeBach->GetTPCFitMap());
594             aodTrack->SetTPCClusterMap(esdCascadeBach->GetTPCClusterMap());
595             aodTrack->SetTPCSharedMap (esdCascadeBach->GetTPCSharedMap());
596             aodTrack->SetChi2perNDF(Chi2perNDF(esdCascadeBach));
597             aodTrack->SetTPCPointsF(esdCascadeBach->GetTPCNclsF());
598             aodTrack->SetTPCNCrossedRows(UShort_t(esdCascadeBach->GetTPCCrossedRows()));
599             aodTrack->SetIntegratedLength(esdCascadeBach->GetIntegratedLength());
600             fAODTrackRefs->AddAt(aodTrack,idxBachFromCascade);
601             
602             if (esdCascadeBach->GetSign() > 0) ++fNumberOfPositiveTracks;
603             aodTrack->ConvertAliPIDtoAODPID();
604             aodTrack->SetFlags(esdCascadeBach->GetStatus());
605       SetAODPID(esdCascadeBach,aodTrack,detpid);
606     }
607     else {
608             aodTrack = static_cast<AliAODTrack*>( fAODTrackRefs->At(idxBachFromCascade) );
609     }
610     
611     vCascade->AddDaughter(aodTrack);
612     
613 //    if (fDebug > 4) {
614 //      printf("---- Cascade / bach dghter : \n");
615 //      aodTrack->Print();
616 //    }
617     
618     
619     // 4 - Add the V0 from the cascade. 
620     // = V0vtx + both pos and neg daughter tracks + the aodV0 itself
621     //
622     
623     if ( !fUsedV0[idxV0FromCascade] ) {
624       // 4.A - if VO structure hasn't been created yet
625       
626       // 4.A.1 - Create the V0 vertex of the cascade
627       
628       esdV0FromCascade->GetXYZ(pos[0], pos[1], pos[2]);
629       esdV0FromCascade->GetPosCov(covVtx);
630       chi2 = esdV0FromCascade->GetChi2V0();  // = chi2/NDF since NDF = 2*2-3 ?
631                         
632       vV0FromCascade = new(verticesArray[fNumberOfVertices++]) AliAODVertex(pos,
633                                                                          covVtx,
634                                                                          chi2,
635                                                                          vCascade,
636                                                                          idxV0FromCascade, //id of ESDv0
637                                                                          AliAODVertex::kV0);
638       // Note:
639       //    one V0 can be used by several cascades.
640       // So, one AOD V0 vtx can have several parent vtx.
641       // This is not directly allowed by AliAODvertex.
642       // Setting the parent vtx (here = param "vCascade") doesn't lead to a crash
643       // but to a problem of consistency within AODEvent.
644       // -> See below paragraph 4.B, for the proposed treatment of such a case.
645       
646       // Add the vV0FromCascade to the aodVOVtxRefs
647       fAODV0VtxRefs->AddAt(vV0FromCascade,idxV0FromCascade);
648       
649       
650       // 4.A.2 - Add the positive tracks from the V0
651       
652       esdCascadePos->GetPxPyPz(momPos);
653       esdCascadePos->GetXYZ(pos);
654       esdCascadePos->GetCovarianceXYZPxPyPz(covTr);
655       esdCascadePos->GetESDpid(pid);
656       
657       
658       if (!fUsedTrack[idxPosFromV0Dghter]) {
659         fUsedTrack[idxPosFromV0Dghter] = kTRUE;
660         
661         UInt_t selectInfo = 0;
662         if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdCascadePos);
663         if(fMChandler) fMChandler->SelectParticle(esdCascadePos->GetLabel());
664         aodTrack = new(tracksArray[fNumberOfTracks++]) 
665         AliAODTrack(  esdCascadePos->GetID(),
666                     esdCascadePos->GetLabel(), 
667                     momPos, 
668                     kTRUE,
669                     pos,
670                     kFALSE, // Why kFALSE for "isDCA" ? FIXME
671                     covTr, 
672                     (Short_t)esdCascadePos->GetSign(),
673                     esdCascadePos->GetITSClusterMap(), 
674                     pid,
675                     vV0FromCascade,
676                     kTRUE,  // usedForVtxFit = kFALSE ? FIXME
677                     vtx->UsesTrack(esdCascadePos->GetID()),
678                     AliAODTrack::kSecondary,
679                     selectInfo);
680         aodTrack->SetTPCFitMap(esdCascadePos->GetTPCFitMap());
681         aodTrack->SetTPCClusterMap(esdCascadePos->GetTPCClusterMap());
682         aodTrack->SetTPCSharedMap (esdCascadePos->GetTPCSharedMap());
683         aodTrack->SetChi2perNDF(Chi2perNDF(esdCascadePos));
684         aodTrack->SetTPCPointsF(esdCascadePos->GetTPCNclsF());
685         aodTrack->SetTPCNCrossedRows(UShort_t(esdCascadePos->GetTPCCrossedRows()));
686         aodTrack->SetIntegratedLength(esdCascadePos->GetIntegratedLength());
687         fAODTrackRefs->AddAt(aodTrack,idxPosFromV0Dghter);
688         
689         if (esdCascadePos->GetSign() > 0) ++fNumberOfPositiveTracks;
690         aodTrack->ConvertAliPIDtoAODPID();
691         aodTrack->SetFlags(esdCascadePos->GetStatus());
692         SetAODPID(esdCascadePos,aodTrack,detpid);
693       }
694       else {
695         aodTrack = static_cast<AliAODTrack*>(fAODTrackRefs->At(idxPosFromV0Dghter));
696       }
697       vV0FromCascade->AddDaughter(aodTrack);
698       
699       
700       // 4.A.3 - Add the negative tracks from the V0
701       
702       esdCascadeNeg->GetPxPyPz(momNeg);
703       esdCascadeNeg->GetXYZ(pos);
704       esdCascadeNeg->GetCovarianceXYZPxPyPz(covTr);
705       esdCascadeNeg->GetESDpid(pid);
706       
707       
708       if (!fUsedTrack[idxNegFromV0Dghter]) {
709         fUsedTrack[idxNegFromV0Dghter] = kTRUE;
710         
711         UInt_t selectInfo = 0;
712         if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdCascadeNeg);
713         if(fMChandler)fMChandler->SelectParticle(esdCascadeNeg->GetLabel());
714         aodTrack = new(tracksArray[fNumberOfTracks++]) AliAODTrack(  esdCascadeNeg->GetID(),
715                                                       esdCascadeNeg->GetLabel(),
716                                                       momNeg,
717                                                       kTRUE,
718                                                       pos,
719                                                       kFALSE, // Why kFALSE for "isDCA" ? FIXME
720                                                       covTr, 
721                                                       (Short_t)esdCascadeNeg->GetSign(),
722                                                       esdCascadeNeg->GetITSClusterMap(), 
723                                                       pid,
724                                                       vV0FromCascade,
725                                                       kTRUE,  // usedForVtxFit = kFALSE ? FIXME
726                                                       vtx->UsesTrack(esdCascadeNeg->GetID()),
727                                                       AliAODTrack::kSecondary,
728                                                       selectInfo);
729         aodTrack->SetTPCFitMap(esdCascadeNeg->GetTPCFitMap());
730         aodTrack->SetTPCClusterMap(esdCascadeNeg->GetTPCClusterMap());
731         aodTrack->SetTPCSharedMap (esdCascadeNeg->GetTPCSharedMap());
732         aodTrack->SetChi2perNDF(Chi2perNDF(esdCascadeNeg));
733         aodTrack->SetTPCPointsF(esdCascadeNeg->GetTPCNclsF());
734         aodTrack->SetTPCNCrossedRows(UShort_t(esdCascadeNeg->GetTPCCrossedRows()));
735         aodTrack->SetIntegratedLength(esdCascadeNeg->GetIntegratedLength());
736         fAODTrackRefs->AddAt(aodTrack,idxNegFromV0Dghter);
737         
738         if (esdCascadeNeg->GetSign() > 0) ++fNumberOfPositiveTracks;
739         aodTrack->ConvertAliPIDtoAODPID();
740         aodTrack->SetFlags(esdCascadeNeg->GetStatus());
741         SetAODPID(esdCascadeNeg,aodTrack,detpid);
742       }
743       else {
744         aodTrack = static_cast<AliAODTrack*>(fAODTrackRefs->At(idxNegFromV0Dghter));
745       }
746       
747       vV0FromCascade->AddDaughter(aodTrack);
748       
749                         
750       // 4.A.4 - Add the V0 from cascade to the V0 array
751       
752       Double_t  dcaV0Daughters      = esdV0FromCascade->GetDcaV0Daughters();
753       Double_t  dcaV0ToPrimVertex   = esdV0FromCascade->GetD( esd.GetPrimaryVertex()->GetX(),
754                                                              esd.GetPrimaryVertex()->GetY(),
755                                                              esd.GetPrimaryVertex()->GetZ() );
756       esdV0FromCascade->GetPPxPyPz( momPosAtV0vtx[0],momPosAtV0vtx[1],momPosAtV0vtx[2] ); 
757       esdV0FromCascade->GetNPxPyPz( momNegAtV0vtx[0],momNegAtV0vtx[1],momNegAtV0vtx[2] ); 
758       
759       Double_t dcaDaughterToPrimVertex[2] = { 999., 999.}; // ..[0] = DCA in (x,y) for Pos and ..[1] = Neg
760       dcaDaughterToPrimVertex[0] = TMath::Abs(esdCascadePos->GetD(      esd.GetPrimaryVertex()->GetX(),
761                                                                   esd.GetPrimaryVertex()->GetY(),
762                                                                   esd.GetMagneticField())        );
763       dcaDaughterToPrimVertex[1] = TMath::Abs(esdCascadeNeg->GetD(      esd.GetPrimaryVertex()->GetX(),
764                                                                   esd.GetPrimaryVertex()->GetY(),
765                                                                   esd.GetMagneticField())        );
766       
767       aodV0 = new(V0s()[fNumberOfV0s++]) AliAODv0( vV0FromCascade, 
768                                                   dcaV0Daughters,
769                                                   dcaV0ToPrimVertex, 
770                                                   momPosAtV0vtx, 
771                                                   momNegAtV0vtx, 
772                                                   dcaDaughterToPrimVertex); 
773       // set the aod v0 on-the-fly status
774       aodV0->SetOnFlyStatus(esdV0FromCascade->GetOnFlyStatus());
775       
776       // Add the aodV0 to the aodVORefs
777       fAODV0Refs->AddAt(aodV0,idxV0FromCascade);
778       
779       fUsedV0[idxV0FromCascade] = kTRUE;
780       
781     } else { 
782       // 4.B - if V0 structure already used
783       
784       // Note :
785       //    one V0 can be used by several cascades (frequent in PbPb evts) : 
786       // same V0 which used but attached to different bachelor tracks
787       // -> aodVORefs and fAODV0VtxRefs are needed.
788       // Goal : avoid a redundancy of the info in "Vertices" and "v0s" clones array.
789       
790       vV0FromCascade = static_cast<AliAODVertex*>( fAODV0VtxRefs->At(idxV0FromCascade) );
791       aodV0          = static_cast<AliAODv0*>    ( fAODV0Refs   ->At(idxV0FromCascade) );
792       
793       // - Treatment of the parent for such a "re-used" V0 :
794       // Insert the cascade that reuses the V0 vertex in the lineage chain
795       // Before : vV0 -> vCascade1 -> vPrimary
796       //  - Hyp : cascade2 uses the same V0 as cascade1
797       //  After :  vV0 -> vCascade2 -> vCascade1 -> vPrimary
798       
799       AliAODVertex *vCascadePreviousParent = static_cast<AliAODVertex*> (vV0FromCascade->GetParent());
800                         vV0FromCascade->SetParent(vCascade);
801                         vCascade      ->SetParent(vCascadePreviousParent);
802       
803 //      if(fDebug > 2)  
804 //        printf("---- Cascade / Lineage insertion\n"
805 //               "Parent of V0 vtx                 = Cascade vtx %p\n"
806 //               "Parent of the cascade vtx        = Cascade vtx %p\n"
807 //               "Parent of the parent cascade vtx = Cascade vtx %p\n", 
808 //               static_cast<void*> (vV0FromCascade->GetParent()),
809 //               static_cast<void*> (vCascade->GetParent()),
810 //               static_cast<void*> (vCascadePreviousParent->GetParent()) );
811       
812     }// end if V0 structure already used
813     
814 //    if (fDebug > 2) {
815 //      printf("---- Cascade / V0 vertex: \n");
816 //      vV0FromCascade->Print();
817 //    }
818 //    
819 //    if (fDebug > 4) {
820 //      printf("---- Cascade / pos dghter : \n");
821 //                      aodTrack->Print();
822 //      printf("---- Cascade / neg dghter : \n");
823 //                      aodTrack->Print();
824 //      printf("---- Cascade / aodV0 : \n");
825 //                      aodV0->Print();
826 //    }
827     
828     // In any case (used V0 or not), add the V0 vertex to the cascade one.
829     vCascade->AddDaughter(vV0FromCascade);      
830     
831                 
832     // 5 - Add the primary track of the cascade (if any)
833     
834     
835     // 6 - Add the cascade to the AOD array of cascades
836     
837     Double_t dcaBachToPrimVertexXY = TMath::Abs(esdCascadeBach->GetD(esd.GetPrimaryVertex()->GetX(),
838                                                                      esd.GetPrimaryVertex()->GetY(),
839                                                                      esd.GetMagneticField())        );
840     
841     Double_t momBachAtCascadeVtx[3]={0.};
842
843     esdCascade->GetBPxPyPz(momBachAtCascadeVtx[0], momBachAtCascadeVtx[1], momBachAtCascadeVtx[2]);
844     
845     aodCascade = new(cascadesArray[fNumberOfCascades++]) AliAODcascade(  vCascade,
846                                                           esdCascade->Charge(),
847                                                           esdCascade->GetDcaXiDaughters(),
848                                                           -999.,
849                                                           // DCAXiToPrimVtx -> needs to be calculated   ----|
850                                                           // doesn't exist at ESD level;
851                                                           // See AODcascade::DcaXiToPrimVertex(Double, Double, Double)
852                                                           dcaBachToPrimVertexXY,
853                                                           momBachAtCascadeVtx,
854                                                           *aodV0);
855     
856     if (fDebug > 10) {
857       printf("---- Cascade / AOD cascade : \n\n");
858       aodCascade->PrintXi(fPrimaryVertex->GetX(), fPrimaryVertex->GetY(), fPrimaryVertex->GetZ());
859     }
860     
861   } // end of the loop on cascades
862   
863   Cascades().Expand(fNumberOfCascades);
864 }
865
866 //______________________________________________________________________________
867 void AliAnalysisTaskESDfilter::ConvertV0s(const AliESDEvent& esd)
868 {
869   // Access to the AOD container of V0s
870   
871   AliCodeTimerAuto("",0);
872
873   //
874   // V0s
875   //
876   
877   Double_t pos[3] = { 0. };      
878   Double_t chi2(0.0);
879   Double_t covVtx[6] = { 0. };
880   Double_t momPos[3]={0.};
881   Double_t covTr[21]={0.};
882   Double_t pid[10]={0.};
883   AliAODTrack* aodTrack(0x0);
884   AliAODPid* detpid(0x0);
885   Double_t momNeg[3]={0.};
886   Double_t momPosAtV0vtx[3]={0.};
887   Double_t momNegAtV0vtx[3]={0.};
888     
889   for (Int_t nV0 = 0; nV0 < esd.GetNumberOfV0s(); ++nV0) 
890   {
891     if (fUsedV0[nV0]) continue; // skip if already added to the AOD
892     
893     AliESDv0 *v0 = esd.GetV0(nV0);
894     Int_t posFromV0 = v0->GetPindex();
895     Int_t negFromV0 = v0->GetNindex();
896     
897     // V0 selection 
898     //
899     AliESDVertex *esdVtx   = new AliESDVertex(*(esd.GetPrimaryVertex()));
900     AliESDtrack  *esdV0Pos = esd.GetTrack(posFromV0);
901     AliESDtrack  *esdV0Neg = esd.GetTrack(negFromV0);
902     TList v0objects;
903     v0objects.AddAt(v0,                      0);
904     v0objects.AddAt(esdV0Pos,                1);
905     v0objects.AddAt(esdV0Neg,                2);
906     v0objects.AddAt(esdVtx,                  3);
907     UInt_t selectV0 = 0;
908     if (fV0Filter) {
909       selectV0 = fV0Filter->IsSelected(&v0objects);
910       // this is a little awkward but otherwise the 
911       // list wants to access the pointer (delete it) 
912       // again when going out of scope
913       delete v0objects.RemoveAt(3); // esdVtx created via copy construct
914       esdVtx = 0;
915       if (!selectV0) 
916         continue;
917     }
918     else{
919       delete v0objects.RemoveAt(3); // esdVtx created via copy construct
920       esdVtx = 0;
921     }
922     
923     v0->GetXYZ(pos[0], pos[1], pos[2]);
924     
925     if (!fOldESDformat) {
926             chi2 = v0->GetChi2V0(); // = chi2/NDF since NDF = 2*2-3
927             v0->GetPosCov(covVtx);
928     } else {
929             chi2 = -999.;
930             for (Int_t i = 0; i < 6; i++)  covVtx[i] = 0.;
931     }
932     
933     
934     AliAODVertex * vV0 = 
935           new(Vertices()[fNumberOfVertices++]) AliAODVertex(pos,
936                                                       covVtx,
937                                                       chi2,
938                                                       fPrimaryVertex,
939                                                       nV0,
940                                                       AliAODVertex::kV0);
941     fPrimaryVertex->AddDaughter(vV0);
942     
943     
944     // Add the positive tracks from the V0
945     
946
947     esdV0Pos->GetPxPyPz(momPos);
948     esdV0Pos->GetXYZ(pos);
949     esdV0Pos->GetCovarianceXYZPxPyPz(covTr);
950     esdV0Pos->GetESDpid(pid);
951     
952     const AliESDVertex *vtx = esd.GetPrimaryVertex();
953
954     if (!fUsedTrack[posFromV0]) {
955             fUsedTrack[posFromV0] = kTRUE;
956             UInt_t selectInfo = 0;
957             if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdV0Pos);
958             if(fMChandler)fMChandler->SelectParticle(esdV0Pos->GetLabel());
959             aodTrack = new(Tracks()[fNumberOfTracks++]) AliAODTrack(esdV0Pos->GetID(),
960                                                     esdV0Pos->GetLabel(), 
961                                                     momPos, 
962                                                     kTRUE,
963                                                     pos,
964                                                     kFALSE,
965                                                     covTr, 
966                                                     (Short_t)esdV0Pos->GetSign(),
967                                                     esdV0Pos->GetITSClusterMap(), 
968                                                     pid,
969                                                     vV0,
970                                                     kTRUE,  // check if this is right
971                                                     vtx->UsesTrack(esdV0Pos->GetID()),
972                                                     AliAODTrack::kSecondary,
973                                                     selectInfo);
974             aodTrack->SetTPCFitMap(esdV0Pos->GetTPCFitMap());
975             aodTrack->SetTPCClusterMap(esdV0Pos->GetTPCClusterMap());
976             aodTrack->SetTPCSharedMap (esdV0Pos->GetTPCSharedMap());
977             aodTrack->SetChi2perNDF(Chi2perNDF(esdV0Pos));
978             aodTrack->SetTPCPointsF(esdV0Pos->GetTPCNclsF());
979             aodTrack->SetTPCNCrossedRows(UShort_t(esdV0Pos->GetTPCCrossedRows()));
980             aodTrack->SetIntegratedLength(esdV0Pos->GetIntegratedLength());
981             fAODTrackRefs->AddAt(aodTrack,posFromV0);
982             //      if (fDebug > 0) printf("-------------------Bo: pos track from original pt %.3f \n",aodTrack->Pt());
983             if (esdV0Pos->GetSign() > 0) ++fNumberOfPositiveTracks;
984             aodTrack->ConvertAliPIDtoAODPID();
985             aodTrack->SetFlags(esdV0Pos->GetStatus());
986       SetAODPID(esdV0Pos,aodTrack,detpid);
987     }
988     else {
989             aodTrack = static_cast<AliAODTrack*>(fAODTrackRefs->At(posFromV0));
990             //      if (fDebug > 0) printf("-------------------Bo pos track from refArray pt %.3f \n",aodTrack->Pt());
991     }
992     vV0->AddDaughter(aodTrack);
993     
994     // Add the negative tracks from the V0
995     
996     esdV0Neg->GetPxPyPz(momNeg);
997     esdV0Neg->GetXYZ(pos);
998     esdV0Neg->GetCovarianceXYZPxPyPz(covTr);
999     esdV0Neg->GetESDpid(pid);
1000     
1001     if (!fUsedTrack[negFromV0]) {
1002             fUsedTrack[negFromV0] = kTRUE;
1003             UInt_t selectInfo = 0;
1004             if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdV0Neg);
1005             if(fMChandler)fMChandler->SelectParticle(esdV0Neg->GetLabel());
1006             aodTrack = new(Tracks()[fNumberOfTracks++]) AliAODTrack(esdV0Neg->GetID(),
1007                                                     esdV0Neg->GetLabel(),
1008                                                     momNeg,
1009                                                     kTRUE,
1010                                                     pos,
1011                                                     kFALSE,
1012                                                     covTr, 
1013                                                     (Short_t)esdV0Neg->GetSign(),
1014                                                     esdV0Neg->GetITSClusterMap(), 
1015                                                     pid,
1016                                                     vV0,
1017                                                     kTRUE,  // check if this is right
1018                                                     vtx->UsesTrack(esdV0Neg->GetID()),
1019                                                     AliAODTrack::kSecondary,
1020                                                     selectInfo);
1021             aodTrack->SetTPCFitMap(esdV0Neg->GetTPCFitMap());
1022             aodTrack->SetTPCClusterMap(esdV0Neg->GetTPCClusterMap());
1023             aodTrack->SetTPCSharedMap (esdV0Neg->GetTPCSharedMap());
1024             aodTrack->SetChi2perNDF(Chi2perNDF(esdV0Neg));
1025             aodTrack->SetTPCPointsF(esdV0Neg->GetTPCNclsF());
1026             aodTrack->SetTPCNCrossedRows(UShort_t(esdV0Neg->GetTPCCrossedRows()));
1027             aodTrack->SetIntegratedLength(esdV0Neg->GetIntegratedLength());
1028             fAODTrackRefs->AddAt(aodTrack,negFromV0);
1029             //      if (fDebug > 0) printf("-------------------Bo: neg track from original pt %.3f \n",aodTrack->Pt());
1030             if (esdV0Neg->GetSign() > 0) ++fNumberOfPositiveTracks;
1031             aodTrack->ConvertAliPIDtoAODPID();
1032             aodTrack->SetFlags(esdV0Neg->GetStatus());
1033       SetAODPID(esdV0Neg,aodTrack,detpid);
1034     }
1035     else {
1036             aodTrack = static_cast<AliAODTrack*>(fAODTrackRefs->At(negFromV0));
1037             //      if (fDebug > 0) printf("-------------------Bo neg track from refArray pt %.3f \n",aodTrack->Pt());
1038     }
1039     vV0->AddDaughter(aodTrack);
1040     
1041     
1042     // Add the V0 the V0 array as well
1043     
1044     Double_t  dcaV0Daughters      = v0->GetDcaV0Daughters();
1045     Double_t  dcaV0ToPrimVertex   = v0->GetD(esd.GetPrimaryVertex()->GetX(),
1046                                              esd.GetPrimaryVertex()->GetY(),
1047                                              esd.GetPrimaryVertex()->GetZ());
1048     
1049     v0->GetPPxPyPz(momPosAtV0vtx[0],momPosAtV0vtx[1],momPosAtV0vtx[2]); 
1050     v0->GetNPxPyPz(momNegAtV0vtx[0],momNegAtV0vtx[1],momNegAtV0vtx[2]); 
1051     
1052     Double_t dcaDaughterToPrimVertex[2] = { 999., 999.}; // ..[0] = DCA in (x,y) for Pos and ..[1] = Neg
1053     dcaDaughterToPrimVertex[0] = TMath::Abs(esdV0Pos->GetD(  esd.GetPrimaryVertex()->GetX(),
1054                                                            esd.GetPrimaryVertex()->GetY(),
1055                                                            esd.GetMagneticField()) );
1056     dcaDaughterToPrimVertex[1] = TMath::Abs(esdV0Neg->GetD(  esd.GetPrimaryVertex()->GetX(),
1057                                                            esd.GetPrimaryVertex()->GetY(),
1058                                                            esd.GetMagneticField()) );
1059     
1060     AliAODv0* aodV0 = new(V0s()[fNumberOfV0s++]) AliAODv0(vV0, 
1061                                                 dcaV0Daughters,
1062                                                 dcaV0ToPrimVertex,
1063                                                 momPosAtV0vtx,
1064                                                 momNegAtV0vtx,
1065                                                 dcaDaughterToPrimVertex);
1066     
1067     // set the aod v0 on-the-fly status
1068     aodV0->SetOnFlyStatus(v0->GetOnFlyStatus());
1069   }//End of loop on V0s 
1070   
1071   V0s().Expand(fNumberOfV0s);    
1072 }
1073
1074 //______________________________________________________________________________
1075 void AliAnalysisTaskESDfilter::ConvertTPCOnlyTracks(const AliESDEvent& esd)
1076 {
1077   // Convert TPC only tracks
1078   // Here we have wo hybrid appraoch to remove fakes
1079   // ******* ITSTPC ********
1080   // Uses a cut on the ITS properties to select global tracks
1081   // which are than marked as HybdridITSTPC for the remainder 
1082   // the TPC only tracks are flagged as HybridITSTPConly. 
1083   // Note, in order not to get fakes back in the TPC cuts, one needs 
1084   // two "ITS" cuts one tight (1) (to throw out fakes) and one lose (2) (to NOT flag the trakcs in the TPC only)
1085   // using cut number (3)
1086   // so fHybridFilterMask == (1)|(2) fTPCFilterMask = (3), Usercode needs to slect with mask = (1)|(3) and track->IsHybridITSTPC()
1087   // ******* TPC ********
1088   // Here, only TPC tracks are flagged that pass the tight ITS cuts and tracks that pass the TPC cuts and NOT the loose ITS cuts
1089   // the ITS cuts neeed to be added to the filter as extra cuts, since here the selections info is reset in the global and put to the TPC only track
1090
1091   AliCodeTimerAuto("",0);
1092   
1093   // Loop over the tracks and extract and mask out all aod tracks that pass the selections for AODt racks
1094   for(int it = 0;it < fNumberOfTracks;++it)
1095   {
1096     AliAODTrack *tr = (AliAODTrack*)(Tracks().At(it));
1097     if(!tr)continue;
1098     UInt_t map = tr->GetFilterMap();
1099     if(map&fTPCConstrainedFilterMask){
1100       // we only reset the track select ionfo, no deletion...
1101       tr->SetFilterMap(map&~fTPCConstrainedFilterMask);
1102     }
1103     if(map&fHybridFilterMaskTPCCG){
1104       // this is one part of the hybrid tracks
1105       // the others not passing the selection will be TPC only selected below
1106       tr->SetIsHybridTPCConstrainedGlobal(kTRUE);
1107     }
1108   }
1109   // Loop over the ESD trcks and pick out the tracks passing TPC only cuts
1110   
1111   
1112   const AliESDVertex *vtxSPD = esd.GetPrimaryVertexSPD();
1113   const AliESDVertex *vtx = esd.GetPrimaryVertex();
1114
1115   Double_t pos[3] = { 0. };      
1116   Double_t covTr[21]={0.};
1117   Double_t pid[10]={0.};  
1118
1119
1120   Double_t p[3] = { 0. };
1121
1122   Double_t pDCA[3] = { 0. }; // momentum at DCA
1123   Double_t rDCA[3] = { 0. }; // position at DCA
1124   Float_t  dDCA[2] = {0.};    // DCA to the vertex d and z
1125   Float_t  cDCA[3] = {0.};    // covariance of impact parameters
1126
1127
1128   AliAODTrack* aodTrack(0x0);
1129   //  AliAODPid* detpid(0x0);
1130
1131   // account for change in pT after the constraint
1132   Float_t ptMax = 1E10;
1133   Float_t ptMin = 0;
1134   for(int i = 0;i<32;i++){
1135     if(fTPCConstrainedFilterMask&(1<<i)){
1136       AliESDtrackCuts*cuts = (AliESDtrackCuts*)fTrackFilter->GetCuts()->At(i);
1137       Float_t tmp1= 0,tmp2 = 0;
1138       cuts->GetPtRange(tmp1,tmp2);
1139       if(tmp1>ptMin)ptMin=tmp1;
1140       if(tmp2<ptMax)ptMax=tmp2;
1141     }
1142   } 
1143
1144   for (Int_t nTrack = 0; nTrack < esd.GetNumberOfTracks(); ++nTrack) 
1145   {
1146     AliESDtrack* esdTrack = esd.GetTrack(nTrack); //carefull do not modify it othwise  need to work with a copy 
1147     
1148     UInt_t selectInfo = 0;
1149     Bool_t isHybridITSTPC = false;
1150     //
1151     // Track selection
1152     if (fTrackFilter) {
1153       selectInfo = fTrackFilter->IsSelected(esdTrack);
1154     }
1155
1156     if(!(selectInfo&fHybridFilterMaskTPCCG)){
1157       // not already selected tracks, use second part of hybrid tracks
1158       isHybridITSTPC = true;
1159       // too save space one could only store these...
1160     }
1161
1162     selectInfo &= fTPCConstrainedFilterMask;
1163     if (!selectInfo)continue;
1164     if (fWriteHybridTPCCOnly&&!isHybridITSTPC)continue; // write only complementary tracks
1165     // create a tpc only tracl
1166     AliESDtrack *track = AliESDtrackCuts::GetTPCOnlyTrack(const_cast<AliESDEvent*>(&esd),esdTrack->GetID());
1167     if(!track) continue;
1168     
1169     if(track->Pt()>0.)
1170     {
1171       // only constrain tracks above threshold
1172       AliExternalTrackParam exParam;
1173       // take the B-field from the ESD, no 3D fieldMap available at this point
1174       Bool_t relate = false;
1175       relate = track->RelateToVertexTPC(vtxSPD,esd.GetMagneticField(),kVeryBig,&exParam);
1176       if(!relate){
1177         delete track;
1178         continue;
1179       }
1180       // fetch the track parameters at the DCA (unconstraint)
1181       if(track->GetTPCInnerParam()){
1182         track->GetTPCInnerParam()->GetPxPyPz(pDCA);
1183         track->GetTPCInnerParam()->GetXYZ(rDCA);
1184       }
1185       // get the DCA to the vertex:
1186       track->GetImpactParametersTPC(dDCA,cDCA);
1187       // set the constrained parameters to the track
1188       track->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
1189     }
1190     
1191     track->GetPxPyPz(p);
1192
1193     Float_t pT = track->Pt();
1194     if(pT<ptMin||pT>ptMax){
1195       delete track;
1196       continue;
1197     }
1198
1199     // 
1200
1201
1202     track->GetXYZ(pos);
1203     track->GetCovarianceXYZPxPyPz(covTr);
1204     esdTrack->GetESDpid(pid);// original PID
1205
1206     if(fMChandler)fMChandler->SelectParticle(esdTrack->GetLabel());
1207     aodTrack = new(Tracks()[fNumberOfTracks++]) AliAODTrack((track->GetID()+1)*-1,
1208                                                             track->GetLabel(),
1209                                                             p,
1210                                                             kTRUE,
1211                                                             pos,
1212                                                             kFALSE,
1213                                                             covTr, 
1214                                                             (Short_t)track->GetSign(),
1215                                                             track->GetITSClusterMap(), 
1216                                                             pid,
1217                                                             fPrimaryVertex,
1218                                                             kTRUE, // check if this is right
1219                                                             vtx->UsesTrack(track->GetID()),
1220                                                             AliAODTrack::kPrimary, 
1221                                                             selectInfo);
1222     aodTrack->SetIsHybridTPCConstrainedGlobal(isHybridITSTPC);    
1223     aodTrack->SetTPCFitMap(track->GetTPCFitMap());
1224     aodTrack->SetTPCClusterMap(track->GetTPCClusterMap());
1225     aodTrack->SetTPCSharedMap (track->GetTPCSharedMap());
1226     aodTrack->SetIsTPCConstrained(kTRUE);    
1227     aodTrack->SetChi2perNDF(Chi2perNDF(esdTrack)); // original track
1228     // set the DCA values to the AOD track
1229     aodTrack->SetPxPyPzAtDCA(pDCA[0],pDCA[1],pDCA[2]);
1230     aodTrack->SetXYAtDCA(rDCA[0],rDCA[1]);
1231     aodTrack->SetDCA(dDCA[0],dDCA[1]);
1232
1233     aodTrack->SetFlags(track->GetStatus());
1234     aodTrack->SetTPCPointsF(track->GetTPCNclsF());
1235     aodTrack->SetTPCNCrossedRows(UShort_t(track->GetTPCCrossedRows()));
1236     aodTrack->SetIntegratedLength(track->GetIntegratedLength());
1237
1238     //Perform progagation of tracks if needed
1239     if(fDoPropagateTrackToEMCal) PropagateTrackToEMCal(esdTrack);
1240     aodTrack->SetTrackPhiEtaPtOnEMCal(esdTrack->GetTrackPhiOnEMCal(),esdTrack->GetTrackEtaOnEMCal(),esdTrack->GetTrackPtOnEMCal());
1241
1242     // do not duplicate PID information 
1243     //    aodTrack->ConvertAliPIDtoAODPID();
1244     //    SetAODPID(esdTrack,aodTrack,detpid);
1245
1246     delete track;
1247   } // end of loop on tracks
1248   
1249 }
1250
1251 //______________________________________________________________________________
1252 void AliAnalysisTaskESDfilter::ConvertGlobalConstrainedTracks(const AliESDEvent& esd)
1253 {
1254
1255   // Here we have the option to store the complement from global constraint information
1256   // to tracks passing tight cuts (1) in order not to get fakes back in, one needs 
1257   // two sets of cuts one tight (1) (to throw out fakes) and one lose (2) (fakes/bad tracks would pass (2) but not (1))
1258   // using cut number (3) selects the tracks that complement (1) e.g. tracks witout ITS refit or cluster requirement
1259
1260
1261   AliCodeTimerAuto("",0);
1262   
1263   // Loop over the tracks and extract and mask out all aod tracks that pass the selections for AODt racks
1264   for(int it = 0;it < fNumberOfTracks;++it)
1265   {
1266     AliAODTrack *tr = (AliAODTrack*)(Tracks().At(it));
1267     if(!tr)continue;
1268     UInt_t map = tr->GetFilterMap();
1269     if(map&fGlobalConstrainedFilterMask){
1270       // we only reset the track select info, no deletion...
1271       // mask reset mask in case track is already taken
1272       tr->SetFilterMap(map&~fGlobalConstrainedFilterMask);
1273     }
1274     if(map&fHybridFilterMaskGCG){
1275       // this is one part of the hybrid tracks
1276       // the others not passing the selection will be the ones selected below
1277       tr->SetIsHybridGlobalConstrainedGlobal(kTRUE);
1278     }
1279   }
1280   // Loop over the ESD trcks and pick out the tracks passing the GlobalConstraint cuts
1281  
1282
1283   Double_t pos[3] = { 0. };      
1284   Double_t covTr[21]={0.};
1285   Double_t pid[10]={0.};  
1286   Double_t p[3] = { 0. };
1287
1288   Double_t pDCA[3] = { 0. }; // momentum at DCA
1289   Double_t rDCA[3] = { 0. }; // position at DCA
1290   Float_t  dDCA[2] = {0.};    // DCA to the vertex d and z
1291   Float_t  cDCA[3] = {0.};    // covariance of impact parameters
1292
1293
1294   AliAODTrack* aodTrack(0x0);
1295   AliAODPid* detpid(0x0);
1296   const AliESDVertex *vtx = esd.GetPrimaryVertex();
1297
1298   // account for change in pT after the constraint
1299   Float_t ptMax = 1E10;
1300   Float_t ptMin = 0;
1301   for(int i = 0;i<32;i++){
1302     if(fGlobalConstrainedFilterMask&(1<<i)){
1303       AliESDtrackCuts*cuts = (AliESDtrackCuts*)fTrackFilter->GetCuts()->At(i);
1304       Float_t tmp1= 0,tmp2 = 0;
1305       cuts->GetPtRange(tmp1,tmp2);
1306       if(tmp1>ptMin)ptMin=tmp1;
1307       if(tmp2<ptMax)ptMax=tmp2;
1308     }
1309   } 
1310
1311
1312
1313  for (Int_t nTrack = 0; nTrack < esd.GetNumberOfTracks(); ++nTrack) 
1314   {
1315     AliESDtrack* esdTrack = esd.GetTrack(nTrack); //carefull do not modify it othwise  need to work with a copy 
1316     const AliExternalTrackParam * exParamGC = esdTrack->GetConstrainedParam();
1317     if(!exParamGC)continue;
1318
1319     UInt_t selectInfo = 0;
1320     Bool_t isHybridGC = false;
1321
1322     //
1323     // Track selection
1324     if (fTrackFilter) {
1325       selectInfo = fTrackFilter->IsSelected(esdTrack);
1326     }
1327
1328
1329     if(!(selectInfo&fHybridFilterMaskGCG))isHybridGC = true;
1330     if (fWriteHybridGCOnly&&!isHybridGC)continue; // write only complementary tracks
1331
1332     selectInfo &= fGlobalConstrainedFilterMask;
1333     if (!selectInfo)continue;
1334     // fetch the track parameters at the DCA (unconstrained)
1335     esdTrack->GetPxPyPz(pDCA);
1336     esdTrack->GetXYZ(rDCA);
1337     // get the DCA to the vertex:
1338     esdTrack->GetImpactParameters(dDCA,cDCA);
1339
1340     if (!esdTrack->GetConstrainedPxPyPz(p)) continue;
1341
1342
1343     Float_t pT = exParamGC->Pt();
1344     if(pT<ptMin||pT>ptMax){
1345       continue;
1346     }
1347
1348
1349     esdTrack->GetConstrainedXYZ(pos);
1350     exParamGC->GetCovarianceXYZPxPyPz(covTr);
1351     esdTrack->GetESDpid(pid);
1352     if(fMChandler)fMChandler->SelectParticle(esdTrack->GetLabel());
1353     aodTrack = new(Tracks()[fNumberOfTracks++]) AliAODTrack((esdTrack->GetID()+1)*-1,
1354                                                             esdTrack->GetLabel(),
1355                                                             p,
1356                                                             kTRUE,
1357                                                             pos,
1358                                                             kFALSE,
1359                                                             covTr, 
1360                                                             (Short_t)esdTrack->GetSign(),
1361                                                             esdTrack->GetITSClusterMap(), 
1362                                                             pid,
1363                                                             fPrimaryVertex,
1364                                                             kTRUE, // check if this is right
1365                                                             vtx->UsesTrack(esdTrack->GetID()),
1366                                                             AliAODTrack::kPrimary, 
1367                                                             selectInfo);
1368     aodTrack->SetIsHybridGlobalConstrainedGlobal(isHybridGC);    
1369     aodTrack->SetIsGlobalConstrained(kTRUE);    
1370     aodTrack->SetTPCFitMap(esdTrack->GetTPCFitMap());
1371     aodTrack->SetTPCClusterMap(esdTrack->GetTPCClusterMap());
1372     aodTrack->SetTPCSharedMap (esdTrack->GetTPCSharedMap());
1373     aodTrack->SetChi2perNDF(Chi2perNDF(esdTrack));
1374
1375
1376     // set the DCA values to the AOD track
1377     aodTrack->SetPxPyPzAtDCA(pDCA[0],pDCA[1],pDCA[2]);
1378     aodTrack->SetXYAtDCA(rDCA[0],rDCA[1]);
1379     aodTrack->SetDCA(dDCA[0],dDCA[1]);
1380
1381     aodTrack->SetFlags(esdTrack->GetStatus());
1382     aodTrack->SetTPCPointsF(esdTrack->GetTPCNclsF());
1383     aodTrack->SetTPCNCrossedRows(UShort_t(esdTrack->GetTPCCrossedRows()));
1384     aodTrack->SetIntegratedLength(esdTrack->GetIntegratedLength());
1385
1386     if(isHybridGC){
1387       // only copy AOD information for hybrid, no duplicate information
1388       aodTrack->ConvertAliPIDtoAODPID();
1389       SetAODPID(esdTrack,aodTrack,detpid);
1390     }
1391
1392     //Perform progagation of tracks if needed
1393     if(fDoPropagateTrackToEMCal) PropagateTrackToEMCal(esdTrack);
1394     aodTrack->SetTrackPhiEtaPtOnEMCal(esdTrack->GetTrackPhiOnEMCal(),esdTrack->GetTrackEtaOnEMCal(),esdTrack->GetTrackPtOnEMCal());
1395   } // end of loop on tracks
1396   
1397 }
1398
1399
1400 //______________________________________________________________________________
1401 void AliAnalysisTaskESDfilter::ConvertTracks(const AliESDEvent& esd)
1402 {
1403   // Tracks (primary and orphan)
1404
1405   AliCodeTimerAuto("",0);
1406   
1407   AliDebug(1,Form("NUMBER OF ESD TRACKS %5d\n", esd.GetNumberOfTracks()));
1408   
1409   const AliESDVertex *vtx = esd.GetPrimaryVertex();
1410   Double_t p[3] = { 0. };
1411   Double_t pos[3] = { 0. };
1412   Double_t covTr[21] = { 0. };
1413   Double_t pid[10] = { 0. };
1414   AliAODTrack* aodTrack(0x0);
1415   AliAODPid* detpid(0x0);
1416   
1417   for (Int_t nTrack = 0; nTrack < esd.GetNumberOfTracks(); ++nTrack) 
1418   {
1419     if (fUsedTrack[nTrack]) continue;
1420     
1421     AliESDtrack *esdTrack = esd.GetTrack(nTrack);
1422     UInt_t selectInfo = 0;
1423     //
1424     // Track selection
1425     if (fTrackFilter) {
1426             selectInfo = fTrackFilter->IsSelected(esdTrack);
1427             if (!selectInfo && !vtx->UsesTrack(esdTrack->GetID())) continue;
1428     }
1429     
1430     
1431     esdTrack->GetPxPyPz(p);
1432     esdTrack->GetXYZ(pos);
1433     esdTrack->GetCovarianceXYZPxPyPz(covTr);
1434     esdTrack->GetESDpid(pid);
1435     if(fMChandler)fMChandler->SelectParticle(esdTrack->GetLabel());
1436     fPrimaryVertex->AddDaughter(aodTrack =
1437                          new(Tracks()[fNumberOfTracks++]) AliAODTrack(esdTrack->GetID(),
1438                                                             esdTrack->GetLabel(),
1439                                                             p,
1440                                                             kTRUE,
1441                                                             pos,
1442                                                             kFALSE,
1443                                                             covTr, 
1444                                                             (Short_t)esdTrack->GetSign(),
1445                                                             esdTrack->GetITSClusterMap(), 
1446                                                             pid,
1447                                                             fPrimaryVertex,
1448                                                             kTRUE, // check if this is right
1449                                                             vtx->UsesTrack(esdTrack->GetID()),
1450                                                             AliAODTrack::kPrimary, 
1451                                                             selectInfo)
1452                          );
1453     aodTrack->SetTPCFitMap(esdTrack->GetTPCFitMap());
1454     aodTrack->SetTPCClusterMap(esdTrack->GetTPCClusterMap());
1455     aodTrack->SetTPCSharedMap (esdTrack->GetTPCSharedMap());
1456     aodTrack->SetChi2perNDF(Chi2perNDF(esdTrack));
1457     aodTrack->SetTPCPointsF(esdTrack->GetTPCNclsF());
1458     aodTrack->SetTPCNCrossedRows(UShort_t(esdTrack->GetTPCCrossedRows()));
1459     aodTrack->SetIntegratedLength(esdTrack->GetIntegratedLength());
1460     if(esdTrack->IsEMCAL()) aodTrack->SetEMCALcluster(esdTrack->GetEMCALcluster());
1461     if(esdTrack->IsPHOS())  aodTrack->SetPHOScluster(esdTrack->GetPHOScluster());
1462
1463     //Perform progagation of tracks if needed
1464     if(fDoPropagateTrackToEMCal) PropagateTrackToEMCal(esdTrack);
1465     aodTrack->SetTrackPhiEtaPtOnEMCal(esdTrack->GetTrackPhiOnEMCal(),esdTrack->GetTrackEtaOnEMCal(),esdTrack->GetTrackPtOnEMCal());
1466
1467     fAODTrackRefs->AddAt(aodTrack, nTrack);
1468     
1469     
1470     if (esdTrack->GetSign() > 0) ++fNumberOfPositiveTracks;
1471     aodTrack->SetFlags(esdTrack->GetStatus());
1472     aodTrack->ConvertAliPIDtoAODPID();
1473     SetAODPID(esdTrack,aodTrack,detpid);
1474   } // end of loop on tracks
1475 }
1476
1477 //______________________________________________________________________________
1478 void AliAnalysisTaskESDfilter::PropagateTrackToEMCal(AliESDtrack *esdTrack)
1479 {
1480   Double_t trkPos[3] = {0.,0.,0.};
1481   Double_t EMCalEta=-999, EMCalPhi=-999, EMCalPt=-999;
1482   Double_t trkphi = esdTrack->Phi()*TMath::RadToDeg();
1483   if(TMath::Abs(esdTrack->Eta())<0.9 && trkphi > 10 && trkphi < 250 )
1484     {
1485       AliExternalTrackParam *trkParam = const_cast<AliExternalTrackParam*>(esdTrack->GetInnerParam());
1486       if(trkParam)
1487         {
1488           AliExternalTrackParam trkParamTmp(*trkParam);
1489           if(AliTrackerBase::PropagateTrackToBxByBz(&trkParamTmp, fEMCalSurfaceDistance, esdTrack->GetMass(), 20, kTRUE, 0.8, -1))
1490             {
1491               trkParamTmp.GetXYZ(trkPos);
1492               TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
1493               EMCalEta = trkPosVec.Eta();
1494               EMCalPhi = trkPosVec.Phi();
1495               EMCalPt = trkParamTmp.Pt();
1496               if(EMCalPhi<0)  EMCalPhi += 2*TMath::Pi();
1497               esdTrack->SetTrackPhiEtaPtOnEMCal(EMCalPhi,EMCalEta,EMCalPt);
1498             }
1499         }
1500     }
1501 }
1502
1503 //______________________________________________________________________________
1504 void AliAnalysisTaskESDfilter::ConvertPmdClusters(const AliESDEvent& esd)
1505 {
1506 // Convert PMD Clusters 
1507   AliCodeTimerAuto("",0);
1508   Int_t jPmdClusters=0;
1509   // Access to the AOD container of PMD clusters
1510   TClonesArray &pmdClusters = *(AODEvent()->GetPmdClusters());
1511   for (Int_t iPmd = 0; iPmd < esd.GetNumberOfPmdTracks(); ++iPmd) {
1512     // file pmd clusters, to be revised!
1513     AliESDPmdTrack *pmdTrack = esd.GetPmdTrack(iPmd);
1514     Int_t nLabel = 0;
1515     Int_t *label = 0x0;
1516     Double_t posPmd[3] = { pmdTrack->GetClusterX(), pmdTrack->GetClusterY(), pmdTrack->GetClusterZ()};
1517     Double_t pidPmd[13] = { 0.}; // to be revised!
1518     // type not set!
1519     // assoc cluster not set
1520     new(pmdClusters[jPmdClusters++]) AliAODPmdCluster(iPmd, nLabel, label, pmdTrack->GetClusterADC(), posPmd, pidPmd);
1521   }
1522 }
1523
1524
1525 //______________________________________________________________________________
1526 void AliAnalysisTaskESDfilter::ConvertCaloClusters(const AliESDEvent& esd)
1527 {
1528 // Convert Calorimeter Clusters
1529   AliCodeTimerAuto("",0);
1530   
1531   // Access to the AOD container of clusters
1532   TClonesArray &caloClusters = *(AODEvent()->GetCaloClusters());
1533   Int_t jClusters(0);
1534   
1535   for (Int_t iClust=0; iClust<esd.GetNumberOfCaloClusters(); ++iClust) {
1536     
1537     AliESDCaloCluster * cluster = esd.GetCaloCluster(iClust);
1538     
1539     Int_t  id        = cluster->GetID();
1540     Int_t  nLabel    = cluster->GetNLabels();
1541     Int_t *labels    = cluster->GetLabels();
1542     if(labels){ 
1543                   for(int i = 0;i < nLabel;++i){
1544                           if(fMChandler)fMChandler->SelectParticle(labels[i]);
1545                   }
1546           }             
1547     
1548     Float_t energy = cluster->E();
1549     Float_t posF[3] = { 0.};
1550     cluster->GetPosition(posF);
1551     
1552     AliAODCaloCluster *caloCluster = new(caloClusters[jClusters++]) AliAODCaloCluster(id,
1553                                                                                       nLabel,
1554                                                                                       labels,
1555                                                                                       energy,
1556                                                                                       posF,
1557                                                                                       NULL,
1558                                                                                       cluster->GetType(),0);
1559     
1560     caloCluster->SetCaloCluster(cluster->GetDistanceToBadChannel(),
1561                                 cluster->GetDispersion(),
1562                                 cluster->GetM20(), cluster->GetM02(),
1563                                 cluster->GetEmcCpvDistance(),  
1564                                 cluster->GetNExMax(),cluster->GetTOF()) ;
1565     
1566     caloCluster->SetPIDFromESD(cluster->GetPID());
1567     caloCluster->SetNCells(cluster->GetNCells());
1568     caloCluster->SetCellsAbsId(cluster->GetCellsAbsId());
1569     caloCluster->SetCellsAmplitudeFraction(cluster->GetCellsAmplitudeFraction());
1570
1571     caloCluster->SetTrackDistance(cluster->GetTrackDx(), cluster->GetTrackDz());
1572     
1573     Int_t nMatchCount = 0;
1574     TArrayI* matchedT =         cluster->GetTracksMatched();
1575     if (fNumberOfTracks>0 && matchedT && cluster->GetTrackMatchedIndex() >= 0) {        
1576       for (Int_t im = 0; im < matchedT->GetSize(); im++) {
1577         Int_t iESDtrack = matchedT->At(im);;
1578         if (fAODTrackRefs->At(iESDtrack) != 0) {
1579           caloCluster->AddTrackMatched((AliAODTrack*)fAODTrackRefs->At(iESDtrack));
1580           nMatchCount++;
1581         }
1582       }
1583     }
1584     if(nMatchCount==0)
1585       caloCluster->SetTrackDistance(-999,-999);
1586     
1587   } 
1588   caloClusters.Expand(jClusters); // resize TObjArray to 'remove' slots for pseudo clusters      
1589 }
1590
1591 //______________________________________________________________________________
1592 void AliAnalysisTaskESDfilter::ConvertCaloTrigger(TString calo, const AliESDEvent& esd)
1593 {
1594         AliCodeTimerAuto("",0);
1595                 
1596         if (calo == "PHOS") 
1597         {
1598           AliAODCaloTrigger &aodTrigger = *(AODEvent()->GetCaloTrigger(calo));
1599           AliESDCaloTrigger &esdTrigger = *(esd.GetCaloTrigger(calo));
1600
1601           aodTrigger.Allocate(esdTrigger.GetEntries());
1602           esdTrigger.Reset();
1603
1604           Float_t a;
1605           Int_t tmod,tabsId;
1606
1607           while (esdTrigger.Next()) {
1608             esdTrigger.GetPosition(tmod,tabsId);
1609             esdTrigger.GetAmplitude(a);
1610             aodTrigger.Add(tmod,tabsId,a,0.,(Int_t*)NULL,0,0,0);
1611           }
1612
1613           return;
1614         }
1615                         
1616         AliAODHandler *aodHandler = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()); 
1617                         
1618         if (aodHandler)
1619         {
1620                 TTree *aodTree = aodHandler->GetTree();
1621                                         
1622                 if (aodTree)
1623                 {
1624                         Int_t *type = esd.GetCaloTriggerType();
1625                                                         
1626                         for (Int_t i = 0; i < 15; i++) 
1627                         {
1628                                 aodTree->GetUserInfo()->Add(new TParameter<int>(Form("EMCALCaloTrigger%d",i), type[i]));
1629                         }
1630                 }
1631         }
1632                                                 
1633         AliAODCaloTrigger &aodTrigger = *(AODEvent()->GetCaloTrigger(calo));
1634                                                 
1635         AliESDCaloTrigger &esdTrigger = *(esd.GetCaloTrigger(calo));
1636                                                 
1637         aodTrigger.Allocate(esdTrigger.GetEntries());
1638                                                 
1639         esdTrigger.Reset();
1640         while (esdTrigger.Next())
1641         {         
1642                 Int_t px, py, ts, nTimes, times[10], b; 
1643                 Float_t a, t;
1644                                                                 
1645                 esdTrigger.GetPosition(px, py);
1646                                                 
1647                 esdTrigger.GetAmplitude(a);
1648                 esdTrigger.GetTime(t);
1649                                                                 
1650                 esdTrigger.GetL0Times(times);
1651                 esdTrigger.GetNL0Times(nTimes);
1652                                                                 
1653                 esdTrigger.GetL1TimeSum(ts);
1654                                                                 
1655                 esdTrigger.GetTriggerBits(b);
1656                                                                 
1657                 aodTrigger.Add(px, py, a, t, times, nTimes, ts, b);
1658         }
1659                                                         
1660         for (int i = 0; i < 4; i++) aodTrigger.SetL1Threshold(i, esdTrigger.GetL1Threshold(i));
1661                                                         
1662         Int_t v0[2] = 
1663         {
1664                 esdTrigger.GetL1V0(0),
1665                 esdTrigger.GetL1V0(1)
1666         };              
1667                                                                 
1668         aodTrigger.SetL1V0(v0); 
1669         aodTrigger.SetL1FrameMask(esdTrigger.GetL1FrameMask());
1670 }
1671
1672 //______________________________________________________________________________
1673 void AliAnalysisTaskESDfilter::ConvertEMCALCells(const AliESDEvent& esd)
1674 {
1675 // Convert EMCAL Cells
1676   AliCodeTimerAuto("",0);
1677   // fill EMCAL cell info
1678   if (esd.GetEMCALCells()) { // protection against missing ESD information
1679     AliESDCaloCells &esdEMcells = *(esd.GetEMCALCells());
1680     Int_t nEMcell = esdEMcells.GetNumberOfCells() ;
1681     
1682     AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
1683     aodEMcells.CreateContainer(nEMcell);
1684     aodEMcells.SetType(AliAODCaloCells::kEMCALCell);
1685     for (Int_t iCell = 0; iCell < nEMcell; iCell++) {      
1686       aodEMcells.SetCell(iCell,esdEMcells.GetCellNumber(iCell),esdEMcells.GetAmplitude(iCell),
1687                          esdEMcells.GetTime(iCell), esdEMcells.GetMCLabel(iCell), esdEMcells.GetEFraction(iCell));
1688     }
1689     aodEMcells.Sort();
1690   }
1691 }
1692
1693 //______________________________________________________________________________
1694 void AliAnalysisTaskESDfilter::ConvertPHOSCells(const AliESDEvent& esd)
1695 {
1696 // Convert PHOS Cells
1697   AliCodeTimerAuto("",0);
1698   // fill PHOS cell info
1699   if (esd.GetPHOSCells()) { // protection against missing ESD information
1700     AliESDCaloCells &esdPHcells = *(esd.GetPHOSCells());
1701     Int_t nPHcell = esdPHcells.GetNumberOfCells() ;
1702     
1703     AliAODCaloCells &aodPHcells = *(AODEvent()->GetPHOSCells());
1704     aodPHcells.CreateContainer(nPHcell);
1705     aodPHcells.SetType(AliAODCaloCells::kPHOSCell);
1706     for (Int_t iCell = 0; iCell < nPHcell; iCell++) {      
1707       aodPHcells.SetCell(iCell,esdPHcells.GetCellNumber(iCell),esdPHcells.GetAmplitude(iCell),
1708                          esdPHcells.GetTime(iCell), esdPHcells.GetMCLabel(iCell), esdPHcells.GetEFraction(iCell));
1709     }
1710     aodPHcells.Sort();
1711   }
1712 }
1713
1714 //______________________________________________________________________________
1715 void AliAnalysisTaskESDfilter::ConvertTracklets(const AliESDEvent& esd)
1716 {
1717   // tracklets    
1718   AliCodeTimerAuto("",0);
1719
1720   AliAODTracklets &SPDTracklets = *(AODEvent()->GetTracklets());
1721   const AliMultiplicity *mult = esd.GetMultiplicity();
1722   if (mult) {
1723     if (mult->GetNumberOfTracklets()>0) {
1724       SPDTracklets.CreateContainer(mult->GetNumberOfTracklets());
1725       
1726       for (Int_t n=0; n<mult->GetNumberOfTracklets(); n++) {
1727         if(fMChandler){
1728           fMChandler->SelectParticle(mult->GetLabel(n, 0));
1729           fMChandler->SelectParticle(mult->GetLabel(n, 1));
1730         }
1731         SPDTracklets.SetTracklet(n, mult->GetTheta(n), mult->GetPhi(n), mult->GetDeltaPhi(n), mult->GetLabel(n, 0),mult->GetLabel(n, 1));
1732       }
1733     }
1734   } else {
1735     //Printf("ERROR: AliMultiplicity could not be retrieved from ESD");
1736   }
1737 }
1738
1739 //______________________________________________________________________________
1740 void AliAnalysisTaskESDfilter::ConvertKinks(const AliESDEvent& esd)
1741 {
1742   AliCodeTimerAuto("",0);
1743   
1744   // Kinks: it is a big mess the access to the information in the kinks
1745   // The loop is on the tracks in order to find the mother and daugther of each kink
1746   
1747   Double_t covTr[21]={0.};
1748   Double_t pid[10]={0.};
1749   AliAODPid* detpid(0x0);
1750   
1751   fNumberOfKinks = esd.GetNumberOfKinks();
1752
1753   const AliESDVertex* vtx = esd.GetPrimaryVertex();
1754   
1755   for (Int_t iTrack=0; iTrack<esd.GetNumberOfTracks(); ++iTrack) 
1756   {
1757     AliESDtrack * esdTrack = esd.GetTrack(iTrack);
1758     
1759     Int_t ikink = esdTrack->GetKinkIndex(0);
1760     
1761     if (ikink && fNumberOfKinks) {
1762             // Negative kink index: mother, positive: daughter
1763             
1764             // Search for the second track of the kink
1765             
1766             for (Int_t jTrack = iTrack+1; jTrack<esd.GetNumberOfTracks(); ++jTrack) {
1767         
1768         AliESDtrack * esdTrack1 = esd.GetTrack(jTrack);
1769         
1770         Int_t jkink = esdTrack1->GetKinkIndex(0);
1771         
1772         if ( TMath::Abs(ikink)==TMath::Abs(jkink) ) {
1773           
1774           // The two tracks are from the same kink
1775           
1776           if (fUsedKink[TMath::Abs(ikink)-1]) continue; // skip used kinks
1777           
1778           Int_t imother = -1;
1779           Int_t idaughter = -1;
1780           
1781           if (ikink<0 && jkink>0) {
1782             
1783             imother = iTrack;
1784             idaughter = jTrack;
1785           }
1786           else if (ikink>0 && jkink<0) {
1787             
1788             imother = jTrack;
1789             idaughter = iTrack;
1790           }
1791           else {
1792             //                  cerr << "Error: Wrong combination of kink indexes: "
1793             //                       << ikink << " " << jkink << endl;
1794             continue;
1795           }
1796           
1797           // Add the mother track if it passed primary track selection cuts
1798           
1799           AliAODTrack * mother = NULL;
1800           
1801           UInt_t selectInfo = 0;
1802           if (fTrackFilter) {
1803             selectInfo = fTrackFilter->IsSelected(esd.GetTrack(imother));
1804             if (!selectInfo) continue;
1805           }
1806           
1807           if (!fUsedTrack[imother]) {
1808             
1809             fUsedTrack[imother] = kTRUE;
1810             
1811             AliESDtrack *esdTrackM = esd.GetTrack(imother);
1812             Double_t p[3] = { 0. };
1813             Double_t pos[3] = { 0. };
1814             esdTrackM->GetPxPyPz(p);
1815             esdTrackM->GetXYZ(pos);
1816             esdTrackM->GetCovarianceXYZPxPyPz(covTr);
1817             esdTrackM->GetESDpid(pid);
1818             if(fMChandler)fMChandler->SelectParticle(esdTrackM->GetLabel());
1819             mother = 
1820             new(Tracks()[fNumberOfTracks++]) AliAODTrack(esdTrackM->GetID(),
1821                                                esdTrackM->GetLabel(),
1822                                                p,
1823                                                kTRUE,
1824                                                pos,
1825                                                kFALSE,
1826                                                covTr, 
1827                                                (Short_t)esdTrackM->GetSign(),
1828                                                esdTrackM->GetITSClusterMap(), 
1829                                                pid,
1830                                                fPrimaryVertex,
1831                                                kTRUE, // check if this is right
1832                                                vtx->UsesTrack(esdTrack->GetID()),
1833                                                AliAODTrack::kPrimary,
1834                                                selectInfo);
1835             mother->SetTPCFitMap(esdTrackM->GetTPCFitMap());
1836             mother->SetTPCClusterMap(esdTrackM->GetTPCClusterMap());
1837             mother->SetTPCSharedMap (esdTrackM->GetTPCSharedMap());
1838             mother->SetChi2perNDF(Chi2perNDF(esdTrackM));
1839             mother->SetTPCPointsF(esdTrackM->GetTPCNclsF());
1840             mother->SetTPCNCrossedRows(UShort_t(esdTrackM->GetTPCCrossedRows()));
1841             mother->SetIntegratedLength(esdTrackM->GetIntegratedLength());
1842
1843             fAODTrackRefs->AddAt(mother, imother);
1844             
1845             if (esdTrackM->GetSign() > 0) ++fNumberOfPositiveTracks;
1846             mother->SetFlags(esdTrackM->GetStatus());
1847             mother->ConvertAliPIDtoAODPID();
1848             fPrimaryVertex->AddDaughter(mother);
1849             mother->ConvertAliPIDtoAODPID();
1850             SetAODPID(esdTrackM,mother,detpid);
1851           }
1852           else {
1853             //                  cerr << "Error: event " << esd.GetEventNumberInFile() << " kink " << TMath::Abs(ikink)-1
1854             //                       << " track " << imother << " has already been used!" << endl;
1855           }
1856           
1857           // Add the kink vertex
1858           AliESDkink * kink = esd.GetKink(TMath::Abs(ikink)-1);
1859           
1860           AliAODVertex * vkink = 
1861           new(Vertices()[fNumberOfVertices++]) AliAODVertex(kink->GetPosition(),
1862                                                   NULL,
1863                                                   0.,
1864                                                   mother,
1865                                                   esdTrack->GetID(),  // This is the track ID of the mother's track!
1866                                                   AliAODVertex::kKink);
1867           // Add the daughter track
1868           
1869           AliAODTrack * daughter = NULL;
1870           
1871           if (!fUsedTrack[idaughter]) {
1872             
1873             fUsedTrack[idaughter] = kTRUE;
1874             
1875             AliESDtrack *esdTrackD = esd.GetTrack(idaughter);
1876             Double_t p[3] = { 0. };
1877             Double_t pos[3] = { 0. };
1878
1879             esdTrackD->GetPxPyPz(p);
1880             esdTrackD->GetXYZ(pos);
1881             esdTrackD->GetCovarianceXYZPxPyPz(covTr);
1882             esdTrackD->GetESDpid(pid);
1883             selectInfo = 0;
1884             if (fTrackFilter) selectInfo = fTrackFilter->IsSelected(esdTrackD);
1885             if(fMChandler)fMChandler->SelectParticle(esdTrackD->GetLabel());
1886             daughter = 
1887             new(Tracks()[fNumberOfTracks++]) AliAODTrack(esdTrackD->GetID(),
1888                                                esdTrackD->GetLabel(),
1889                                                p,
1890                                                kTRUE,
1891                                                pos,
1892                                                kFALSE,
1893                                                covTr, 
1894                                                (Short_t)esdTrackD->GetSign(),
1895                                                esdTrackD->GetITSClusterMap(), 
1896                                                pid,
1897                                                vkink,
1898                                                kTRUE, // check if this is right
1899                                                vtx->UsesTrack(esdTrack->GetID()),
1900                                                AliAODTrack::kSecondary,
1901                                                selectInfo);
1902             daughter->SetTPCFitMap(esdTrackD->GetTPCFitMap());
1903             daughter->SetTPCClusterMap(esdTrackD->GetTPCClusterMap());
1904             daughter->SetTPCSharedMap (esdTrackD->GetTPCSharedMap());
1905             daughter->SetTPCPointsF(esdTrackD->GetTPCNclsF());
1906             daughter->SetTPCNCrossedRows(UShort_t(esdTrackD->GetTPCCrossedRows()));
1907             daughter->SetIntegratedLength(esdTrackD->GetIntegratedLength());
1908             fAODTrackRefs->AddAt(daughter, idaughter);
1909             
1910             if (esdTrackD->GetSign() > 0) ++fNumberOfPositiveTracks;
1911             daughter->SetFlags(esdTrackD->GetStatus());
1912             daughter->ConvertAliPIDtoAODPID();
1913             vkink->AddDaughter(daughter);
1914             daughter->ConvertAliPIDtoAODPID();
1915             SetAODPID(esdTrackD,daughter,detpid);
1916           }
1917           else {
1918             //                  cerr << "Error: event " << esd.GetEventNumberInFile() << " kink " << TMath::Abs(ikink)-1
1919             //                       << " track " << idaughter << " has already been used!" << endl;
1920           }
1921         }
1922             }
1923     }      
1924   }
1925 }
1926
1927 //______________________________________________________________________________
1928 void AliAnalysisTaskESDfilter::ConvertPrimaryVertices(const AliESDEvent& esd)
1929 {
1930   AliCodeTimerAuto("",0);
1931   
1932   // Access to the AOD container of vertices
1933   fNumberOfVertices = 0;
1934   
1935   Double_t pos[3] = { 0. };
1936   Double_t covVtx[6] = { 0. };
1937
1938   // Add primary vertex. The primary tracks will be defined
1939   // after the loops on the composite objects (V0, cascades, kinks)
1940   const AliESDVertex *vtx = esd.GetPrimaryVertex();
1941   
1942   vtx->GetXYZ(pos); // position
1943   vtx->GetCovMatrix(covVtx); //covariance matrix
1944   
1945   fPrimaryVertex = new(Vertices()[fNumberOfVertices++])
1946   AliAODVertex(pos, covVtx, vtx->GetChi2toNDF(), NULL, -1, AliAODVertex::kPrimary);
1947   fPrimaryVertex->SetName(vtx->GetName());
1948   fPrimaryVertex->SetTitle(vtx->GetTitle());
1949   fPrimaryVertex->SetBC(vtx->GetBC());
1950   
1951   TString vtitle = vtx->GetTitle();
1952   if (!vtitle.Contains("VertexerTracks")) 
1953     fPrimaryVertex->SetNContributors(vtx->GetNContributors());
1954   
1955   if (fDebug > 0) fPrimaryVertex->Print();  
1956   
1957   // Add SPD "main" vertex 
1958   const AliESDVertex *vtxS = esd.GetPrimaryVertexSPD();
1959   vtxS->GetXYZ(pos); // position
1960   vtxS->GetCovMatrix(covVtx); //covariance matrix
1961   AliAODVertex * mVSPD = new(Vertices()[fNumberOfVertices++])
1962   AliAODVertex(pos, covVtx, vtxS->GetChi2toNDF(), NULL, -1, AliAODVertex::kMainSPD);
1963   mVSPD->SetName(vtxS->GetName());
1964   mVSPD->SetTitle(vtxS->GetTitle());
1965   mVSPD->SetNContributors(vtxS->GetNContributors()); 
1966   
1967   // Add SPD pileup vertices
1968   for(Int_t iV=0; iV<esd.GetNumberOfPileupVerticesSPD(); ++iV)
1969   {
1970     const AliESDVertex *vtxP = esd.GetPileupVertexSPD(iV);
1971     vtxP->GetXYZ(pos); // position
1972     vtxP->GetCovMatrix(covVtx); //covariance matrix
1973     AliAODVertex * pVSPD =  new(Vertices()[fNumberOfVertices++])
1974     AliAODVertex(pos, covVtx, vtxP->GetChi2toNDF(), NULL, -1, AliAODVertex::kPileupSPD);
1975     pVSPD->SetName(vtxP->GetName());
1976     pVSPD->SetTitle(vtxP->GetTitle());
1977     pVSPD->SetNContributors(vtxP->GetNContributors()); 
1978     pVSPD->SetBC(vtxP->GetBC());
1979   }
1980   
1981   // Add TRK pileup vertices
1982   for(Int_t iV=0; iV<esd.GetNumberOfPileupVerticesTracks(); ++iV)
1983   {
1984     const AliESDVertex *vtxP = esd.GetPileupVertexTracks(iV);
1985     vtxP->GetXYZ(pos); // position
1986     vtxP->GetCovMatrix(covVtx); //covariance matrix
1987     AliAODVertex * pVTRK = new(Vertices()[fNumberOfVertices++])
1988     AliAODVertex(pos, covVtx, vtxP->GetChi2toNDF(), NULL, -1, AliAODVertex::kPileupTracks);
1989     pVTRK->SetName(vtxP->GetName());
1990     pVTRK->SetTitle(vtxP->GetTitle());
1991     pVTRK->SetNContributors(vtxP->GetNContributors());
1992     pVTRK->SetBC(vtxP->GetBC());
1993   }
1994
1995   // Add TPC "main" vertex 
1996   const AliESDVertex *vtxT = esd.GetPrimaryVertexTPC();
1997   vtxT->GetXYZ(pos); // position
1998   vtxT->GetCovMatrix(covVtx); //covariance matrix
1999   AliAODVertex * mVTPC = new(Vertices()[fNumberOfVertices++])
2000   AliAODVertex(pos, covVtx, vtxT->GetChi2toNDF(), NULL, -1, AliAODVertex::kMainTPC);
2001   mVTPC->SetName(vtxT->GetName());
2002   mVTPC->SetTitle(vtxT->GetTitle());
2003   mVTPC->SetNContributors(vtxT->GetNContributors()); 
2004
2005
2006 }
2007
2008 //______________________________________________________________________________
2009 void AliAnalysisTaskESDfilter::ConvertVZERO(const AliESDEvent& esd)
2010 {
2011   // Convert VZERO data
2012   AliAODVZERO* vzeroData = AODEvent()->GetVZEROData();
2013   *vzeroData = *(esd.GetVZEROData());
2014 }
2015
2016 //______________________________________________________________________________
2017 void AliAnalysisTaskESDfilter::ConvertTZERO(const AliESDEvent& esd)
2018 {
2019   // Convert TZERO data
2020   const AliESDTZERO* esdTzero = esd.GetESDTZERO(); 
2021   AliAODTZERO* aodTzero = AODEvent()->GetTZEROData();
2022
2023   for (Int_t icase=0; icase<3; icase++){ 
2024     aodTzero->SetT0TOF(    icase, esdTzero->GetT0TOF(icase));
2025     aodTzero->SetT0TOFbest(icase, esdTzero->GetT0TOFbest(icase)); 
2026   }
2027   aodTzero->SetBackgroundFlag(esdTzero->GetBackgroundFlag());
2028   aodTzero->SetPileupFlag(esdTzero->GetPileupFlag());
2029   aodTzero->SetSatelliteFlag(esdTzero->GetSatellite()); 
2030
2031   Float_t rawTime[24];
2032   for(Int_t ipmt=0; ipmt<24; ipmt++)
2033     rawTime[ipmt] = esdTzero->GetTimeFull(ipmt,0);
2034    
2035   Int_t idxOfFirstPmtA = -1,       idxOfFirstPmtC = -1;
2036   Float_t timeOfFirstPmtA = 9999, timeOfFirstPmtC = 9999;
2037   for(int ipmt=0;  ipmt<12; ipmt++){
2038     if( rawTime[ipmt] > -200 && rawTime[ipmt] < timeOfFirstPmtC && rawTime[ipmt]!=0){
2039       timeOfFirstPmtC = rawTime[ipmt];
2040       idxOfFirstPmtC  = ipmt;
2041     }
2042   }
2043   for(int ipmt=12; ipmt<24; ipmt++){
2044     if( rawTime[ipmt] > -200 && rawTime[ipmt] < timeOfFirstPmtA && rawTime[ipmt]!=0 ){
2045       timeOfFirstPmtA = rawTime[ipmt];
2046       idxOfFirstPmtA  = ipmt;
2047     }
2048   }
2049
2050   if(idxOfFirstPmtA != -1 && idxOfFirstPmtC != -1){
2051     //speed of light in cm/ns   TMath::C()*1e-7 
2052     Float_t vertexraw = TMath::C()*1e-7 * (rawTime[idxOfFirstPmtA] - rawTime[idxOfFirstPmtC])/2;
2053     aodTzero->SetT0VertexRaw( vertexraw );
2054   }else{
2055     aodTzero->SetT0VertexRaw(99999);
2056   }
2057
2058   aodTzero->SetT0zVertex(esdTzero->GetT0zVertex());
2059 }
2060
2061
2062 //______________________________________________________________________________
2063 void AliAnalysisTaskESDfilter::ConvertZDC(const AliESDEvent& esd)
2064 {
2065   // Convert ZDC data
2066   AliESDZDC* esdZDC = esd.GetZDCData();
2067   
2068   const Double_t zem1Energy = esdZDC->GetZEM1Energy();
2069   const Double_t zem2Energy = esdZDC->GetZEM2Energy();
2070    
2071   const Double_t *towZNC = esdZDC->GetZNCTowerEnergy();
2072   const Double_t *towZPC = esdZDC->GetZPCTowerEnergy();
2073   const Double_t *towZNA = esdZDC->GetZNATowerEnergy();
2074   const Double_t *towZPA = esdZDC->GetZPATowerEnergy();
2075   const Double_t *towZNCLG = esdZDC->GetZNCTowerEnergyLR();
2076   const Double_t *towZNALG = esdZDC->GetZNATowerEnergyLR();
2077   
2078   AliAODZDC* zdcAOD = AODEvent()->GetZDCData();
2079
2080   zdcAOD->SetZEM1Energy(zem1Energy);
2081   zdcAOD->SetZEM2Energy(zem2Energy);
2082   zdcAOD->SetZNCTowers(towZNC, towZNCLG);
2083   zdcAOD->SetZNATowers(towZNA, towZNALG);
2084   zdcAOD->SetZPCTowers(towZPC);
2085   zdcAOD->SetZPATowers(towZPA);
2086   
2087   zdcAOD->SetZDCParticipants(esdZDC->GetZDCParticipants(), esdZDC->GetZDCPartSideA(), esdZDC->GetZDCPartSideC());
2088   zdcAOD->SetZDCImpactParameter(esdZDC->GetImpactParameter(), esdZDC->GetImpactParamSideA(), 
2089         esdZDC->GetImpactParamSideC());
2090   zdcAOD->SetZDCTDCSum(esdZDC->GetZNTDCSum(0)); 
2091   zdcAOD->SetZDCTDCDiff(esdZDC->GetZNTDCDiff(0));       
2092   if(esdZDC->IsZNChit()) zdcAOD->SetZNCTDC(esdZDC->GetZDCTDCCorrected(10,0));
2093   if(esdZDC->IsZNAhit()) zdcAOD->SetZNATDC(esdZDC->GetZDCTDCCorrected(12,0));
2094 }
2095
2096 //_______________________________________________________________________________________________________________________________________
2097 Int_t AliAnalysisTaskESDfilter::ConvertHMPID(const AliESDEvent& esd) // clm
2098 {
2099   //
2100   // Convtert ESD HMPID info to AOD and return the number of good tracks with HMPID signal.
2101   // We need to return an int since there is no signal counter in the ESD.
2102   //
2103   
2104   AliCodeTimerAuto("",0);
2105   
2106   Int_t cntHmpidGoodTracks = 0;
2107   
2108   Float_t  xMip = 0;
2109   Float_t  yMip = 0;
2110   Int_t    qMip = 0;
2111   Int_t    nphMip = 0;
2112   
2113   Float_t xTrk = 0;
2114   Float_t yTrk = 0;
2115   Float_t thetaTrk = 0;
2116   Float_t phiTrk = 0;
2117   
2118   Double_t hmpPid[5]={0};
2119   Double_t hmpMom[3]={0};
2120   
2121   TClonesArray &hmpidRings = *(AODEvent()->GetHMPIDrings());
2122   
2123  for (Int_t iTrack=0; iTrack<esd.GetNumberOfTracks(); ++iTrack) 
2124   {
2125     if(! esd.GetTrack(iTrack) ) continue;
2126       
2127     if(esd.GetTrack(iTrack)->GetHMPIDsignal() > -20 ) {                  // 
2128        
2129       (esd.GetTrack(iTrack))->GetHMPIDmip(xMip, yMip, qMip, nphMip);     // Get MIP properties
2130       (esd.GetTrack(iTrack))->GetHMPIDtrk(xTrk,yTrk,thetaTrk,phiTrk);
2131       (esd.GetTrack(iTrack))->GetHMPIDpid(hmpPid);
2132       if((esd.GetTrack(iTrack))->GetOuterHmpParam()) (esd.GetTrack(iTrack))->GetOuterHmpPxPyPz(hmpMom);
2133       
2134      if(esd.GetTrack(iTrack)->GetHMPIDsignal() == 0 && thetaTrk == 0 && qMip == 0 && nphMip ==0 ) continue;  //
2135       
2136      new(hmpidRings[cntHmpidGoodTracks++]) AliAODHMPIDrings(
2137                                                                  (esd.GetTrack(iTrack))->GetID(),             // Unique track id to attach the ring to
2138                                                                  1000000*nphMip+qMip,                         // MIP charge and number of photons
2139                                                                  (esd.GetTrack(iTrack))->GetHMPIDcluIdx(),    // 1000000*chamber id + cluster idx of the assigned MIP cluster  
2140                                                                  thetaTrk,                                    // track inclination angle theta
2141                                                                  phiTrk,                                      // track inclination angle phi
2142                                                                  (esd.GetTrack(iTrack))->GetHMPIDsignal(),    // Cherenkov angle
2143                                                                  (esd.GetTrack(iTrack))->GetHMPIDoccupancy(), // Occupancy claculated for the given chamber 
2144                                                                  (esd.GetTrack(iTrack))->GetHMPIDchi2(),      // Ring resolution squared
2145                                                                   xTrk,                                       // Track x coordinate (LORS)
2146                                                                   yTrk,                                       // Track y coordinate (LORS)
2147                                                                   xMip,                                       // MIP x coordinate (LORS)
2148                                                                   yMip,                                       // MIP y coordinate (LORS)
2149                                                                   hmpPid,                                     // PID probablities from ESD, remove later once it is in CombinedPid
2150                                                                   hmpMom                                      // Track momentum in HMPID at ring reconstruction
2151                                                                );  
2152      
2153       //  Printf(Form("+++++++++ yes/no: %d  %lf %lf %lf %lf %lf %lf ",(esd.GetTrack(iTrack))->IsHMPID(),thetaTrk, (esd.GetTrack(iTrack))->GetHMPIDchi2(),xTrk, yTrk , xMip, yMip));
2154      
2155                                                                 
2156    }// HMPID signal > -20
2157   }//___esd track loop
2158   
2159   return cntHmpidGoodTracks;
2160 }
2161
2162 //______________________________________________________________________________
2163 void AliAnalysisTaskESDfilter::ConvertESDtoAOD() 
2164 {
2165   // ESD Filter analysis task executed for each event
2166   
2167   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
2168   
2169   if(!esd)return;
2170
2171   AliCodeTimerAuto("",0);
2172   
2173   fOldESDformat = ( esd->GetAliESDOld() != 0x0 );
2174  
2175       // Reconstruct cascades and V0 here
2176   if (fIsV0CascadeRecoEnabled) {
2177     esd->ResetCascades();
2178     esd->ResetV0s();
2179
2180     AliV0vertexer lV0vtxer;
2181     AliCascadeVertexer lCascVtxer;
2182
2183     lV0vtxer.SetCuts(fV0Cuts);
2184     lCascVtxer.SetCuts(fCascadeCuts);
2185
2186
2187     lV0vtxer.Tracks2V0vertices(esd);
2188     lCascVtxer.V0sTracks2CascadeVertices(esd);
2189   }
2190
2191  
2192   fNumberOfTracks = 0;
2193   fNumberOfPositiveTracks = 0;
2194   fNumberOfV0s = 0;
2195   fNumberOfVertices = 0;
2196   fNumberOfCascades = 0;
2197   fNumberOfKinks = 0;
2198     
2199   AliAODHeader* header = ConvertHeader(*esd);
2200
2201   if ( fIsVZEROEnabled ) ConvertVZERO(*esd);
2202   if ( fIsTZEROEnabled ) ConvertTZERO(*esd);
2203   
2204   // Fetch Stack for debuggging if available 
2205   fMChandler=0x0;
2206   if(MCEvent())
2207   {
2208     fMChandler = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler()); 
2209   }
2210   
2211   // loop over events and fill them
2212   // Multiplicity information needed by the header (to be revised!)
2213   Int_t nTracks    = esd->GetNumberOfTracks();
2214   for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) esd->GetTrack(iTrack)->SetESDEvent(esd);
2215
2216   // Update the header
2217
2218   Int_t nV0s      = esd->GetNumberOfV0s();
2219   Int_t nCascades = esd->GetNumberOfCascades();
2220   Int_t nKinks    = esd->GetNumberOfKinks();
2221   Int_t nVertices = nV0s + nCascades /*V0 wihtin cascade already counted*/+ nKinks + 1 /* = prim. vtx*/;
2222   Int_t nPileSPDVertices=1+esd->GetNumberOfPileupVerticesSPD(); // also SPD main vertex
2223   Int_t nPileTrkVertices=esd->GetNumberOfPileupVerticesTracks();
2224   nVertices+=nPileSPDVertices;
2225   nVertices+=nPileTrkVertices;
2226   Int_t nJets     = 0;
2227   Int_t nCaloClus = esd->GetNumberOfCaloClusters();
2228   Int_t nFmdClus  = 0;
2229   Int_t nPmdClus  = esd->GetNumberOfPmdTracks();
2230   Int_t nHmpidRings = 0;  
2231     
2232   AliDebug(1,Form("   NV0=%d  NCASCADES=%d  NKINKS=%d", nV0s, nCascades, nKinks));
2233        
2234   AODEvent()->ResetStd(nTracks, nVertices, nV0s, nCascades, nJets, nCaloClus, nFmdClus, nPmdClus,nHmpidRings);
2235
2236   if (nV0s > 0) 
2237   {
2238     // RefArray to store a mapping between esd V0 number and newly created AOD-Vertex V0
2239     fAODV0VtxRefs = new TRefArray(nV0s);
2240     // RefArray to store the mapping between esd V0 number and newly created AOD-V0
2241     fAODV0Refs = new TRefArray(nV0s); 
2242     // Array to take into account the V0s already added to the AOD (V0 within cascades)
2243     fUsedV0 = new Bool_t[nV0s];
2244     for (Int_t iV0=0; iV0<nV0s; ++iV0) fUsedV0[iV0]=kFALSE;
2245   }
2246   
2247   if (nTracks>0) 
2248   {
2249     // RefArray to store the mapping between esd track number and newly created AOD-Track
2250     
2251     fAODTrackRefs = new TRefArray(nTracks);
2252
2253     // Array to take into account the tracks already added to the AOD    
2254     fUsedTrack = new Bool_t[nTracks];
2255     for (Int_t iTrack=0; iTrack<nTracks; ++iTrack) fUsedTrack[iTrack]=kFALSE;
2256   }
2257   
2258   // Array to take into account the kinks already added to the AOD
2259   if (nKinks>0) 
2260   {
2261     fUsedKink = new Bool_t[nKinks];
2262     for (Int_t iKink=0; iKink<nKinks; ++iKink) fUsedKink[iKink]=kFALSE;
2263   }
2264     
2265   ConvertPrimaryVertices(*esd);
2266
2267   //setting best TOF PID
2268   AliESDInputHandler* esdH = dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
2269   if (esdH)
2270       fESDpid = esdH->GetESDpid();
2271
2272   if (fIsPidOwner && fESDpid){
2273     delete fESDpid;
2274     fESDpid = 0;
2275   }
2276   if(!fESDpid)
2277   { //in case of no Tender attached 
2278     fESDpid = new AliESDpid;
2279     fIsPidOwner = kTRUE;
2280   }
2281   
2282   if(!esd->GetTOFHeader())
2283   { //protection in case the pass2 LHC10b,c,d have been processed without tender. 
2284     Float_t t0spread[10];
2285     Float_t intrinsicTOFres=100; //ps ok for LHC10b,c,d pass2!! 
2286     for (Int_t i=0; i<10; i++) t0spread[i] = (TMath::Sqrt(esd->GetSigma2DiamondZ()))/0.03; //0.03 to convert from cm to ps
2287     fESDpid->GetTOFResponse().SetT0resolution(t0spread);
2288     fESDpid->GetTOFResponse().SetTimeResolution(intrinsicTOFres);
2289       //    fESDpid->SetTOFResponse(esd, (AliESDpid::EStartTimeType_t)fTimeZeroType);    
2290     AliTOFHeader tmpTOFHeader(0,t0spread[0],0,NULL,NULL,NULL,intrinsicTOFres,t0spread[0]);   
2291     AODEvent()->SetTOFHeader(&tmpTOFHeader);         // write dummy TOF header in AOD
2292   } else {
2293     AODEvent()->SetTOFHeader(esd->GetTOFHeader());    // write TOF header in AOD
2294   }
2295   
2296   //  if(esd->GetTOFHeader() && fIsPidOwner) fESDpid->SetTOFResponse(esd, (AliESDpid::EStartTimeType_t)fTimeZeroType); //in case of AOD production strating form LHC10e without Tender. 
2297   
2298   if ( fAreCascadesEnabled ) ConvertCascades(*esd);
2299
2300   if ( fAreV0sEnabled ) ConvertV0s(*esd);
2301   
2302   if ( fAreKinksEnabled ) ConvertKinks(*esd);
2303   
2304   if ( fAreTracksEnabled ) ConvertTracks(*esd);
2305   
2306   // Update number of AOD tracks in header at the end of track loop (M.G.)
2307   header->SetRefMultiplicity(fNumberOfTracks);
2308   header->SetRefMultiplicityPos(fNumberOfPositiveTracks);
2309   header->SetRefMultiplicityNeg(fNumberOfTracks - fNumberOfPositiveTracks);
2310
2311   if ( fTPCConstrainedFilterMask ) ConvertTPCOnlyTracks(*esd);
2312   if( fGlobalConstrainedFilterMask) ConvertGlobalConstrainedTracks(*esd);  
2313
2314   if ( fArePmdClustersEnabled ) ConvertPmdClusters(*esd);
2315   
2316   if ( fAreCaloClustersEnabled ) ConvertCaloClusters(*esd);
2317   
2318   if ( fAreEMCALCellsEnabled )ConvertEMCALCells(*esd);
2319   
2320   if ( fArePHOSCellsEnabled )ConvertPHOSCells(*esd);
2321         
2322         if ( fAreEMCALTriggerEnabled )ConvertCaloTrigger(TString("EMCAL"), *esd);
2323
2324         if ( fArePHOSTriggerEnabled )ConvertCaloTrigger(TString("PHOS"), *esd);
2325   
2326   if ( fAreTrackletsEnabled ) ConvertTracklets(*esd);
2327   if ( fIsZDCEnabled ) ConvertZDC(*esd);
2328   
2329   if(fIsHMPIDEnabled) nHmpidRings = ConvertHMPID(*esd); 
2330
2331   delete fAODTrackRefs; fAODTrackRefs=0x0;
2332   delete fAODV0VtxRefs; fAODV0VtxRefs=0x0;
2333   delete fAODV0Refs; fAODV0Refs=0x0;
2334   
2335   delete[] fUsedTrack; fUsedTrack=0x0;
2336   delete[] fUsedV0; fUsedV0=0x0;
2337   delete[] fUsedKink; fUsedKink=0x0;
2338
2339   if ( fIsPidOwner){
2340     delete fESDpid;
2341     fESDpid = 0x0;
2342   }
2343
2344
2345 }
2346
2347
2348 //______________________________________________________________________________
2349 void AliAnalysisTaskESDfilter::SetAODPID(AliESDtrack *esdtrack, AliAODTrack *aodtrack, AliAODPid *detpid)
2350 {
2351   //
2352   // Setter for the raw PID detector signals
2353   //
2354
2355   // Save PID object for candidate electrons
2356     Bool_t pidSave = kFALSE;
2357     if (fTrackFilter) {
2358         Bool_t selectInfo = fTrackFilter->IsSelected((char*) "Electrons");
2359         if (selectInfo)  pidSave = kTRUE;
2360     }
2361
2362
2363     // Tracks passing pt cut 
2364     if(esdtrack->Pt()>fHighPthreshold) {
2365         pidSave = kTRUE;
2366     } else {
2367         if(fPtshape){
2368             if(esdtrack->Pt()> fPtshape->GetXmin()){
2369                 Double_t y = fPtshape->Eval(esdtrack->Pt())/fPtshape->Eval(fHighPthreshold);
2370                 if(gRandom->Rndm(0)<1./y){
2371                     pidSave = kTRUE;
2372                 }//end rndm
2373             }//end if p < pmin
2374         }//end if p function
2375     }// end else
2376
2377     if (pidSave) {
2378       if(!aodtrack->GetDetPid()){// prevent memory leak when calling SetAODPID twice for the same track
2379         detpid = new AliAODPid();
2380         SetDetectorRawSignals(detpid,esdtrack);
2381         aodtrack->SetDetPID(detpid);
2382       }
2383     }
2384 }
2385
2386 //______________________________________________________________________________
2387 void AliAnalysisTaskESDfilter::SetDetectorRawSignals(AliAODPid *aodpid, AliESDtrack *track)
2388 {
2389 //
2390 //assignment of the detector signals (AliXXXesdPID inspired)
2391 //
2392  if(!track){
2393  AliInfo("no ESD track found. .....exiting");
2394  return;
2395  }
2396  // TPC momentum
2397  const AliExternalTrackParam *in=track->GetInnerParam();
2398  if (in) {
2399    aodpid->SetTPCmomentum(in->GetP());
2400  }else{
2401    aodpid->SetTPCmomentum(-1.);
2402  }
2403
2404
2405  aodpid->SetITSsignal(track->GetITSsignal());
2406  Double_t itsdedx[4]; // dE/dx samples for individual ITS layers
2407  track->GetITSdEdxSamples(itsdedx);
2408  aodpid->SetITSdEdxSamples(itsdedx);
2409
2410  aodpid->SetTPCsignal(track->GetTPCsignal());
2411  aodpid->SetTPCsignalN(track->GetTPCsignalN());
2412  if(track->GetTPCdEdxInfo()) aodpid->SetTPCdEdxInfo(track->GetTPCdEdxInfo());
2413
2414  //n TRD planes = 6
2415  Int_t nslices = track->GetNumberOfTRDslices()*6;
2416  TArrayD trdslices(nslices);
2417  for(Int_t iSl =0; iSl < track->GetNumberOfTRDslices(); iSl++) {
2418    for(Int_t iPl =0; iPl<6; iPl++) trdslices[iPl*track->GetNumberOfTRDslices()+iSl] = track->GetTRDslice(iPl,iSl);
2419  }
2420  
2421 //TRD momentum
2422  for(Int_t iPl=0;iPl<6;iPl++){
2423    Double_t trdmom=track->GetTRDmomentum(iPl);
2424    aodpid->SetTRDmomentum(iPl,trdmom);
2425  }
2426
2427  aodpid->SetTRDslices(track->GetNumberOfTRDslices()*6,trdslices.GetArray());
2428  aodpid->SetTRDsignal(track->GetTRDsignal());
2429
2430  //TRD clusters and tracklets
2431  aodpid->SetTRDncls(track->GetTRDncls());
2432  aodpid->SetTRDntrackletsPID(track->GetTRDntrackletsPID());
2433  
2434  aodpid->SetTRDChi2(track->GetTRDchi2());
2435
2436  //TOF PID  
2437  Double_t times[AliPID::kSPECIES]; track->GetIntegratedTimes(times);
2438  aodpid->SetIntegratedTimes(times);
2439
2440    //  Float_t tzeroTrack = fESDpid->GetTOFResponse().GetStartTime(track->P());
2441    //  aodpid->SetTOFsignal(track->GetTOFsignal()-tzeroTrack);
2442    aodpid->SetTOFsignal(track->GetTOFsignal());
2443   
2444   Double_t tofRes[5];
2445   for (Int_t iMass=0; iMass<5; iMass++){
2446     //    tofRes[iMass]=(Double_t)fESDpid->GetTOFResponse().GetExpectedSigma(track->P(), times[iMass], AliPID::ParticleMass(iMass));
2447     tofRes[iMass]=0; //backward compatibility
2448   }
2449   aodpid->SetTOFpidResolution(tofRes);
2450
2451 //  aodpid->SetHMPIDsignal(0); // set to zero for compression but it will be removed later
2452
2453 }
2454
2455 Double_t  AliAnalysisTaskESDfilter::Chi2perNDF(AliESDtrack* track)
2456 {
2457     // Calculate chi2 per ndf for track
2458     Int_t  nClustersTPC = track->GetTPCNcls();
2459
2460     if ( nClustersTPC > 5) {
2461        return (track->GetTPCchi2()/Float_t(nClustersTPC - 5));
2462     } else {
2463        return (-1.);
2464     }
2465  }
2466
2467
2468 //______________________________________________________________________________
2469 void AliAnalysisTaskESDfilter::Terminate(Option_t */*option*/)
2470 {
2471 // Terminate analysis
2472 //
2473     if (fDebug > 1) printf("AnalysisESDfilter: Terminate() \n");
2474 }
2475
2476 //______________________________________________________________________________
2477 void  AliAnalysisTaskESDfilter::PrintMCInfo(AliStack *pStack,Int_t label){
2478 // Print MC info
2479   if(!pStack)return;
2480   label = TMath::Abs(label);
2481   TParticle *part = pStack->Particle(label);
2482   Printf("########################");
2483   Printf("%s:%d %d UniqueID %d PDG %d P %3.3f",(char*)__FILE__,__LINE__,label,part->GetUniqueID(),part->GetPdgCode(),part->P());
2484   part->Print();
2485   TParticle* mother = part;
2486   Int_t imo = part->GetFirstMother();
2487   Int_t nprim = pStack->GetNprimary();
2488   //  while((imo >= nprim) && (mother->GetUniqueID() == 4)) {
2489   while((imo >= nprim)) {
2490     mother =  pStack->Particle(imo);
2491     Printf("Mother %s:%d Label %d UniqueID %d PDG %d P %3.3f",(char*)__FILE__,__LINE__,imo,mother->GetUniqueID(),mother->GetPdgCode(),mother->P());
2492     mother->Print();
2493     imo =  mother->GetFirstMother();
2494   }
2495   Printf("########################");
2496 }
2497
2498 //______________________________________________________
2499