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