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