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