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