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