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