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