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