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