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