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