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