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