including possibility of selecting negatively charged kaons and pions
[u/mrichter/AliRoot.git] / PWGCF / FEMTOSCOPY / AliFemto / AliFemtoEventReaderESDChainKine.cxx
CommitLineData
76ce4b5b 1////////////////////////////////////////////////////////////////////////////////
2// //
3// AliFemtoEventReaderESDChainKine - the reader class for the Alice ESD and //
4// the model Kinematics information tailored for the Task framework and the //
5// Reads in AliESDfriend to create shared hit/quality information //
6// Authors: Adam Kisiel kisiel@mps.ohio-state.edu //
7// //
8////////////////////////////////////////////////////////////////////////////////
9#include "AliFemtoEventReaderESDChainKine.h"
10
11#include "TFile.h"
12#include "TTree.h"
13#include "TList.h"
14#include "AliESDEvent.h"
15#include "AliESDtrack.h"
16#include "AliESDVertex.h"
17
18#include "AliMultiplicity.h"
19#include "AliCentrality.h"
20#include "AliEventplane.h"
21
22#include "AliFmPhysicalHelixD.h"
23#include "AliFmThreeVectorF.h"
24
25#include "SystemOfUnits.h"
26
27#include "AliFemtoEvent.h"
28
29#include "TParticle.h"
30#include "AliFemtoModelHiddenInfo.h"
31#include "AliFemtoModelGlobalHiddenInfo.h"
32#include "AliGenHijingEventHeader.h"
33#include "AliGenCocktailEventHeader.h"
34
35#include "AliVertexerTracks.h"
36
37#include "AliPID.h"
38
39ClassImp(AliFemtoEventReaderESDChainKine)
40
41#if !(ST_NO_NAMESPACES)
a68b227d 42using namespace units;
76ce4b5b 43#endif
44
45using namespace std;
46//____________________________
47AliFemtoEventReaderESDChainKine::AliFemtoEventReaderESDChainKine():
48 fFileName(" "),
49 fConstrained(true),
50 fReadInner(false),
51 fUseTPCOnly(false),
52 fNumberofEvent(0),
53 fCurEvent(0),
54 fCurFile(0),
55 fEvent(0x0),
56 fStack(0x0),
57 fGenHeader(0x0),
58 fTrackType(kGlobal),
1938eadf 59 fEstEventMult(kReferenceITSTPC),
76ce4b5b 60 fRotateToEventPlane(0),
61 fESDpid(0),
ce7b3d98 62 fIsPidOwner(0),
63 fMagFieldSign(0),
48178c0e 64 fReadV0(0),
1dd4ea76 65 isKaonAnalysis(kFALSE),
d2eafac6 66 isProtonAnalysis(kFALSE),
89033e57 67 isPionAnalysis(kFALSE),
1dd4ea76 68 fOnlyPrimaries(kFALSE)
76ce4b5b 69
70{
48178c0e 71 //constructor with 0 parameters , look at default settings
76ce4b5b 72}
73
74//__________________
75AliFemtoEventReaderESDChainKine::AliFemtoEventReaderESDChainKine(const AliFemtoEventReaderESDChainKine& aReader):
76 AliFemtoEventReader(aReader),
77 fFileName(" "),
78 fConstrained(true),
79 fReadInner(false),
80 fUseTPCOnly(false),
81 fNumberofEvent(0),
82 fCurEvent(0),
83 fCurFile(0),
84 fEvent(0x0),
85 fStack(0x0),
86 fGenHeader(0x0),
87 fTrackType(kGlobal),
1938eadf 88 fEstEventMult(kReferenceITSTPC),
76ce4b5b 89 fRotateToEventPlane(0),
90 fESDpid(0),
ce7b3d98 91 fIsPidOwner(0),
92 fMagFieldSign(0),
48178c0e 93 fReadV0(0),
1dd4ea76 94 isKaonAnalysis(kFALSE),
d2eafac6 95 isProtonAnalysis(kFALSE),
89033e57 96 isPionAnalysis(kFALSE),
1dd4ea76 97 fOnlyPrimaries(kFALSE)
98
48178c0e 99
76ce4b5b 100{
101 // Copy constructor
102 fConstrained = aReader.fConstrained;
103 fReadInner = aReader.fReadInner;
104 fUseTPCOnly = aReader.fUseTPCOnly;
105 fNumberofEvent = aReader.fNumberofEvent;
106 fCurEvent = aReader.fCurEvent;
107 fCurFile = aReader.fCurFile;
108 fEvent = new AliESDEvent();
109 fStack = aReader.fStack;
110 fTrackType = aReader.fTrackType;
111 fEstEventMult = aReader.fEstEventMult;
112 fRotateToEventPlane = aReader.fRotateToEventPlane;
113 fESDpid = aReader.fESDpid;
114 fIsPidOwner = aReader.fIsPidOwner;
ce7b3d98 115 fMagFieldSign = aReader.fMagFieldSign;
116 fReadV0 = aReader.fReadV0;
48178c0e 117 isKaonAnalysis = aReader.isKaonAnalysis;
d2eafac6 118 isProtonAnalysis = aReader.isProtonAnalysis;
89033e57 119 isPionAnalysis = aReader.isPionAnalysis;
1dd4ea76 120 fOnlyPrimaries = aReader.fOnlyPrimaries;
48178c0e 121
76ce4b5b 122}
123//__________________
124AliFemtoEventReaderESDChainKine::~AliFemtoEventReaderESDChainKine()
125{
126 //Destructor
127 delete fEvent;
128}
129
130//__________________
131AliFemtoEventReaderESDChainKine& AliFemtoEventReaderESDChainKine::operator=(const AliFemtoEventReaderESDChainKine& aReader)
132{
133 // Assignment operator
134 if (this == &aReader)
135 return *this;
136
137 fConstrained = aReader.fConstrained;
138 fUseTPCOnly = aReader.fUseTPCOnly;
139 fNumberofEvent = aReader.fNumberofEvent;
140 fCurEvent = aReader.fCurEvent;
141 fCurFile = aReader.fCurFile;
142 if (fEvent) delete fEvent;
143 fEvent = new AliESDEvent();
144 fStack = aReader.fStack;
145 fGenHeader = aReader.fGenHeader;
146 fRotateToEventPlane = aReader.fRotateToEventPlane;
147 fESDpid = aReader.fESDpid;
148 fIsPidOwner = aReader.fIsPidOwner;
ce7b3d98 149 fMagFieldSign = aReader.fMagFieldSign;
150 fReadV0 = aReader.fReadV0;
48178c0e 151 isKaonAnalysis = aReader.isKaonAnalysis;
d2eafac6 152 isProtonAnalysis = aReader.isProtonAnalysis;
89033e57 153 isPionAnalysis = aReader.isPionAnalysis;
1dd4ea76 154 fOnlyPrimaries = aReader.fOnlyPrimaries;
76ce4b5b 155
156 return *this;
157}
158//__________________
159// Simple report
160AliFemtoString AliFemtoEventReaderESDChainKine::Report()
161{
162 AliFemtoString temp = "\n This is the AliFemtoEventReaderESDChainKine\n";
163 return temp;
164}
165
166//__________________
167void AliFemtoEventReaderESDChainKine::SetConstrained(const bool constrained)
168{
169 // Select whether to read constrained or not constrained momentum
170 fConstrained=constrained;
171}
172//__________________
173bool AliFemtoEventReaderESDChainKine::GetConstrained() const
174{
175 // Check whether we read constrained or not constrained momentum
176 return fConstrained;
177}
178//__________________
179void AliFemtoEventReaderESDChainKine::SetReadTPCInner(const bool readinner)
180{
181 fReadInner=readinner;
182}
183
184bool AliFemtoEventReaderESDChainKine::GetReadTPCInner() const
185{
186 return fReadInner;
187}
188
189//__________________
190void AliFemtoEventReaderESDChainKine::SetUseTPCOnly(const bool usetpconly)
191{
192 fUseTPCOnly=usetpconly;
193}
194
195bool AliFemtoEventReaderESDChainKine::GetUseTPCOnly() const
196{
197 return fUseTPCOnly;
198}
199void AliFemtoEventReaderESDChainKine::SetUseMultiplicity(EstEventMult aType)
200{
201 fEstEventMult = aType;
202}
203//__________________
204AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent()
205{
206 // Get the event, read all the relevant information
207 // and fill the AliFemtoEvent class
208 // Returns a valid AliFemtoEvent
209 AliFemtoEvent *hbtEvent = 0;
210 string tFriendFileName;
211
212 // Get the friend information
213 cout << "AliFemtoEventReaderESDChainKine::Starting to read event: "<<fCurEvent<<endl;
214 // fEvent->SetESDfriend(fEventFriend);
48178c0e 215
76ce4b5b 216 hbtEvent = new AliFemtoEvent;
217 //setting basic things
218 // hbtEvent->SetEventNumber(fEvent->GetEventNumber());
219 hbtEvent->SetRunNumber(fEvent->GetRunNumber());
220 //hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
221 hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
222 hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
223 hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
224 hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
225 hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
226 hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
227 hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
228 hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
229 //hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
230
231 if ((fEvent->IsTriggerClassFired("CINT1WU-B-NOPF-ALL")) ||
232 (fEvent->IsTriggerClassFired("CINT1B-ABCE-NOPF-ALL")) ||
233 (fEvent->IsTriggerClassFired("CINT1-B-NOPF-ALLNOTRD")) ||
234 (fEvent->IsTriggerClassFired("CINT1-B-NOPF-FASTNOTRD")))
235 hbtEvent->SetTriggerCluster(1);
236 else if ((fEvent->IsTriggerClassFired("CSH1WU-B-NOPF-ALL")) ||
a68b227d 237 (fEvent->IsTriggerClassFired("CSH1-B-NOPF-ALLNOTRD")))
76ce4b5b 238 hbtEvent->SetTriggerCluster(2);
48178c0e 239 else
76ce4b5b 240 hbtEvent->SetTriggerCluster(0);
241
242 //Vertex
243 double fV1[3];
244 double fVCov[6];
245 if (fUseTPCOnly) {
246 fEvent->GetPrimaryVertexTPC()->GetXYZ(fV1);
247 fEvent->GetPrimaryVertexTPC()->GetCovMatrix(fVCov);
248 if (!fEvent->GetPrimaryVertexTPC()->GetStatus())
249 fVCov[4] = -1001.0;
250 }
251 else {
252 fEvent->GetPrimaryVertex()->GetXYZ(fV1);
253 fEvent->GetPrimaryVertex()->GetCovMatrix(fVCov);
254 if (!fEvent->GetPrimaryVertex()->GetStatus())
255 fVCov[4] = -1001.0;
256 }
257
258 AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
259 hbtEvent->SetPrimVertPos(vertex);
260 hbtEvent->SetPrimVertCov(fVCov);
261
262 Double_t tReactionPlane = 0;
263
48178c0e 264 AliGenHijingEventHeader *hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader);
76ce4b5b 265 if (!hdh) {
266 AliGenCocktailEventHeader *cdh = dynamic_cast<AliGenCocktailEventHeader *> (fGenHeader);
267 if (cdh) {
268 TList *tGenHeaders = cdh->GetHeaders();
269 for (int ihead = 0; ihead<tGenHeaders->GetEntries(); ihead++) {
a68b227d 270 hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader);
271 if (hdh) break;
76ce4b5b 272 }
273 }
274 }
275
276 if (hdh)
a68b227d 277 {
278 tReactionPlane = hdh->ReactionPlaneAngle();
279 cout << "Got reaction plane " << tReactionPlane << endl;
280 }
76ce4b5b 281
282 hbtEvent->SetReactionPlaneAngle(tReactionPlane);
283
284 Int_t spdetaonecount = 0;
48178c0e 285
286 for (int iter=0; iter<fEvent->GetMultiplicity()->GetNumberOfTracklets(); iter++)
76ce4b5b 287 if (fabs(fEvent->GetMultiplicity()->GetEta(iter)) < 1.0)
288 spdetaonecount++;
289
290 // hbtEvent->SetSPDMult(fEvent->GetMultiplicity()->GetNumberOfTracklets());
291 hbtEvent->SetSPDMult(spdetaonecount);
292
293 //starting to reading tracks
294 int nofTracks=0; //number of reconstructed tracks in event
295 nofTracks=fEvent->GetNumberOfTracks();
296 int realnofTracks=0;//number of track which we use ina analysis
297
298 int tNormMult = 0;
299 int tNormMultPos = 0;
300 int tNormMultNeg = 0;
301
302 Float_t tTotalPt = 0.0;
303
304 Float_t b[2];
305 Float_t bCov[3];
306
1938eadf 307 Int_t tTracklet=0, tITSTPC=0;
48178c0e 308
1938eadf 309 //fEvent->EstimateMultiplicity(tTracklet, tITSTPC, tITSPure, 1.2);
48178c0e 310
76ce4b5b 311 hbtEvent->SetMultiplicityEstimateITSTPC(tITSTPC);
312 hbtEvent->SetMultiplicityEstimateTracklets(tTracklet);
313 // hbtEvent->SetMultiplicityEstimateITSPure(tITSPure);
314 hbtEvent->SetMultiplicityEstimateITSPure(fEvent->GetMultiplicity()->GetNumberOfITSClusters(1));
315
316 Int_t *motherids;
317 if (fStack) {
318 motherids = new Int_t[fStack->GetNtrack()];
319 for (int ip=0; ip<fStack->GetNtrack(); ip++) motherids[ip] = 0;
320
321 // Read in mother ids
322 TParticle *motherpart;
323 for (int ip=0; ip<fStack->GetNtrack(); ip++) {
324 motherpart = fStack->Particle(ip);
325 if (motherpart->GetDaughter(0) > 0)
a68b227d 326 motherids[motherpart->GetDaughter(0)] = ip;
76ce4b5b 327 if (motherpart->GetDaughter(1) > 0)
a68b227d 328 motherids[motherpart->GetDaughter(1)] = ip;
48178c0e 329
76ce4b5b 330// if (motherpart->GetPdgCode() == 211) {
48178c0e 331// cout << "Mother " << ip << " has daughters "
332// << motherpart->GetDaughter(0) << " "
333// << motherpart->GetDaughter(1) << " "
334// << motherpart->Vx() << " "
335// << motherpart->Vy() << " "
336// << motherpart->Vz() << " "
76ce4b5b 337// << endl;
48178c0e 338
76ce4b5b 339// }
340 }
341 }
342 else {
343 printf ("No Stack ???\n");
344 delete hbtEvent;
345 return 0;
346 }
347
348 for (int i=0;i<nofTracks;i++)
a68b227d 349 {
350 // cout << "Reading track " << i << endl;
351 bool tGoodMomentum=true; //flaga to chcek if we can read momentum of this track
352
353
354 const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
355 // const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
356
357
358 if ((esdtrack->GetStatus() & AliESDtrack::kTPCrefit) &&
359 (esdtrack->GetStatus() & AliESDtrack::kITSrefit)) {
360 if (esdtrack->GetTPCNcls() > 70)
361 if (esdtrack->GetTPCchi2()/esdtrack->GetTPCNcls() < 4.0) {
362 if (TMath::Abs(esdtrack->Eta()) < 1.2) {
363 esdtrack->GetImpactParameters(b,bCov);
364 if ((b[0]<0.2) && (b[1] < 0.25)) {
365 tNormMult++;
366 tTotalPt += esdtrack->Pt();
367 }
368 }
369 }
370 }
371 else if (esdtrack->GetStatus() & AliESDtrack::kTPCrefit) {
372 if (esdtrack->GetTPCNcls() > 100)
373 if (esdtrack->GetTPCchi2()/esdtrack->GetTPCNcls() < 4.0) {
374 if (TMath::Abs(esdtrack->Eta()) < 1.2) {
375 esdtrack->GetImpactParameters(b,bCov);
376 if ((b[0]<2.4) && (b[1] < 3.2)) {
377 tNormMult++;
378 tTotalPt += esdtrack->Pt();
379 }
380 }
381 }
382 }
48178c0e 383
a68b227d 384 hbtEvent->SetZDCEMEnergy(tTotalPt);
385 // if (esdtrack->GetStatus() & AliESDtrack::kTPCrefit)
386 // if (esdtrack->GetTPCNcls() > 80)
387 // if (esdtrack->GetTPCchi2()/esdtrack->GetTPCNcls() < 6.0)
388 // if (esdtrack->GetConstrainedParam())
389 // if (fabs(esdtrack->GetConstrainedParam()->Eta()) < 0.5)
390 // if (esdtrack->GetConstrainedParam()->Pt() < 1.0) {
391 // if (esdtrack->GetSign() > 0)
392 // tNormMultPos++;
393 // else if (esdtrack->GetSign() < 0)
394 // tNormMultNeg--;
395 // }
396
397 // If reading ITS-only tracks, reject all with TPC
398 if (fTrackType == kITSOnly) {
399 if (esdtrack->GetStatus() & AliESDtrack::kTPCrefit) continue;
400 if (!(esdtrack->GetStatus() & AliESDtrack::kITSrefit)) continue;
401 if (esdtrack->GetStatus() & AliESDtrack::kTPCin) continue;
402 UChar_t iclm = esdtrack->GetITSClusterMap();
403 Int_t incls = 0;
404 for (int iter=0; iter<6; iter++) if (iclm&(1<<iter)) incls++;
405 if (incls<=3) {
406 if (Debug()>1) cout << "Rejecting track with " << incls << " clusters" << endl;
407 continue;
76ce4b5b 408 }
a68b227d 409 }
76ce4b5b 410
411
412
a68b227d 413 AliFemtoTrack* trackCopy = new AliFemtoTrack();
414 trackCopy->SetCharge((short)esdtrack->GetSign());
76ce4b5b 415
a68b227d 416 //in aliroot we have AliPID
417 //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
418 //we use only 5 first
419 double esdpid[5];
420 esdtrack->GetESDpid(esdpid);
421 trackCopy->SetPidProbElectron(esdpid[0]);
422 trackCopy->SetPidProbMuon(esdpid[1]);
423 trackCopy->SetPidProbPion(esdpid[2]);
424 trackCopy->SetPidProbKaon(esdpid[3]);
425 trackCopy->SetPidProbProton(esdpid[4]);
76ce4b5b 426
a68b227d 427 esdpid[0] = -100000.0;
428 esdpid[1] = -100000.0;
429 esdpid[2] = -100000.0;
430 esdpid[3] = -100000.0;
431 esdpid[4] = -100000.0;
76ce4b5b 432
a68b227d 433 double tTOF = 0.0;
76ce4b5b 434
a68b227d 435 if (esdtrack->GetStatus()&AliESDtrack::kTOFpid) {
436 tTOF = esdtrack->GetTOFsignal();
437 esdtrack->GetIntegratedTimes(esdpid);
438 }
76ce4b5b 439
a68b227d 440 trackCopy->SetTofExpectedTimes(tTOF-esdpid[2], tTOF-esdpid[3], tTOF-esdpid[4]);
76ce4b5b 441
442
a68b227d 443 ////// TPC ////////////////////////////////////////////
444 float nsigmaTPCK=-1000.;
445 float nsigmaTPCPi=-1000.;
446 float nsigmaTPCP=-1000.;
76ce4b5b 447
a68b227d 448 if ((fESDpid) && (esdtrack->IsOn(AliESDtrack::kTPCpid))){
449 nsigmaTPCK = fESDpid->NumberOfSigmasTPC(esdtrack,AliPID::kKaon);
450 nsigmaTPCPi = fESDpid->NumberOfSigmasTPC(esdtrack,AliPID::kPion);
451 nsigmaTPCP = fESDpid->NumberOfSigmasTPC(esdtrack,AliPID::kProton);
76ce4b5b 452
a68b227d 453 }
454 trackCopy->SetNSigmaTPCPi(nsigmaTPCPi);
455 trackCopy->SetNSigmaTPCK(nsigmaTPCK);
456 trackCopy->SetNSigmaTPCP(nsigmaTPCP);
76ce4b5b 457
a68b227d 458 ///// TOF ///////////////////////////////////////////////
76ce4b5b 459
a68b227d 460 float vp=-1000.;
461 float nsigmaTOFPi=-1000.;
462 float nsigmaTOFK=-1000.;
463 float nsigmaTOFP=-1000.;
76ce4b5b 464
a68b227d 465 if ((esdtrack->GetStatus()&AliESDtrack::kTOFpid) &&
466 (esdtrack->GetStatus()&AliESDtrack::kTOFout) &&
467 (esdtrack->GetStatus()&AliESDtrack::kTIME))
76ce4b5b 468 {
469
470 //if ((esdtrack->GetStatus()&AliESDtrack::kTOFpid) &&
471 //(esdtrack->GetStatus()&AliESDtrack::kTOFout) &&
472 //(esdtrack->GetStatus()&AliESDtrack::kTIME)){
473 // collect info from ESDpid class
48178c0e 474
a68b227d 475 if ((fESDpid) && (esdtrack->IsOn(AliESDtrack::kTOFpid))) {
48178c0e 476
477
76ce4b5b 478 double tZero = fESDpid->GetTOFResponse().GetStartTime(esdtrack->P());
48178c0e 479
76ce4b5b 480 nsigmaTOFPi = fESDpid->NumberOfSigmasTOF(esdtrack,AliPID::kPion,tZero);
481 nsigmaTOFK = fESDpid->NumberOfSigmasTOF(esdtrack,AliPID::kKaon,tZero);
482 nsigmaTOFP = fESDpid->NumberOfSigmasTOF(esdtrack,AliPID::kProton,tZero);
48178c0e 483
76ce4b5b 484 Double_t len=esdtrack->GetIntegratedLength();
485 Double_t tof=esdtrack->GetTOFsignal();
486 if(tof > 0.) vp=len/tof/0.03;
487 }
488 }
489
a68b227d 490 trackCopy->SetVTOF(vp);
491 trackCopy->SetNSigmaTOFPi(nsigmaTOFPi);
492 trackCopy->SetNSigmaTOFK(nsigmaTOFK);
493 trackCopy->SetNSigmaTOFP(nsigmaTOFP);
76ce4b5b 494
48178c0e 495
a68b227d 496 double pxyz[3];
497 double rxyz[3];
498 double impact[2];
499 double covimpact[3];
48178c0e 500
a68b227d 501 if (fUseTPCOnly) {
502 if (!esdtrack->GetTPCInnerParam()) {
503 cout << "No TPC inner param !" << endl;
504 delete trackCopy;
505 continue;
506 }
76ce4b5b 507
a68b227d 508 AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
509 param->GetXYZ(rxyz);
510 param->PropagateToDCA(fEvent->GetPrimaryVertexTPC(), (fEvent->GetMagneticField()), 10000, impact, covimpact);
511 param->GetPxPyPz(pxyz);//reading noconstarined momentum
512
513 if (fReadInner == true) {
514 AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
515 tInfo->SetPDGPid(211);
516 tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]);
517 tInfo->SetMass(0.13957);
518 // tInfo->SetEmissionPoint(rxyz[0], rxyz[1], rxyz[2], 0.0);
519 // tInfo->SetEmissionPoint(fV1[0], fV1[1], fV1[2], 0.0);
520 tInfo->SetEmissionPoint(rxyz[0]-fV1[0], rxyz[1]-fV1[1], rxyz[2]-fV1[2], 0.0);
521 trackCopy->SetHiddenInfo(tInfo);
76ce4b5b 522 }
76ce4b5b 523
a68b227d 524 if (fRotateToEventPlane) {
525 double tPhi = TMath::ATan2(pxyz[1], pxyz[0]);
526 double tRad = TMath::Hypot(pxyz[0], pxyz[1]);
76ce4b5b 527
a68b227d 528 pxyz[0] = tRad*TMath::Cos(tPhi - tReactionPlane);
529 pxyz[1] = tRad*TMath::Sin(tPhi - tReactionPlane);
530 }
76ce4b5b 531
a68b227d 532 AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
533 if (v.Mag() < 0.0001) {
534 // cout << "Found 0 momentum ???? " << pxyz[0] << " " << pxyz[1] << " " << pxyz[2] << endl;
535 delete trackCopy;
536 continue;
537 }
538 trackCopy->SetP(v);//setting momentum
539 trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
540
541 const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
542 const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
543 //setting helix I do not if it is ok
544 AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
545 trackCopy->SetHelix(helix);
546
547 //some stuff which could be useful
548 trackCopy->SetImpactD(impact[0]);
549 trackCopy->SetImpactZ(impact[1]);
550 trackCopy->SetCdd(covimpact[0]);
551 trackCopy->SetCdz(covimpact[1]);
552 trackCopy->SetCzz(covimpact[2]);
553 trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact, covimpact));
554
555 delete param;
556 }
557 else {
558 if (fReadInner == true) {
559
560 if (esdtrack->GetTPCInnerParam()) {
561 AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
562 trackCopy->SetInnerMomentum(esdtrack->GetTPCmomentum());
563 param->GetXYZ(rxyz);
564 // param->PropagateToDCA(fEvent->GetPrimaryVertex(), (fEvent->GetMagneticField()), 10000);
565 param->GetPxPyPz(pxyz);//reading noconstarined momentum
566 delete param;
567
568 AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
569 tInfo->SetPDGPid(211);
570 tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]);
571 tInfo->SetMass(0.13957);
572 // tInfo->SetEmissionPoint(rxyz[0], rxyz[1], rxyz[2], 0.0);
573 //tInfo->SetEmissionPoint(fV1[0], fV1[1], fV1[2], 0.0);
574 tInfo->SetEmissionPoint(rxyz[0]-fV1[0], rxyz[1]-fV1[1], rxyz[2]-fV1[2], 0.0);
575 trackCopy->SetHiddenInfo(tInfo);
576 }
76ce4b5b 577 }
48178c0e 578
48178c0e 579
a68b227d 580 if(fTrackType == kGlobal) {
581 if (fConstrained==true)
582 tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
583 else
584 tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
585 }
586 else if (fTrackType == kTPCOnly) {
587 if (esdtrack->GetTPCInnerParam())
588 esdtrack->GetTPCInnerParam()->GetPxPyPz(pxyz);
589 else {
590 delete trackCopy;
591 continue;
592 }
593 }
594 else if (fTrackType == kITSOnly) {
595 if (fConstrained==true)
596 tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
597 else
598 tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
599 }
600
601 if (fRotateToEventPlane) {
602 double tPhi = TMath::ATan2(pxyz[1], pxyz[0]);
603 double tRad = TMath::Hypot(pxyz[0], pxyz[1]);
76ce4b5b 604
a68b227d 605 pxyz[0] = tRad*TMath::Cos(tPhi - tReactionPlane);
606 pxyz[1] = tRad*TMath::Sin(tPhi - tReactionPlane);
607 }
76ce4b5b 608
a68b227d 609 AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
610 if (v.Mag() < 0.0001) {
611 // cout << "Found 0 momentum ???? " << pxyz[0] << " " << pxyz[1] << " " << pxyz[2] << endl;
612 delete trackCopy;
613 continue;
614 }
76ce4b5b 615
a68b227d 616 trackCopy->SetP(v);//setting momentum
617 trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
618 const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
619 const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
620 //setting helix I do not if it is ok
621 AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
622 trackCopy->SetHelix(helix);
623
624 //some stuff which could be useful
625 float imp[2];
626 float cim[3];
627 esdtrack->GetImpactParameters(imp,cim);
628
629 impact[0] = imp[0];
630 impact[1] = imp[1];
631 covimpact[0] = cim[0];
632 covimpact[1] = cim[1];
633 covimpact[2] = cim[2];
634
635 trackCopy->SetImpactD(impact[0]);
636 trackCopy->SetImpactZ(impact[1]);
637 trackCopy->SetCdd(covimpact[0]);
638 trackCopy->SetCdz(covimpact[1]);
639 trackCopy->SetCzz(covimpact[2]);
640 trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact,covimpact));
641 }
76ce4b5b 642
a68b227d 643 trackCopy->SetTrackId(esdtrack->GetID());
644 trackCopy->SetFlags(esdtrack->GetStatus());
645 trackCopy->SetLabel(esdtrack->GetLabel());
76ce4b5b 646
a68b227d 647 trackCopy->SetITSchi2(esdtrack->GetITSchi2());
648 trackCopy->SetITSncls(esdtrack->GetNcls(0));
649 trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());
650 trackCopy->SetTPCncls(esdtrack->GetTPCNcls());
651 trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());
652 trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class
653 trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma());
654 trackCopy->SetTPCsignal(esdtrack->GetTPCsignal());
76ce4b5b 655
a68b227d 656 trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap());
657 trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap());
76ce4b5b 658
a68b227d 659 double xtpc[3];
660 esdtrack->GetInnerXYZ(xtpc);
661 xtpc[2] -= fV1[2];
662 trackCopy->SetNominalTPCEntrancePoint(xtpc);
76ce4b5b 663
a68b227d 664 esdtrack->GetOuterXYZ(xtpc);
665 xtpc[2] -= fV1[2];
666 trackCopy->SetNominalTPCExitPoint(xtpc);
48178c0e 667
a68b227d 668 int indexes[3];
669 for (int ik=0; ik<3; ik++) {
670 indexes[ik] = esdtrack->GetKinkIndex(ik);
671 }
672 trackCopy->SetKinkIndexes(indexes);
48178c0e 673
a68b227d 674 for (int ii=0; ii<6; ii++){
675 trackCopy->SetITSHitOnLayer(ii,esdtrack->HasPointOnITSLayer(ii));
676 }
48178c0e 677
48178c0e 678
a68b227d 679 // Fill the hidden information with the simulated data
680 if (TMath::Abs(esdtrack->GetLabel()) < fStack->GetNtrack()) {
681 TParticle *tPart = fStack->Particle(TMath::Abs(esdtrack->GetLabel()));
48178c0e 682
a68b227d 683 // Check the mother information
1dd4ea76 684
a68b227d 685 // Using the new way of storing the freeze-out information
686 // Final state particle is stored twice on the stack
687 // one copy (mother) is stored with original freeze-out information
688 // and is not tracked
689 // the other one (daughter) is stored with primary vertex position
690 // and is tracked
1dd4ea76 691
a68b227d 692 // Freeze-out coordinates
693 double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
694 fpx = tPart->Vx() - fV1[0];
695 fpy = tPart->Vy() - fV1[1];
696 fpz = tPart->Vz() - fV1[2];
697 fpt = tPart->T();
7eb9f5d0 698
a68b227d 699 AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
700 tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
701 trackCopy->SetGlobalEmissionPoint(fpx, fpy, fpz);
7eb9f5d0 702
a68b227d 703 fpx *= 1e13;
704 fpy *= 1e13;
705 fpz *= 1e13;
706 fpt *= 1e13;
707
708
709 if (fOnlyPrimaries && !fStack->IsPhysicalPrimary(TMath::Abs(esdtrack->GetLabel()))){
710 delete trackCopy;
711 continue;
712 }
713
714 // fillDCA
715
716 if (TMath::Abs(impact[0]) > 0.001) {
717 if (fStack->IsPhysicalPrimary(TMath::Abs(esdtrack->GetLabel()))){
718 trackCopy->SetImpactDprim(impact[0]);
719 //cout << "prim" << endl;
720
721 }
722 else if (fStack->IsSecondaryFromWeakDecay(TMath::Abs(esdtrack->GetLabel()))) {
723 trackCopy->SetImpactDweak(impact[0]);
724 //cout << "wea" << endl;
725 }
726 else if (fStack->IsSecondaryFromMaterial(TMath::Abs(esdtrack->GetLabel()))) {
727 trackCopy->SetImpactDmat(impact[0]);
728 //cout << "mat" << endl;
729 }
730 }
731 // end fillDCA
732
733 // cout << "Looking for mother ids " << endl;
734 if (motherids[TMath::Abs(esdtrack->GetLabel())]>0) {
735 // cout << "Got mother id" << endl;
736 TParticle *mother = fStack->Particle(motherids[TMath::Abs(esdtrack->GetLabel())]);
737 // Check if this is the same particle stored twice on the stack
738 if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
739 // It is the same particle
740 // Read in the original freeze-out information
741 // and convert it from to [fm]
742
743 // EPOS style
744 fpx = mother->Vx()*1e13*0.197327;
745 fpy = mother->Vy()*1e13*0.197327;
746 fpz = mother->Vz()*1e13*0.197327;
747 fpt = mother->T() *1e13*0.197327;
748
749
750 // Therminator style
76ce4b5b 751// fpx = mother->Vx()*1e13;
752// fpy = mother->Vy()*1e13;
753// fpz = mother->Vz()*1e13;
754// fpt = mother->T() *1e13*3e10;
48178c0e 755
a68b227d 756 }
76ce4b5b 757 }
48178c0e 758
a68b227d 759 if (fRotateToEventPlane) {
760 double tPhi = TMath::ATan2(fpy, fpx);
761 double tRad = TMath::Hypot(fpx, fpy);
48178c0e 762
a68b227d 763 fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
764 fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
48178c0e 765 }
766
a68b227d 767 tInfo->SetPDGPid(tPart->GetPdgCode());
768 trackCopy->SetPDGPid(tPart->GetPdgCode());
48178c0e 769
a68b227d 770 if (fRotateToEventPlane) {
771 double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
772 double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
773
774 tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
775 tRad*TMath::Sin(tPhi - tReactionPlane),
776 tPart->Pz());
777 trackCopy->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
778 tRad*TMath::Sin(tPhi - tReactionPlane),
779 tPart->Pz());
780 }
781 else
782 {
783 tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
784 trackCopy->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
785 }
786 Double_t mass2 = (tPart->Energy() *tPart->Energy() -
787 tPart->Px()*tPart->Px() -
788 tPart->Py()*tPart->Py() -
789 tPart->Pz()*tPart->Pz());
790 if (mass2>0.0)
791 {
792 tInfo->SetMass(TMath::Sqrt(mass2));
793 trackCopy->SetMass(TMath::Sqrt(mass2));
794 }
76ce4b5b 795 else
a68b227d 796 {
797 tInfo->SetMass(0.0);
798 trackCopy->SetMass(0.0);
799 }
800
801 tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
802 trackCopy->SetEmissionPoint(fpx, fpy, fpz, fpt);
803 trackCopy->SetHiddenInfo(tInfo);
804 }
805 else {
806 AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
807 tInfo->SetMass(0.0);
808 double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
809 fpx = fV1[0]*1e13;
810 fpy = fV1[1]*1e13;
811 fpz = fV1[2]*1e13;
812 fpt = 0.0;
813 tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
814
815 tInfo->SetTrueMomentum(pxyz[0],pxyz[1],pxyz[2]);
816
817 trackCopy->SetHiddenInfo(tInfo);
818 }
819 // cout << "Got freeze-out " << fpx << " " << fpy << " " << fpz << " " << fpt << " " << mass2 << " " << tPart->GetPdgCode() << endl;
820
821 //decision if we want this track
822 //if we using diffrent labels we want that this label was use for first time
823 //if we use hidden info we want to have match between sim data and ESD
824
825 bool trackAccept = true;
c0869e74 826 if (isKaonAnalysis == true && TMath::Abs(trackCopy->GetPDGPid()) != 321) {
a68b227d 827 trackAccept = false;
828 }
d2eafac6 829 if (isProtonAnalysis == true && TMath::Abs(trackCopy->GetPDGPid()) != 2212) {
830 trackAccept = false;
831 }
c0869e74 832 if (isPionAnalysis == true && TMath::Abs(trackCopy->GetPDGPid()) != 211) {
89033e57 833 trackAccept = false;
834 }
48178c0e 835
a68b227d 836 if (tGoodMomentum==true && trackAccept == true)
837 {
838
839 hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
840 realnofTracks++;//real number of tracks
841 // delete trackCopy;
842 }
843 else
844 {
845 delete trackCopy;
76ce4b5b 846 }
847
a68b227d 848 }
849
76ce4b5b 850 delete [] motherids;
48178c0e 851
852 hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
76ce4b5b 853
854 AliCentrality *cent = fEvent->GetCentrality();
855 if (cent) {
856 hbtEvent->SetCentralityV0(cent->GetCentralityPercentile("V0M"));
18757d69 857 hbtEvent->SetCentralityV0A(cent->GetCentralityPercentile("V0A"));
858 hbtEvent->SetCentralityV0C(cent->GetCentralityPercentile("V0C"));
c863b24a 859 hbtEvent->SetCentralityZNA(cent->GetCentralityPercentile("ZNA"));
18757d69 860 hbtEvent->SetCentralityZNC(cent->GetCentralityPercentile("ZNC"));
c863b24a 861 hbtEvent->SetCentralityCL1(cent->GetCentralityPercentile("CL1"));
18757d69 862 hbtEvent->SetCentralityCL0(cent->GetCentralityPercentile("CL0"));
863 hbtEvent->SetCentralityTKL(cent->GetCentralityPercentile("TKL"));
864 hbtEvent->SetCentralityFMD(cent->GetCentralityPercentile("FMD"));
865 hbtEvent->SetCentralityFMD(cent->GetCentralityPercentile("NPA"));
866 // hbtEvent->SetCentralitySPD1(cent->GetCentralityPercentile("CL1"));
867 hbtEvent->SetCentralityTrk(cent->GetCentralityPercentile("TRK"));
868 hbtEvent->SetCentralityCND(cent->GetCentralityPercentile("CND"));
76ce4b5b 869
870 if (Debug()>1) printf(" FemtoReader Got Event with %f %f %f %f\n", cent->GetCentralityPercentile("V0M"), 0.0, cent->GetCentralityPercentile("CL1"), 0.0);
871 }
872
48178c0e 873 if (fEstEventMult == kGlobalCount)
76ce4b5b 874 hbtEvent->SetNormalizedMult(tNormMult);
7eb9f5d0 875 else if (fEstEventMult == kReferenceITSTPC)
876 hbtEvent->SetNormalizedMult(AliESDtrackCuts::GetReferenceMultiplicity(fEvent,AliESDtrackCuts::kTrackletsITSTPC,1.2));
877 else if(fEstEventMult == kReferenceITSSA)
878 hbtEvent->SetNormalizedMult(AliESDtrackCuts::GetReferenceMultiplicity(fEvent,AliESDtrackCuts::kTrackletsITSSA,1.2));
879 else if(fEstEventMult == kReferenceTracklets)
880 hbtEvent->SetNormalizedMult(AliESDtrackCuts::GetReferenceMultiplicity(fEvent,AliESDtrackCuts::kTracklets,1.2));
76ce4b5b 881 else if (fEstEventMult == kSPDLayer1)
882 hbtEvent->SetNormalizedMult(fEvent->GetMultiplicity()->GetNumberOfITSClusters(1));
1938eadf 883 else if (fEstEventMult == kVZERO)
a68b227d 884 {
885 Float_t multV0 = 0;
886 for (Int_t i=0; i<64; i++)
887 multV0 += fEvent->GetVZEROData()->GetMultiplicity(i);
888 hbtEvent->SetNormalizedMult(multV0);
889 }
1938eadf 890 else if (fEstEventMult == kCentrality) {
76ce4b5b 891 // centrality between 0 (central) and 1 (very peripheral)
892
893 if (cent) {
894 if (cent->GetCentralityPercentile("V0M") < 0.00001)
a68b227d 895 hbtEvent->SetNormalizedMult(-1);
76ce4b5b 896 else
a68b227d 897 hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("V0M")));
48178c0e 898 if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
a68b227d 899 10.0*cent->GetCentralityPercentile("V0M"), lrint(10.0*cent->GetCentralityPercentile("V0M")));
76ce4b5b 900 }
901 }
18757d69 902 else if (fEstEventMult == kCentralityV0A) {
903 // centrality between 0 (central) and 1 (very peripheral)
904
905 if (cent) {
906 if (cent->GetCentralityPercentile("V0A") < 0.00001)
a68b227d 907 hbtEvent->SetNormalizedMult(-1);
18757d69 908 else
a68b227d 909 hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("V0A")));
48178c0e 910 if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
a68b227d 911 10.0*cent->GetCentralityPercentile("V0A"), lrint(10.0*cent->GetCentralityPercentile("V0A")));
18757d69 912 }
913 }
914 else if (fEstEventMult == kCentralityV0C) {
915 // centrality between 0 (central) and 1 (very peripheral)
916
917 if (cent) {
918 if (cent->GetCentralityPercentile("V0C") < 0.00001)
a68b227d 919 hbtEvent->SetNormalizedMult(-1);
18757d69 920 else
a68b227d 921 hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("V0C")));
48178c0e 922 if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
a68b227d 923 10.0*cent->GetCentralityPercentile("V0C"), lrint(10.0*cent->GetCentralityPercentile("V0C")));
18757d69 924 }
925 }
c863b24a 926 else if (fEstEventMult == kCentralityZNA) {
927 // centrality between 0 (central) and 1 (very peripheral)
928
929 if (cent) {
930 if (cent->GetCentralityPercentile("ZNA") < 0.00001)
a68b227d 931 hbtEvent->SetNormalizedMult(-1);
c863b24a 932 else
a68b227d 933 hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("ZNA")));
48178c0e 934 if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
a68b227d 935 10.0*cent->GetCentralityPercentile("ZNA"), lrint(10.0*cent->GetCentralityPercentile("ZNA")));
c863b24a 936 }
937 }
18757d69 938 else if (fEstEventMult == kCentralityZNC) {
939 // centrality between 0 (central) and 1 (very peripheral)
940
941 if (cent) {
942 if (cent->GetCentralityPercentile("ZNC") < 0.00001)
a68b227d 943 hbtEvent->SetNormalizedMult(-1);
18757d69 944 else
a68b227d 945 hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("ZNC")));
48178c0e 946 if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
a68b227d 947 10.0*cent->GetCentralityPercentile("ZNC"), lrint(10.0*cent->GetCentralityPercentile("ZNC")));
18757d69 948 }
949 }
c863b24a 950 else if (fEstEventMult == kCentralityCL1) {
951 // centrality between 0 (central) and 1 (very peripheral)
952
953 if (cent) {
954 if (cent->GetCentralityPercentile("CL1") < 0.00001)
a68b227d 955 hbtEvent->SetNormalizedMult(-1);
c863b24a 956 else
a68b227d 957 hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("CL1")));
48178c0e 958 if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
a68b227d 959 10.0*cent->GetCentralityPercentile("CL1"), lrint(10.0*cent->GetCentralityPercentile("CL1")));
c863b24a 960 }
961 }
18757d69 962 else if (fEstEventMult == kCentralityCL0) {
963 // centrality between 0 (central) and 1 (very peripheral)
964
965 if (cent) {
966 if (cent->GetCentralityPercentile("CL0") < 0.00001)
a68b227d 967 hbtEvent->SetNormalizedMult(-1);
18757d69 968 else
a68b227d 969 hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("CL0")));
48178c0e 970 if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
a68b227d 971 10.0*cent->GetCentralityPercentile("CL0"), lrint(10.0*cent->GetCentralityPercentile("CL0")));
18757d69 972 }
973 }
607aa748 974 else if (fEstEventMult == kCentralityTRK) {
975 // centrality between 0 (central) and 1 (very peripheral)
976
977 if (cent) {
978 if (cent->GetCentralityPercentile("TRK") < 0.00001)
a68b227d 979 hbtEvent->SetNormalizedMult(-1);
607aa748 980 else
a68b227d 981 hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("TRK")));
48178c0e 982 if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
a68b227d 983 10.0*cent->GetCentralityPercentile("TRK"), lrint(10.0*cent->GetCentralityPercentile("TRK")));
607aa748 984 }
985 }
18757d69 986 else if (fEstEventMult == kCentralityTKL) {
987 // centrality between 0 (central) and 1 (very peripheral)
988
989 if (cent) {
990 if (cent->GetCentralityPercentile("TKL") < 0.00001)
a68b227d 991 hbtEvent->SetNormalizedMult(-1);
18757d69 992 else
a68b227d 993 hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("TKL")));
48178c0e 994 if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
a68b227d 995 10.0*cent->GetCentralityPercentile("TKL"), lrint(10.0*cent->GetCentralityPercentile("TKL")));
18757d69 996 }
997 }
998 else if (fEstEventMult == kCentralityNPA) {
999 // centrality between 0 (central) and 1 (very peripheral)
1000
1001 if (cent) {
1002 if (cent->GetCentralityPercentile("NPA") < 0.00001)
a68b227d 1003 hbtEvent->SetNormalizedMult(-1);
18757d69 1004 else
a68b227d 1005 hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("NPA")));
48178c0e 1006 if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
a68b227d 1007 10.0*cent->GetCentralityPercentile("NPA"), lrint(10.0*cent->GetCentralityPercentile("NPA")));
18757d69 1008 }
1009 }
607aa748 1010 else if (fEstEventMult == kCentralityCND) {
1011 // centrality between 0 (central) and 1 (very peripheral)
1012
1013 if (cent) {
1014 if (cent->GetCentralityPercentile("CND") < 0.00001)
a68b227d 1015 hbtEvent->SetNormalizedMult(-1);
607aa748 1016 else
a68b227d 1017 hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("CND")));
48178c0e 1018 if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
a68b227d 1019 10.0*cent->GetCentralityPercentile("CND"), lrint(10.0*cent->GetCentralityPercentile("CND")));
607aa748 1020 }
1021 }
18757d69 1022 else if (fEstEventMult == kCentralityFMD) {
1023 // centrality between 0 (central) and 1 (very peripheral)
1024
1025 if (cent) {
1026 if (cent->GetCentralityPercentile("FMD") < 0.00001)
a68b227d 1027 hbtEvent->SetNormalizedMult(-1);
18757d69 1028 else
a68b227d 1029 hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("FMD")));
48178c0e 1030 if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
a68b227d 1031 10.0*cent->GetCentralityPercentile("FMD"), lrint(10.0*cent->GetCentralityPercentile("FMD")));
18757d69 1032 }
1033 }
607aa748 1034
76ce4b5b 1035
1036 if (tNormMultPos > tNormMultNeg)
1037 hbtEvent->SetZDCParticipants(tNormMultPos);
1038 else
1039 hbtEvent->SetZDCParticipants(tNormMultNeg);
48178c0e 1040
76ce4b5b 1041 AliEventplane* ep = fEvent->GetEventplane();
1042 if (ep) {
1043 hbtEvent->SetEP(ep);
1044 hbtEvent->SetReactionPlaneAngle(ep->GetEventplane("Q"));
1045 }
1046
ce7b3d98 1047 if(fReadV0)
a68b227d 1048 {
1049 for (Int_t i = 0; i < fEvent->GetNumberOfV0s(); i++) {
1050 AliESDv0* esdv0 = fEvent->GetV0(i);
1051 if (!esdv0) continue;
1052 //if(esdv0->GetNDaughters()>2) continue;
1053 //if(esdv0->GetNProngs()>2) continue;
1054 if(esdv0->Charge()!=0) continue;
1055 AliESDtrack *trackPos = fEvent->GetTrack(esdv0->GetPindex());
1056 if(!trackPos) continue;
1057 AliESDtrack *trackNeg = fEvent->GetTrack(esdv0->GetNindex());
1058 if(!trackNeg) continue;
1059 if(trackPos->Charge()==trackNeg->Charge()) continue;
1060 AliFemtoV0* trackCopyV0 = new AliFemtoV0();
1061 CopyESDtoFemtoV0(esdv0, trackCopyV0, fEvent);
1062 hbtEvent->V0Collection()->push_back(trackCopyV0);
1063 //cout<<"Pushback v0 to v0collection"<<endl;
ce7b3d98 1064 }
a68b227d 1065 }
ce7b3d98 1066
76ce4b5b 1067
48178c0e 1068 fCurEvent++;
76ce4b5b 1069 // cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
48178c0e 1070 return hbtEvent;
76ce4b5b 1071}
1072//___________________
1073void AliFemtoEventReaderESDChainKine::SetESDSource(AliESDEvent *aESD)
1074{
1075 // The chain loads the ESD for us
1076 // You must provide the address where it can be found
1077 fEvent = aESD;
1078}
1079//___________________
1080void AliFemtoEventReaderESDChainKine::SetStackSource(AliStack *aStack)
1081{
1082 // The chain loads the stack for us
1083 // You must provide the address where it can be found
1084 fStack = aStack;
1085}
1086//___________________
1087void AliFemtoEventReaderESDChainKine::SetGenEventHeader(AliGenEventHeader *aGenHeader)
1088{
1089 // The chain loads the generator event header for us
1090 // You must provide the address where it can be found
1091 fGenHeader = aGenHeader;
1092}
1093void AliFemtoEventReaderESDChainKine::SetRotateToEventPlane(short dorotate)
1094{
1095 fRotateToEventPlane=dorotate;
1096}
1097Float_t AliFemtoEventReaderESDChainKine::GetSigmaToVertex(double *impact, double *covar)
1098{
1099 // Calculates the number of sigma to the vertex.
1100
1101 Float_t b[2];
1102 Float_t bRes[2];
1103 Float_t bCov[3];
1104
1105 b[0] = impact[0];
1106 b[1] = impact[1];
1107 bCov[0] = covar[0];
1108 bCov[1] = covar[1];
1109 bCov[2] = covar[2];
1110
1111 bRes[0] = TMath::Sqrt(bCov[0]);
1112 bRes[1] = TMath::Sqrt(bCov[2]);
1113
1114 // -----------------------------------
1115 // How to get to a n-sigma cut?
1116 //
1117 // The accumulated statistics from 0 to d is
1118 //
1119 // -> Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
1120 // -> 1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
1121 //
1122 // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
1123 // Can this be expressed in a different way?
1124
1125 if (bRes[0] == 0 || bRes[1] ==0)
1126 return -1;
1127
1128 Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
1129
1130 // stupid rounding problem screws up everything:
1131 // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
1132 if (TMath::Exp(-d * d / 2) < 1e-10)
1133 return 1000;
1134
1135 d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
1136 return d;
1137}
1138
1139void AliFemtoEventReaderESDChainKine::SetReadTrackType(ReadTrackType aType)
1140{
1141 fTrackType = aType;
1142}
1143
ce7b3d98 1144void AliFemtoEventReaderESDChainKine::SetReadV0(bool a)
1145{
1146 fReadV0 = a;
1147}
1148
48178c0e 1149void AliFemtoEventReaderESDChainKine::SetKaonAnalysis(Bool_t a)
1150{
1151 isKaonAnalysis = a;
1152}
1153
d2eafac6 1154void AliFemtoEventReaderESDChainKine::SetProtonAnalysis(Bool_t a)
1155{
1156 isProtonAnalysis = a;
1157}
1158
89033e57 1159void AliFemtoEventReaderESDChainKine::SetPionAnalysis(Bool_t a)
1160{
1161 isPionAnalysis = a;
1162}
1163
1dd4ea76 1164void AliFemtoEventReaderESDChainKine::SetOnlyPrimaries(Bool_t a)
1165{
1166 fOnlyPrimaries = a;
1167}
1168
ce7b3d98 1169void AliFemtoEventReaderESDChainKine::CopyESDtoFemtoV0(AliESDv0 *tESDv0, AliFemtoV0 *tFemtoV0, AliESDEvent *tESDevent)
1170{
1171 double fPrimaryVtxPosition[3];
1172 tESDevent->GetPrimaryVertex()->GetXYZ(fPrimaryVtxPosition);
1173 tFemtoV0->SetdcaV0ToPrimVertex(tESDv0->GetD(fPrimaryVtxPosition[0],fPrimaryVtxPosition[1],fPrimaryVtxPosition[2]));
1174
1175 //tFemtoV0->SetdecayLengthV0(tESDv0->DecayLengthV0()); //wrocic do tego
1176 //tFemtoV0->SetdecayVertexV0X(tESDv0->DecayVertexV0X());
1177 //tFemtoV0->SetdecayVertexV0Y(tESDv0->DecayVertexV0Y());
1178 //tFemtoV0->SetdecayVertexV0Z(tESDv0->DecayVertexV0Z()); //nie ma w AliESDv0
1179 //AliFemtoThreeVector decayvertex(tESDv0->DecayVertexV0X(),tESDv0->DecayVertexV0Y(),tESDv0->DecayVertexV0Z());
1180 //tFemtoV0->SetdecayVertexV0(decayvertex);
1181 tFemtoV0->SetdcaV0Daughters(tESDv0->GetDcaV0Daughters());
1182 tFemtoV0->SetmomV0X(tESDv0->Px());
1183 tFemtoV0->SetmomV0Y(tESDv0->Py());
1184 tFemtoV0->SetmomV0Z(tESDv0->Pz());
1185 AliFemtoThreeVector momv0(tESDv0->Px(),tESDv0->Py(),tESDv0->Pz());
1186 tFemtoV0->SetmomV0(momv0);
1187 tFemtoV0->SetalphaV0(tESDv0->AlphaV0());
1188 tFemtoV0->SetptArmV0(tESDv0->PtArmV0());
1189 //tFemtoV0->SeteLambda(tESDv0->ELambda());
1190 //tFemtoV0->SeteK0Short(tESDv0->EK0Short());
1191 //tFemtoV0->SetePosProton(tESDv0->EPosProton());
1192 //tFemtoV0->SeteNegProton(tESDv0->ENegProton());
1193 tFemtoV0->SetmassLambda(tESDv0->GetEffMass(4,2));
1194 tFemtoV0->SetmassAntiLambda(tESDv0->GetEffMass(2,4));
1195 tFemtoV0->SetmassK0Short(tESDv0->GetEffMass(2,2));
1196 //tFemtoV0->SetrapLambda(tESDv0->RapLambda());
1197 //tFemtoV0->SetrapK0Short(tESDv0->RapK0Short());
1198 tFemtoV0->SetptV0(tESDv0->Pt());
1199 tFemtoV0->SetptotV0(tESDv0->P());
1200 tFemtoV0->SetEtaV0(tESDv0->Eta());
1201 tFemtoV0->SetPhiV0(tESDv0->Phi());
1202 tFemtoV0->SetCosPointingAngle(tESDv0->GetV0CosineOfPointingAngle(fPrimaryVtxPosition[0],fPrimaryVtxPosition[1], fPrimaryVtxPosition[2]));
1203 tFemtoV0->SetYV0(tESDv0->Y());
1204
1205 AliESDtrack *trackpos = tESDevent->GetTrack(tESDv0->GetPindex()); //AliAODTrack *trackpos = (AliAODTrack*)tESDv0->GetDaughter(0);
1206 AliESDtrack *trackneg = tESDevent->GetTrack(tESDv0->GetNindex()); //AliAODTrack *trackneg = (AliAODTrack*)tESDv0->GetDaughter(1);
1207
1208 if(trackpos && trackneg)
a68b227d 1209 {
1210 tFemtoV0->SetdcaPosToPrimVertex(TMath::Abs(trackpos->GetD(fPrimaryVtxPosition[0],fPrimaryVtxPosition[1],tESDevent->GetMagneticField())));
1211 tFemtoV0->SetdcaNegToPrimVertex(TMath::Abs(trackneg->GetD(fPrimaryVtxPosition[0],fPrimaryVtxPosition[1],tESDevent->GetMagneticField())));
1212 double MomPos[3];
1213 trackpos->PxPyPz(MomPos);
1214 tFemtoV0->SetmomPosX(MomPos[0]);
1215 tFemtoV0->SetmomPosY(MomPos[1]);
1216 tFemtoV0->SetmomPosZ(MomPos[2]);
1217 AliFemtoThreeVector mompos(MomPos[0],MomPos[1],MomPos[2]);
1218 tFemtoV0->SetmomPos(mompos);
1219
1220 double MomNeg[3];
1221 trackneg->PxPyPz(MomNeg);
1222 tFemtoV0->SetmomNegX(MomNeg[0]);
1223 tFemtoV0->SetmomNegY(MomNeg[1]);
1224 tFemtoV0->SetmomNegZ(MomNeg[2]);
1225 AliFemtoThreeVector momneg(MomNeg[0],MomNeg[1],MomNeg[2]);
1226 tFemtoV0->SetmomNeg(momneg);
1227
1228 tFemtoV0->SetptPos(trackpos->Pt());
1229 tFemtoV0->SetptotPos(trackpos->P());
1230 tFemtoV0->SetptNeg(trackneg->Pt());
1231 tFemtoV0->SetptotNeg(trackneg->P());
1232
1233 tFemtoV0->SetidNeg(trackneg->GetID());
1234 //cout<<"tESDv0->GetNegID(): "<<tESDv0->GetNegID()<<endl;
1235 //cout<<"tFemtoV0->IdNeg(): "<<tFemtoV0->IdNeg()<<endl;
1236 tFemtoV0->SetidPos(trackpos->GetID());
1237
1238 tFemtoV0->SetEtaPos(trackpos->Eta());
1239 tFemtoV0->SetEtaNeg(trackneg->Eta());
1240
1241 tFemtoV0->SetEtaPos(trackpos->Eta()); //tESDv0->PseudoRapPos()
1242 tFemtoV0->SetEtaNeg(trackneg->Eta()); //tESDv0->PseudoRapNeg()
1243 tFemtoV0->SetTPCNclsPos(trackpos->GetTPCNcls());
1244 tFemtoV0->SetTPCNclsNeg(trackneg->GetTPCNcls());
1245 tFemtoV0->SetTPCclustersPos(trackpos->GetTPCClusterMap());
1246 tFemtoV0->SetTPCclustersNeg(trackneg->GetTPCClusterMap());
1247 tFemtoV0->SetTPCsharingPos(trackpos->GetTPCSharedMap());
1248 tFemtoV0->SetTPCsharingNeg(trackneg->GetTPCSharedMap());
1249 tFemtoV0->SetNdofPos(trackpos->GetTPCchi2()/trackpos->GetTPCNcls());
1250 tFemtoV0->SetNdofNeg(trackneg->GetTPCchi2()/trackneg->GetTPCNcls());
1251 tFemtoV0->SetStatusPos(trackpos->GetStatus());
1252 tFemtoV0->SetStatusNeg(trackneg->GetStatus());
1253
1254 float bfield = 5*fMagFieldSign;
1255 float globalPositionsAtRadiiPos[9][3];
1256 GetGlobalPositionAtGlobalRadiiThroughTPC(trackpos,bfield,globalPositionsAtRadiiPos);
1257 double tpcEntrancePos[3]={globalPositionsAtRadiiPos[0][0],globalPositionsAtRadiiPos[0][1],globalPositionsAtRadiiPos[0][2]};
1258 double tpcExitPos[3]={globalPositionsAtRadiiPos[8][0],globalPositionsAtRadiiPos[8][1],globalPositionsAtRadiiPos[8][2]};
1259
1260 float globalPositionsAtRadiiNeg[9][3];
1261 GetGlobalPositionAtGlobalRadiiThroughTPC(trackneg,bfield,globalPositionsAtRadiiNeg);
1262 double tpcEntranceNeg[3]={globalPositionsAtRadiiNeg[0][0],globalPositionsAtRadiiNeg[0][1],globalPositionsAtRadiiNeg[0][2]};
1263 double tpcExitNeg[3]={globalPositionsAtRadiiNeg[8][0],globalPositionsAtRadiiNeg[8][1],globalPositionsAtRadiiNeg[8][2]};
1264
1265 AliFemtoThreeVector tmpVec;
1266 tmpVec.SetX(tpcEntrancePos[0]); tmpVec.SetX(tpcEntrancePos[1]); tmpVec.SetX(tpcEntrancePos[2]);
1267 tFemtoV0->SetNominalTpcEntrancePointPos(tmpVec);
1268
1269 tmpVec.SetX(tpcExitPos[0]); tmpVec.SetX(tpcExitPos[1]); tmpVec.SetX(tpcExitPos[2]);
1270 tFemtoV0->SetNominalTpcExitPointPos(tmpVec);
1271
1272 tmpVec.SetX(tpcEntranceNeg[0]); tmpVec.SetX(tpcEntranceNeg[1]); tmpVec.SetX(tpcEntranceNeg[2]);
1273 tFemtoV0->SetNominalTpcEntrancePointNeg(tmpVec);
1274
1275 tmpVec.SetX(tpcExitNeg[0]); tmpVec.SetX(tpcExitNeg[1]); tmpVec.SetX(tpcExitNeg[2]);
1276 tFemtoV0->SetNominalTpcExitPointNeg(tmpVec);
1277
1278 AliFemtoThreeVector vecTpcPos[9];
1279 AliFemtoThreeVector vecTpcNeg[9];
1280 for(int i=0;i<9;i++)
ce7b3d98 1281 {
a68b227d 1282 vecTpcPos[i].SetX(globalPositionsAtRadiiPos[i][0]); vecTpcPos[i].SetY(globalPositionsAtRadiiPos[i][1]); vecTpcPos[i].SetZ(globalPositionsAtRadiiPos[i][2]);
1283 vecTpcNeg[i].SetX(globalPositionsAtRadiiNeg[i][0]); vecTpcNeg[i].SetY(globalPositionsAtRadiiNeg[i][1]); vecTpcNeg[i].SetZ(globalPositionsAtRadiiNeg[i][2]);
1284 }
1285 tFemtoV0->SetNominalTpcPointPos(vecTpcPos);
1286 tFemtoV0->SetNominalTpcPointNeg(vecTpcNeg);
ce7b3d98 1287
a68b227d 1288 tFemtoV0->SetTPCMomentumPos(trackpos->GetTPCInnerParam()->P()); //trackpos->GetTPCmomentum();
1289 tFemtoV0->SetTPCMomentumNeg(trackneg->GetTPCInnerParam()->P()); //trackneg->GetTPCmomentum();
ce7b3d98 1290
a68b227d 1291 tFemtoV0->SetdedxPos(trackpos->GetTPCsignal());
1292 tFemtoV0->SetdedxNeg(trackneg->GetTPCsignal());
1293
1294
1295 if (fESDpid) {
1296 tFemtoV0->SetPosNSigmaTPCK(fESDpid->NumberOfSigmasTPC(trackpos,AliPID::kKaon));
1297 tFemtoV0->SetNegNSigmaTPCK(fESDpid->NumberOfSigmasTPC(trackneg,AliPID::kKaon));
1298 tFemtoV0->SetPosNSigmaTPCP(fESDpid->NumberOfSigmasTPC(trackpos,AliPID::kProton));
1299 tFemtoV0->SetNegNSigmaTPCP(fESDpid->NumberOfSigmasTPC(trackneg,AliPID::kProton));
1300 tFemtoV0->SetPosNSigmaTPCPi(fESDpid->NumberOfSigmasTPC(trackpos,AliPID::kPion));
1301 tFemtoV0->SetNegNSigmaTPCPi(fESDpid->NumberOfSigmasTPC(trackneg,AliPID::kPion));
1302 }
1303 else {
1304 tFemtoV0->SetPosNSigmaTPCK(-1000);
1305 tFemtoV0->SetNegNSigmaTPCK(-1000);
1306 tFemtoV0->SetPosNSigmaTPCP(-1000);
1307 tFemtoV0->SetNegNSigmaTPCP(-1000);
1308 tFemtoV0->SetPosNSigmaTPCPi(-1000);
1309 tFemtoV0->SetNegNSigmaTPCPi(-1000);
1310 }
ce7b3d98 1311
a68b227d 1312 if((tFemtoV0->StatusPos()&AliESDtrack::kTOFpid)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTOFout)==0)
1313 {
1314 if((tFemtoV0->StatusNeg()&AliESDtrack::kTOFpid)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTOFout)==0)
ce7b3d98 1315 {
1316 tFemtoV0->SetPosNSigmaTOFK(-1000);
1317 tFemtoV0->SetNegNSigmaTOFK(-1000);
1318 tFemtoV0->SetPosNSigmaTOFP(-1000);
1319 tFemtoV0->SetNegNSigmaTOFP(-1000);
1320 tFemtoV0->SetPosNSigmaTOFPi(-1000);
1321 tFemtoV0->SetNegNSigmaTOFPi(-1000);
1322 }
a68b227d 1323 }
1324 else
1325 {
1326 if (fESDpid) {
1327 tFemtoV0->SetPosNSigmaTOFK(fESDpid->NumberOfSigmasTOF(trackpos,AliPID::kKaon));
1328 tFemtoV0->SetNegNSigmaTOFK(fESDpid->NumberOfSigmasTOF(trackneg,AliPID::kKaon));
1329 tFemtoV0->SetPosNSigmaTOFP(fESDpid->NumberOfSigmasTOF(trackpos,AliPID::kProton));
1330 tFemtoV0->SetNegNSigmaTOFP(fESDpid->NumberOfSigmasTOF(trackneg,AliPID::kProton));
1331 tFemtoV0->SetPosNSigmaTOFPi(fESDpid->NumberOfSigmasTOF(trackpos,AliPID::kPion));
1332 tFemtoV0->SetNegNSigmaTOFPi(fESDpid->NumberOfSigmasTOF(trackneg,AliPID::kPion));
1333 }
1334 else {
1335 tFemtoV0->SetPosNSigmaTOFK(-1000);
1336 tFemtoV0->SetNegNSigmaTOFK(-1000);
1337 tFemtoV0->SetPosNSigmaTOFP(-1000);
1338 tFemtoV0->SetNegNSigmaTOFP(-1000);
1339 tFemtoV0->SetPosNSigmaTOFPi(-1000);
1340 tFemtoV0->SetNegNSigmaTOFPi(-1000);
1341 }
1342 }
7eb9f5d0 1343
a68b227d 1344 if ((TMath::Abs(trackpos->GetLabel()) < fStack->GetNtrack()) && (TMath::Abs(trackneg->GetLabel()) < fStack->GetNtrack())) {
7eb9f5d0 1345
a68b227d 1346 //cout<<"tESDv0->GetPdgCode(): "<<tESDv0->GetPdgCode()<<endl;
1347 // cout<<"Labels: "<<trackpos->GetLabel()<<" "<<trackneg->GetLabel()<<endl;
1348 AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
1349 //TParticle *tPart = fStack->Particle(tESDv0->GetLabel()); //zle
7eb9f5d0 1350
a68b227d 1351 int labelpos = TMath::Abs(trackpos->GetLabel());
1352 int labelneg = TMath::Abs(trackneg->GetLabel());
1353 TParticle *tPartPos = fStack->Particle(labelpos);
1354 TParticle *tPartNeg = fStack->Particle(labelneg);
48178c0e 1355
7eb9f5d0 1356
a68b227d 1357 double impactpos[2];
1358 double impactneg[2];
1359 double covimpact[3];
7eb9f5d0 1360
a68b227d 1361 AliExternalTrackParam *parampos = new AliExternalTrackParam(*trackpos->GetTPCInnerParam());
1362 parampos->PropagateToDCA(fEvent->GetPrimaryVertexTPC(), (fEvent->GetMagneticField()), 10000, impactpos, covimpact);
1363 AliExternalTrackParam *paramneg = new AliExternalTrackParam(*trackneg->GetTPCInnerParam());
1364 paramneg->PropagateToDCA(fEvent->GetPrimaryVertexTPC(), (fEvent->GetMagneticField()), 10000, impactneg, covimpact);
1365
1366
1367 // fillDCA
1368 if (TMath::Abs(impactpos[0]) > 0.001) {
1369 if (fStack->IsPhysicalPrimary(labelpos)){
1370 tFemtoV0->SetImpactDprimPos(impactpos[0]);
1371 }
1372 else if (fStack->IsSecondaryFromWeakDecay(labelpos)) {
1373 tFemtoV0->SetImpactDweakPos(impactpos[0]);
1374 //cout << "wea" << endl;
1375 }
1376 else if (fStack->IsSecondaryFromMaterial(labelpos)) {
1377 tFemtoV0->SetImpactDmatPos(impactpos[0]);
1378 //cout << "mat" << endl;
1379 }
1380 }
1381 if (TMath::Abs(impactneg[0]) > 0.001) {
1382 if (fStack->IsPhysicalPrimary(labelneg)){
1383 tFemtoV0->SetImpactDprimNeg(impactneg[0]);
1384 //cout << "prim" << endl;
1385 }
1386 else if (fStack->IsSecondaryFromWeakDecay(labelneg)) {
1387 tFemtoV0->SetImpactDweakNeg(impactneg[0]);
1388 //cout << "wea" << endl;
1389 }
1390 else if (fStack->IsSecondaryFromMaterial(labelneg)) {
1391 tFemtoV0->SetImpactDmatNeg(impactneg[0]);
1392 //cout << "mat" << endl;
1393 }
7eb9f5d0 1394
a68b227d 1395 }
1396 // end fillDCA
48178c0e 1397
7eb9f5d0 1398
a68b227d 1399 //tInfo->SetPDGPid();
1400 //tInfo->SetMass();
1401 //tInfo->SetTrueMomentum();
1402 //tInfo->SetEmissionPoint();
7eb9f5d0 1403
a68b227d 1404 // Freeze-out coordinates
1405 double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
7eb9f5d0 1406
a68b227d 1407 fpx = tPartPos->Vx() - fPrimaryVtxPosition[0];
1408 fpy = tPartPos->Vy() - fPrimaryVtxPosition[1];
1409 fpz = tPartPos->Vz() - fPrimaryVtxPosition[2];
1410 fpt = tPartPos->T();
48178c0e 1411
a68b227d 1412 fpx *= 1e13;
1413 fpy *= 1e13;
1414 fpz *= 1e13;
1415 fpt *= 1e13;
7eb9f5d0 1416
a68b227d 1417 tInfo->SetPDGPidPos(tPartPos->GetPdgCode());
1418 tInfo->SetMassPos(tPartPos->GetMass());
1419 tInfo->SetTrueMomentumPos(tPartPos->Px(),tPartPos->Py(),tPartPos->Pz());
1420 tInfo->SetEmissionPointPos(fpx,fpy,fpz,fpt);
7eb9f5d0 1421
a68b227d 1422 fpx = tPartNeg->Vx() - fPrimaryVtxPosition[0];
1423 fpy = tPartNeg->Vy() - fPrimaryVtxPosition[1];
1424 fpz = tPartNeg->Vz() - fPrimaryVtxPosition[2];
1425 fpt = tPartNeg->T();
48178c0e 1426
a68b227d 1427 fpx *= 1e13;
1428 fpy *= 1e13;
1429 fpz *= 1e13;
1430 fpt *= 1e13;
ce7b3d98 1431
a68b227d 1432 tInfo->SetPDGPidNeg(tPartNeg->GetPdgCode());
1433 tInfo->SetMassNeg(tPartNeg->GetMass());
1434 tInfo->SetTrueMomentumNeg(tPartNeg->Px(),tPartNeg->Py(),tPartNeg->Pz());
1435 tInfo->SetEmissionPointNeg(fpx,fpy,fpz,fpt);
ce7b3d98 1436
a68b227d 1437 tFemtoV0->SetHiddenInfo(tInfo);
ce7b3d98 1438 }
a68b227d 1439 }
ce7b3d98 1440 else
a68b227d 1441 {
1442 tFemtoV0->SetStatusPos(999);
1443 tFemtoV0->SetStatusNeg(999);
1444 }
ce7b3d98 1445 tFemtoV0->SetOnFlyStatusV0(tESDv0->GetOnFlyStatus());
1446}
1447
1448void AliFemtoEventReaderESDChainKine::SetMagneticFieldSign(int s)
1449{
1450 if(s>0)
1451 fMagFieldSign = 1;
1452 else if(s<0)
1453 fMagFieldSign = -1;
1454 else
1455 fMagFieldSign = 0;
1456}
1457
1458void AliFemtoEventReaderESDChainKine::GetGlobalPositionAtGlobalRadiiThroughTPC(AliESDtrack *track, Float_t bfield, Float_t globalPositionsAtRadii[9][3])
1459{
1460 // Gets the global position of the track at nine different radii in the TPC
1461 // track is the track you want to propagate
1462 // bfield is the magnetic field of your event
1463 // globalPositionsAtRadii is the array of global positions in the radii and xyz
1464
1465 // Initialize the array to something indicating there was no propagation
1466 for(Int_t i=0;i<9;i++){
1467 for(Int_t j=0;j<3;j++){
1468 globalPositionsAtRadii[i][j]=-9999.;
1469 }
1470 }
1471
1472 // Make a copy of the track to not change parameters of the track
1473 AliExternalTrackParam etp; etp.CopyFromVTrack(track);
1474 //printf("\nAfter CopyFromVTrack\n");
1475 //etp.Print();
1476
1477 // The global position of the the track
1478 Double_t xyz[3]={-9999.,-9999.,-9999.};
1479
1480 // Counter for which radius we want
1481 Int_t iR=0;
1482 // The radii at which we get the global positions
1483 // IROC (OROC) from 84.1 cm to 132.1 cm (134.6 cm to 246.6 cm)
1484 Float_t Rwanted[9]={85.,105.,125.,145.,165.,185.,205.,225.,245.};
1485 // The global radius we are at
1486 Float_t globalRadius=0;
1487
1488 // Propagation is done in local x of the track
1489 for (Float_t x = etp.GetX();x<247.;x+=1.){ // GetX returns local coordinates
1490 // Starts at the tracks fX and goes outwards. x = 245 is the outer radial limit
1491 // of the TPC when the track is straight, i.e. has inifinite pt and doesn't get bent.
1492 // If the track's momentum is smaller than infinite, it will develop a y-component, which
1493 // adds to the global radius
1494
1495 // Stop if the propagation was not succesful. This can happen for low pt tracks
1496 // that don't reach outer radii
1497 if(!etp.PropagateTo(x,bfield))break;
1498 etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1499 globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
1500
1501 // Roughly reached the radius we want
1502 if(globalRadius > Rwanted[iR]){
1503
1504 // Bigger loop has bad precision, we're nearly one centimeter too far, go back in small steps.
1505 while (globalRadius>Rwanted[iR]){
a68b227d 1506 x-=.1;
1507 // printf("propagating to x %5.2f\n",x);
1508 if(!etp.PropagateTo(x,bfield))break;
1509 etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1510 globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
ce7b3d98 1511 }
1512 //printf("At Radius:%05.2f (local x %5.2f). Setting position to x %4.1f y %4.1f z %4.1f\n",globalRadius,x,xyz[0],xyz[1],xyz[2]);
1513 globalPositionsAtRadii[iR][0]=xyz[0];
1514 globalPositionsAtRadii[iR][1]=xyz[1];
1515 globalPositionsAtRadii[iR][2]=xyz[2];
1516 // Indicate we want the next radius
1517 iR+=1;
1518 }
1519 if(iR>=8){
1520 // TPC edge reached
1521 return;
1522 }
1523 }
1524}