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