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