]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/FEMTOSCOPY/AliFemto/AliFemtoEventReaderStandard.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Add Standard ALICE event reader
[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemto / AliFemtoEventReaderStandard.cxx
CommitLineData
c43ed0f6 1////////////////////////////////////////////////////////////////////////////////
2// //
3// AliFemtoEventReaderStandard - the reader class for the Alice ESD, AOD //
4// the model Kinematics information tailored for the Task framework //
5// Authors: Adam Kisiel Adam.Kisiel@cern.ch //
6// //
7////////////////////////////////////////////////////////////////////////////////
8#include "AliFemtoEventReaderStandard.h"
9
10#include "TFile.h"
11#include "TTree.h"
12#include "TList.h"
13#include "TBits.h"
14#include "AliESDEvent.h"
15#include "AliESDtrack.h"
16#include "AliESDVertex.h"
17
18#include "AliFmPhysicalHelixD.h"
19#include "AliFmThreeVectorF.h"
20
21#include "SystemOfUnits.h"
22
23#include "AliFemtoEvent.h"
24
25#include "TParticle.h"
26#include "AliFemtoModelHiddenInfo.h"
27#include "AliFemtoModelGlobalHiddenInfo.h"
28#include "AliGenHijingEventHeader.h"
29#include "AliGenCocktailEventHeader.h"
30
31#include "AliAODMCHeader.h"
32#include "AliAODMCParticle.h"
33
34#include "AliVertexerTracks.h"
35
36ClassImp(AliFemtoEventReaderStandard)
37
38#if !(ST_NO_NAMESPACES)
39 using namespace units;
40#endif
41
42using namespace std;
43//____________________________
44AliFemtoEventReaderStandard::AliFemtoEventReaderStandard():
45 AliFemtoEventReader(),
46 fFileName(" "),
47 fNumberofEvent(0),
48 fCurEvent(0),
49 fCurFile(0),
50 fESDEvent(0x0),
51 fAODEvent(0x0),
52 fStack(0x0),
53 fGenHeader(0x0),
54 fInputType(kUnknown),
55 fUsePhysicsSel(kFALSE),
56 fSelect(0),
57 fTrackCuts(0x0)
58{
59 //constructor with 0 parameters , look at default settings
60}
61
62//__________________
63AliFemtoEventReaderStandard::AliFemtoEventReaderStandard(const AliFemtoEventReaderStandard& aReader):
64 AliFemtoEventReader(aReader),
65 fFileName(" "),
66 fNumberofEvent(0),
67 fCurEvent(0),
68 fCurFile(0),
69 fESDEvent(0x0),
70 fAODEvent(0x0),
71 fStack(0x0),
72 fGenHeader(0x0),
73 fInputType(kUnknown),
74 fUsePhysicsSel(kFALSE),
75 fSelect(0),
76 fTrackCuts(0x0)
77{
78 // Copy constructor
79 fCurEvent = aReader.fCurEvent;
80 fCurFile = aReader.fCurFile;
81 fESDEvent = new AliESDEvent();
82 fAODEvent = new AliAODEvent();
83 fStack = aReader.fStack;
84 fInputType = aReader.fInputType;
85 fUsePhysicsSel = aReader.fUsePhysicsSel;
86 if (fUsePhysicsSel) fSelect = new AliPhysicsSelection();
87 fTrackCuts = new AliESDtrackCuts(*(aReader.fTrackCuts));
88}
89//__________________
90AliFemtoEventReaderStandard::~AliFemtoEventReaderStandard()
91{
92 //Destructor
93 delete fESDEvent;
94 delete fAODEvent;
95 if (fTrackCuts) delete fTrackCuts;
96}
97
98//__________________
99AliFemtoEventReaderStandard& AliFemtoEventReaderStandard::operator=(const AliFemtoEventReaderStandard& aReader)
100{
101 // Assignment operator
102 if (this == &aReader)
103 return *this;
104
105 fCurEvent = aReader.fCurEvent;
106 fCurFile = aReader.fCurFile;
107 if (fESDEvent) delete fESDEvent;
108 fESDEvent = new AliESDEvent();
109 if (fAODEvent) delete fAODEvent;
110 fAODEvent = new AliAODEvent();
111 fStack = aReader.fStack;
112 fGenHeader = aReader.fGenHeader;
113 fInputType = aReader.fInputType;
114 fUsePhysicsSel = aReader.fUsePhysicsSel;
115 if (fUsePhysicsSel) fSelect = new AliPhysicsSelection();
116 if (fTrackCuts) delete fTrackCuts;
117 fTrackCuts = new AliESDtrackCuts(*(aReader.fTrackCuts));
118 return *this;
119}
120//__________________
121// Simple report
122AliFemtoString AliFemtoEventReaderStandard::Report()
123{
124 AliFemtoString temp = "\n This is the AliFemtoEventReaderStandard\n";
125 return temp;
126}
127//__________________
128AliFemtoEvent* AliFemtoEventReaderStandard::ReturnHbtEvent()
129{
130 // Get the event, read all the relevant information
131 // and fill the AliFemtoEvent class
132 // Returns a valid AliFemtoEvent
133 AliFemtoEvent *hbtEvent = 0;
134 string tFriendFileName;
135
136 hbtEvent = new AliFemtoEvent;
137
138 // Get the friend information
139 cout<<"starting to read event "<<fCurEvent<<endl;
140 if ((fInputType == kESD) || (fInputType == kESDKine)) {
141 if(fESDEvent->GetAliESDOld())fESDEvent->CopyFromOldESD();
142
143 if (fUsePhysicsSel) {
144 hbtEvent->SetIsCollisionCandidate(fSelect->IsCollisionCandidate(fESDEvent));
145 if (!(fSelect->IsCollisionCandidate(fESDEvent)))
146 printf("Event not a collision candidate\n");
147 }
148 else
149 hbtEvent->SetIsCollisionCandidate(kTRUE);
150 }
151 else {
152 hbtEvent->SetIsCollisionCandidate(kTRUE);
153 }
154
155 hbtEvent = new AliFemtoEvent;
156
157 double fV1[3];
158
159 //setting basic things
160 if ((fInputType == kESD) || (fInputType == kESDKine)) {
161 hbtEvent->SetRunNumber(fESDEvent->GetRunNumber());
162 hbtEvent->SetMagneticField(fESDEvent->GetMagneticField()*kilogauss);//to check if here is ok
163 hbtEvent->SetZDCN1Energy(fESDEvent->GetZDCN1Energy());
164 hbtEvent->SetZDCP1Energy(fESDEvent->GetZDCP1Energy());
165 hbtEvent->SetZDCN2Energy(fESDEvent->GetZDCN2Energy());
166 hbtEvent->SetZDCP2Energy(fESDEvent->GetZDCP2Energy());
167 hbtEvent->SetZDCEMEnergy(fESDEvent->GetZDCEMEnergy());
168 hbtEvent->SetZDCParticipants(fESDEvent->GetZDCParticipants());
169 hbtEvent->SetTriggerMask(fESDEvent->GetTriggerMask());
170 hbtEvent->SetTriggerCluster(fESDEvent->GetTriggerCluster());
171
172 //Vertex
173 double fVCov[6];
174 // if (fUseTPCOnly) {
175 // fESDEvent->GetPrimaryVertexTPC()->GetXYZ(fV1);
176 // fESDEvent->GetPrimaryVertexTPC()->GetCovMatrix(fVCov);
177 // if (!fESDEvent->GetPrimaryVertexTPC()->GetStatus())
178 // fVCov[4] = -1001.0;
179 // }
180 // else {
181 if (fESDEvent->GetPrimaryVertex()) {
182 fESDEvent->GetPrimaryVertex()->GetXYZ(fV1);
183 fESDEvent->GetPrimaryVertex()->GetCovMatrix(fVCov);
184
185 if (!fESDEvent->GetPrimaryVertex()->GetStatus()) {
186 // Get the vertex from SPD
187 fESDEvent->GetPrimaryVertexSPD()->GetXYZ(fV1);
188 fESDEvent->GetPrimaryVertexSPD()->GetCovMatrix(fVCov);
189
190
191 if (!fESDEvent->GetPrimaryVertexSPD()->GetStatus())
192 fVCov[4] = -1001.0;
193 else {
194 fESDEvent->GetPrimaryVertexSPD()->GetXYZ(fV1);
195 fESDEvent->GetPrimaryVertexSPD()->GetCovMatrix(fVCov);
196 }
197 }
198 }
199 else {
200 if (fESDEvent->GetPrimaryVertexSPD()) {
201 fESDEvent->GetPrimaryVertexSPD()->GetXYZ(fV1);
202 fESDEvent->GetPrimaryVertexSPD()->GetCovMatrix(fVCov);
203 }
204 }
205 if ((!fESDEvent->GetPrimaryVertex()) && (!fESDEvent->GetPrimaryVertexSPD()))
206 {
207 cout << "No vertex found !!!" << endl;
208 fV1[0] = 10000.0;
209 fV1[1] = 10000.0;
210 fV1[2] = 10000.0;
211 fVCov[4] = -1001.0;
212 }
213
214 AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
215
216 hbtEvent->SetPrimVertPos(vertex);
217 hbtEvent->SetPrimVertCov(fVCov);
218 }
219
220 if ((fInputType == kAOD) || (fInputType == kAODKine)) {
221 hbtEvent->SetRunNumber(fAODEvent->GetRunNumber());
222 hbtEvent->SetMagneticField(fAODEvent->GetMagneticField()*kilogauss);//to check if here is ok
223 hbtEvent->SetZDCN1Energy(fAODEvent->GetZDCN1Energy());
224 hbtEvent->SetZDCP1Energy(fAODEvent->GetZDCP1Energy());
225 hbtEvent->SetZDCN2Energy(fAODEvent->GetZDCN2Energy());
226 hbtEvent->SetZDCP2Energy(fAODEvent->GetZDCP2Energy());
227 hbtEvent->SetZDCEMEnergy(fAODEvent->GetZDCEMEnergy(0));
228 hbtEvent->SetZDCParticipants(0);
229 hbtEvent->SetTriggerMask(fAODEvent->GetTriggerMask());
230 hbtEvent->SetTriggerCluster(fAODEvent->GetTriggerCluster());
231
232 // Primary Vertex position
233 fAODEvent->GetPrimaryVertex()->GetPosition(fV1);
234
235 AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
236 hbtEvent->SetPrimVertPos(vertex);
237 }
238
239 if ((fInputType == kESDKine) || (fInputType == kAODKine)) {
240 Double_t tReactionPlane = 0;
241
242 AliGenHijingEventHeader *hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader);
243 if (!hdh) {
244 AliGenCocktailEventHeader *cdh = dynamic_cast<AliGenCocktailEventHeader *> (fGenHeader);
245 if (cdh) {
246 TList *tGenHeaders = cdh->GetHeaders();
247 for (int ihead = 0; ihead<tGenHeaders->GetEntries(); ihead++) {
248 hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader);
249 if (hdh) break;
250 }
251 }
252 }
253 if (hdh)
254 {
255 tReactionPlane = hdh->ReactionPlaneAngle();
256 cout << "Got reaction plane " << tReactionPlane << endl;
257 }
258
259 hbtEvent->SetReactionPlaneAngle(tReactionPlane);
260 }
261
262 //starting to reading tracks
263 int nofTracks=0; //number of reconstructed tracks in event
264 if ((fInputType == kESD) || (fInputType == kESDKine))
265 nofTracks = fESDEvent->GetNumberOfTracks();
266 else if ((fInputType == kAOD) || (fInputType == kAODKine))
267 nofTracks = fAODEvent->GetNumberOfTracks();
268
269 int realnofTracks=0;//number of track which we use ina analysis
270
271 TClonesArray *mcP = 0;
272 Int_t *motherids=0;
273
274 if (fInputType == kAODKine) {
275 // Attempt to access MC header
276 AliAODMCHeader *mcH;
277 mcH = (AliAODMCHeader *) fAODEvent->FindListObject(AliAODMCHeader::StdBranchName());
278 if (!mcH) {
279 cout << "AOD MC information requested, but no header found!" << endl;
280 }
281
282 mcP = (TClonesArray *) fAODEvent->FindListObject(AliAODMCParticle::StdBranchName());
283 if (!mcP) {
284 cout << "AOD MC information requested, but no particle array found!" << endl;
285 }
286
287 hbtEvent->SetReactionPlaneAngle(fAODEvent->GetHeader()->GetQTheta(0)/2.0);
288
289 if (mcP) {
290 motherids = new Int_t[((AliAODMCParticle *) mcP->At(mcP->GetEntries()-1))->GetLabel()];
291 for (int ip=0; ip<mcP->GetEntries(); ip++) motherids[ip] = 0;
292
293 // Read in mother ids
294 AliAODMCParticle *motherpart;
295 for (int ip=0; ip<mcP->GetEntries(); ip++) {
296 motherpart = (AliAODMCParticle *) mcP->At(ip);
297 if (motherpart->GetDaughter(0) > 0)
298 motherids[motherpart->GetDaughter(0)] = ip;
299 if (motherpart->GetDaughter(1) > 0)
300 motherids[motherpart->GetDaughter(1)] = ip;
301 }
302 }
303 }
304
305 if (fInputType == kESDKine) {
306 motherids = new Int_t[fStack->GetNtrack()];
307 for (int ip=0; ip<fStack->GetNtrack(); ip++) motherids[ip] = 0;
308
309 // Read in mother ids
310 TParticle *motherpart;
311 for (int ip=0; ip<fStack->GetNtrack(); ip++) {
312 motherpart = fStack->Particle(ip);
313 if (motherpart->GetDaughter(0) > 0)
314 motherids[motherpart->GetDaughter(0)] = ip;
315 if (motherpart->GetDaughter(1) > 0)
316 motherids[motherpart->GetDaughter(1)] = ip;
317 }
318 }
319
320 for (int i=0;i<nofTracks;i++)
321 {
322 // cout << "Reading track " << i << endl;
323 bool tGoodMomentum=true; //flaga to chcek if we can read momentum of this track
324
325 AliFemtoTrack* trackCopy = new AliFemtoTrack();
326
327 if ((fInputType == kESD) || (fInputType == kESDKine)) {
328
329 AliESDtrack *esdtrack=fESDEvent->GetTrack(i);//getting next track
330 if (fTrackCuts->AcceptTrack(esdtrack)) {
331
332 trackCopy->SetCharge((short)esdtrack->GetSign());
333
334 //in aliroot we have AliPID
335 //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
336 //we use only 5 first
337 double esdpid[5];
338 esdtrack->GetESDpid(esdpid);
339 trackCopy->SetPidProbElectron(esdpid[0]);
340 trackCopy->SetPidProbMuon(esdpid[1]);
341 trackCopy->SetPidProbPion(esdpid[2]);
342 trackCopy->SetPidProbKaon(esdpid[3]);
343 trackCopy->SetPidProbProton(esdpid[4]);
344
345 double pxyz[3];
346 double impact[2];
347 double covimpact[3];
348
349 // if (fUseTPCOnly) {
350 // if (!esdtrack->GetTPCInnerParam()) {
351 // cout << "No TPC inner param !" << endl;
352 // delete trackCopy;
353 // continue;
354 // }
355
356 // AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
357 // param->GetXYZ(rxyz);
358 // param->PropagateToDCA(fESDEvent->GetPrimaryVertexTPC(), (fESDEvent->GetMagneticField()), 10000, impact, covimpact);
359 // param->GetPxPyPz(pxyz);//reading noconstarined momentum
360
361 // if (fRotateToEventPlane) {
362 // double tPhi = TMath::ATan2(pxyz[1], pxyz[0]);
363 // double tRad = TMath::Hypot(pxyz[0], pxyz[1]);
364
365 // pxyz[0] = tRad*TMath::Cos(tPhi - tReactionPlane);
366 // pxyz[1] = tRad*TMath::Sin(tPhi - tReactionPlane);
367 // }
368
369 // AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
370 // if (v.mag() < 0.0001) {
371 // // cout << "Found 0 momentum ???? " << pxyz[0] << " " << pxyz[1] << " " << pxyz[2] << endl;
372 // delete trackCopy;
373 // continue;
374 // }
375
376 // trackCopy->SetP(v);//setting momentum
377 // trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
378
379 // const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
380 // const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
381 // //setting helix I do not if it is ok
382 // AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fESDEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
383 // trackCopy->SetHelix(helix);
384
385 // //some stuff which could be useful
386 // trackCopy->SetImpactD(impact[0]);
387 // trackCopy->SetImpactZ(impact[1]);
388 // trackCopy->SetCdd(covimpact[0]);
389 // trackCopy->SetCdz(covimpact[1]);
390 // trackCopy->SetCzz(covimpact[2]);
391 // trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact, covimpact));
392
393 // delete param;
394 // }
395 // else {
396 tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
397
398 AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
399 if (v.mag() < 0.0001) {
400
401 delete trackCopy;
402 continue;
403 }
404
405 trackCopy->SetP(v);//setting momentum
406 trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
407 const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
408 const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
409 //setting helix I do not if it is ok
410 AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fESDEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
411 trackCopy->SetHelix(helix);
412
413 //some stuff which could be useful
414 float imp[2];
415 float cim[3];
416 esdtrack->GetImpactParameters(imp,cim);
417
418 impact[0] = imp[0];
419 impact[1] = imp[1];
420 covimpact[0] = cim[0];
421 covimpact[1] = cim[1];
422 covimpact[2] = cim[2];
423
424 trackCopy->SetImpactD(impact[0]);
425 trackCopy->SetImpactZ(impact[1]);
426 trackCopy->SetCdd(covimpact[0]);
427 trackCopy->SetCdz(covimpact[1]);
428 trackCopy->SetCzz(covimpact[2]);
429 trackCopy->SetSigmaToVertex(AliESDtrackCuts::GetSigmaToVertex(esdtrack));
430
431 trackCopy->SetTrackId(esdtrack->GetID());
432 trackCopy->SetFlags(esdtrack->GetStatus());
433 trackCopy->SetLabel(esdtrack->GetLabel());
434
435 trackCopy->SetITSchi2(esdtrack->GetITSchi2());
436 trackCopy->SetITSncls(esdtrack->GetNcls(0));
437 trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());
438 trackCopy->SetTPCncls(esdtrack->GetTPCNcls());
439 trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());
440 trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class
441 trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma());
442
443 trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap());
444 trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap());
445
446 double xtpc[3];
447 esdtrack->GetInnerXYZ(xtpc);
448 xtpc[2] -= fV1[2];
449 trackCopy->SetNominalTPCEntrancePoint(xtpc);
450
451 esdtrack->GetOuterXYZ(xtpc);
452 xtpc[2] -= fV1[2];
453 trackCopy->SetNominalTPCExitPoint(xtpc);
454
455 int indexes[3];
456 for (int ik=0; ik<3; ik++) {
457 indexes[ik] = esdtrack->GetKinkIndex(ik);
458 }
459 trackCopy->SetKinkIndexes(indexes);
460
461 // Fill the hidden information with the simulated data
462 if (fInputType == kESDKine) {
463 if (TMath::Abs(esdtrack->GetLabel()) < fStack->GetNtrack()) {
464 TParticle *tPart = fStack->Particle(TMath::Abs(esdtrack->GetLabel()));
465
466 // Check the mother information
467
468 // Using the new way of storing the freeze-out information
469 // Final state particle is stored twice on the stack
470 // one copy (mother) is stored with original freeze-out information
471 // and is not tracked
472 // the other one (daughter) is stored with primary vertex position
473 // and is tracked
474
475 // Freeze-out coordinates
476 double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
477 fpx = tPart->Vx() - fV1[0];
478 fpy = tPart->Vy() - fV1[1];
479 fpz = tPart->Vz() - fV1[2];
480 fpt = tPart->T();
481
482 AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
483 tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
484
485 fpx *= 1e13;
486 fpy *= 1e13;
487 fpz *= 1e13;
488 fpt *= 1e13;
489
490 if (motherids[TMath::Abs(esdtrack->GetLabel())]>0) {
491 TParticle *mother = fStack->Particle(motherids[TMath::Abs(esdtrack->GetLabel())]);
492 // Check if this is the same particle stored twice on the stack
493 if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
494 // It is the same particle
495 // Read in the original freeze-out information
496 // and convert it from to [fm]
497
498 // EPOS style
499 fpx = mother->Vx()*1e13*0.197327;
500 fpy = mother->Vy()*1e13*0.197327;
501 fpz = mother->Vz()*1e13*0.197327;
502 fpt = mother->T() *1e13*0.197327;
503
504
505 // Therminator style
506 // fpx = mother->Vx()*1e13;
507 // fpy = mother->Vy()*1e13;
508 // fpz = mother->Vz()*1e13;
509 // fpt = mother->T() *1e13*3e10;
510
511 }
512 }
513
514 tInfo->SetPDGPid(tPart->GetPdgCode());
515
516 tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
517 Double_t mass2 = (tPart->Energy() *tPart->Energy() -
518 tPart->Px()*tPart->Px() -
519 tPart->Py()*tPart->Py() -
520 tPart->Pz()*tPart->Pz());
521 if (mass2>0.0)
522 tInfo->SetMass(TMath::Sqrt(mass2));
523 else
524 tInfo->SetMass(0.0);
525
526 tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
527 trackCopy->SetHiddenInfo(tInfo);
528 }
529 else {
530 AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
531 tInfo->SetMass(0.0);
532 double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
533 fpx = fV1[0]*1e13;
534 fpy = fV1[1]*1e13;
535 fpz = fV1[2]*1e13;
536 fpt = 0.0;
537 tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
538
539 tInfo->SetTrueMomentum(pxyz[0],pxyz[1],pxyz[2]);
540
541 trackCopy->SetHiddenInfo(tInfo);
542 }
543 }
544 // cout << "Got freeze-out " << fpx << " " << fpy << " " << fpz << " " << fpt << " " << mass2 << " " << tPart->GetPdgCode() << endl;
545 }
546 else
547 tGoodMomentum = false;
548 }
549
550 if ((fInputType == kAOD) || (fInputType == kAODKine)) {
551 // Read in the normal AliAODTracks
552 const AliAODTrack *aodtrack=fAODEvent->GetTrack(i); // getting the AODtrack directly
553
554 // if (!aodtrack->TestFilterBit(fFilterBit))
555 // continue;
556
557 CopyAODtoFemtoTrack(aodtrack, trackCopy);
558
559 if (mcP) {
560 // Fill the hidden information with the simulated data
561 // Int_t pLabel = aodtrack->GetLabel();
562 AliAODMCParticle *tPart = GetParticleWithLabel(mcP, (TMath::Abs(aodtrack->GetLabel())));
563
564 AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
565 double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
566 if (!tPart) {
567 fpx = fV1[0];
568 fpy = fV1[1];
569 fpz = fV1[2];
570 tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
571 tInfo->SetPDGPid(0);
572 tInfo->SetTrueMomentum(0.0, 0.0, 0.0);
573 tInfo->SetEmissionPoint(0.0, 0.0, 0.0, 0.0);
574 tInfo->SetMass(0);
575 }
576 else {
577 // Check the mother information
578
579 // Using the new way of storing the freeze-out information
580 // Final state particle is stored twice on the stack
581 // one copy (mother) is stored with original freeze-out information
582 // and is not tracked
583 // the other one (daughter) is stored with primary vertex position
584 // and is tracked
585
586 // Freeze-out coordinates
587 fpx = tPart->Xv() - fV1[0];
588 fpy = tPart->Yv() - fV1[1];
589 fpz = tPart->Zv() - fV1[2];
590 // fpt = tPart->T();
591
592 tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
593
594 fpx *= 1e13;
595 fpy *= 1e13;
596 fpz *= 1e13;
597 // fpt *= 1e13;
598
599 // cout << "Looking for mother ids " << endl;
600 if (motherids[TMath::Abs(aodtrack->GetLabel())]>0) {
601 // cout << "Got mother id" << endl;
602 AliAODMCParticle *mother = GetParticleWithLabel(mcP, motherids[TMath::Abs(aodtrack->GetLabel())]);
603 // Check if this is the same particle stored twice on the stack
604 if (mother) {
605 if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
606 // It is the same particle
607 // Read in the original freeze-out information
608 // and convert it from to [fm]
609
610 // EPOS style
611 // fpx = mother->Xv()*1e13*0.197327;
612 // fpy = mother->Yv()*1e13*0.197327;
613 // fpz = mother->Zv()*1e13*0.197327;
614 // fpt = mother->T() *1e13*0.197327*0.5;
615
616
617 // Therminator style
618 fpx = mother->Xv()*1e13;
619 fpy = mother->Yv()*1e13;
620 fpz = mother->Zv()*1e13;
621 // fpt = mother->T() *1e13*3e10;
622
623 }
624 }
625 }
626
627 // if (fRotateToEventPlane) {
628 // double tPhi = TMath::ATan2(fpy, fpx);
629 // double tRad = TMath::Hypot(fpx, fpy);
630
631 // fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
632 // fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
633 // }
634
635 tInfo->SetPDGPid(tPart->GetPdgCode());
636
637 // if (fRotateToEventPlane) {
638 // double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
639 // double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
640
641 // tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
642 // tRad*TMath::Sin(tPhi - tReactionPlane),
643 // tPart->Pz());
644 // }
645 // else
646 tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
647 Double_t mass2 = (tPart->E() *tPart->E() -
648 tPart->Px()*tPart->Px() -
649 tPart->Py()*tPart->Py() -
650 tPart->Pz()*tPart->Pz());
651 if (mass2>0.0)
652 tInfo->SetMass(TMath::Sqrt(mass2));
653 else
654 tInfo->SetMass(0.0);
655
656 tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
657 }
658 trackCopy->SetHiddenInfo(tInfo);
659 }
660
661 double pxyz[3];
662 aodtrack->PxPyPz(pxyz);//reading noconstarined momentum
663 const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
664 // Check the sanity of the tracks - reject zero momentum tracks
665 if (ktP.mag() == 0) {
666 delete trackCopy;
667 continue;
668 }
669
670
671 }
672
673
674 //decision if we want this track
675 //if we using diffrent labels we want that this label was use for first time
676 //if we use hidden info we want to have match between sim data and ESD
677 if (tGoodMomentum==true)
678 {
679 hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
680 realnofTracks++;//real number of tracks
681 // delete trackCopy;
682 }
683 else
684 {
685 delete trackCopy;
686 }
687
688 }
689
690 if (motherids)
691 delete motherids;
692
693 hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
694 fCurEvent++;
695 // cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
696 return hbtEvent;
697}
698//___________________
699void AliFemtoEventReaderStandard::SetESDSource(AliESDEvent *aESD)
700{
701 // The chain loads the ESD for us
702 // You must provide the address where it can be found
703 fESDEvent = aESD;
704}
705//___________________
706void AliFemtoEventReaderStandard::SetAODSource(AliAODEvent *aAOD)
707{
708 // The chain loads the ESD for us
709 // You must provide the address where it can be found
710 fAODEvent = aAOD;
711}
712//___________________
713void AliFemtoEventReaderStandard::SetStackSource(AliStack *aStack)
714{
715 // The chain loads the stack for us
716 // You must provide the address where it can be found
717 fStack = aStack;
718}
719void AliFemtoEventReaderStandard::SetInputType(AliFemtoInputType aInput)
720{
721 // Set the proper input type
722 fInputType = aInput;
723}
724//___________________
725void AliFemtoEventReaderStandard::SetGenEventHeader(AliGenEventHeader *aGenHeader)
726{
727 // The chain loads the generator event header for us
728 // You must provide the address where it can be found
729 fGenHeader = aGenHeader;
730}
731
732void AliFemtoEventReaderStandard::SetUsePhysicsSelection(const bool usephysics)
733{
734 fUsePhysicsSel = usephysics;
735 if (!fSelect) fSelect = new AliPhysicsSelection();
736}
737
738void AliFemtoEventReaderStandard::SetESDTrackCuts(AliESDtrackCuts *esdcuts)
739{
740 // Set external ESD track cuts
741 fTrackCuts = esdcuts;
742}
743
744void AliFemtoEventReaderStandard::CopyAODtoFemtoTrack(const AliAODTrack *tAodTrack,
745 AliFemtoTrack *tFemtoTrack)
746{
747 // Copy the track information from the AOD into the internal AliFemtoTrack
748 // If it exists, use the additional information from the PWG2 AOD
749
750 // Primary Vertex position
751 double fV1[3];
752 fAODEvent->GetPrimaryVertex()->GetPosition(fV1);
753
754 tFemtoTrack->SetCharge(tAodTrack->Charge());
755
756 //in aliroot we have AliPID
757 //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
758 //we use only 5 first
759
760 // AOD pid has 10 components
761 double aodpid[10];
762 tAodTrack->GetPID(aodpid);
763 tFemtoTrack->SetPidProbElectron(aodpid[0]);
764 tFemtoTrack->SetPidProbMuon(aodpid[1]);
765 tFemtoTrack->SetPidProbPion(aodpid[2]);
766 tFemtoTrack->SetPidProbKaon(aodpid[3]);
767 tFemtoTrack->SetPidProbProton(aodpid[4]);
768
769 double pxyz[3];
770 tAodTrack->PxPyPz(pxyz);//reading noconstrained momentum
771 AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
772 tFemtoTrack->SetP(v);//setting momentum
773 tFemtoTrack->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
774 const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
775 //setting track helix
776 const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
777 AliFmPhysicalHelixD helix(ktP,kOrigin,(double)(fAODEvent->GetMagneticField())*kilogauss,(double)(tFemtoTrack->Charge()));
778 tFemtoTrack->SetHelix(helix);
779
780 // Flags
781 tFemtoTrack->SetTrackId(tAodTrack->GetID());
782 tFemtoTrack->SetFlags(1);
783 tFemtoTrack->SetLabel(tAodTrack->GetLabel());
784
785 // Track quality information
786 float covmat[6];
787 tAodTrack->GetCovMatrix(covmat);
788 tFemtoTrack->SetImpactD(covmat[0]);
789 tFemtoTrack->SetImpactZ(covmat[2]);
790 tFemtoTrack->SetCdd(covmat[3]);
791 tFemtoTrack->SetCdz(covmat[4]);
792 tFemtoTrack->SetCzz(covmat[5]);
793 // This information is only available in the ESD
794 // We put in fake values or reasonable estimates
795 tFemtoTrack->SetITSchi2(tAodTrack->Chi2perNDF());
796 tFemtoTrack->SetITSncls(1);
797 tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF());
798 tFemtoTrack->SetTPCncls(1);
799 tFemtoTrack->SetTPCnclsF(1);
800 tFemtoTrack->SetTPCsignalN(1);
801 tFemtoTrack->SetTPCsignalS(1);
802
803 TBits tAllTrue;
804 TBits tAllFalse;
805 tAllTrue.ResetAllBits(kTRUE);
806 tAllFalse.ResetAllBits(kFALSE);
807
808
809 // If not use dummy values
810 tFemtoTrack->SetTPCClusterMap(tAllTrue);
811 tFemtoTrack->SetTPCSharedMap(tAllFalse);
812
813 double xtpc[3] = {0,0,0};
814 tFemtoTrack->SetNominalTPCEntrancePoint(xtpc);
815 tFemtoTrack->SetNominalTPCExitPoint(xtpc);
816
817 int indexes[3];
818 for (int ik=0; ik<3; ik++) {
819 indexes[ik] = 0;
820 }
821 tFemtoTrack->SetKinkIndexes(indexes);
822}
823
824AliAODMCParticle* AliFemtoEventReaderStandard::GetParticleWithLabel(TClonesArray *mcP, Int_t aLabel)
825{
826 if (aLabel < 0) return 0;
827 AliAODMCParticle *aodP;
828 Int_t posstack = 0;
829 if (aLabel > mcP->GetEntries())
830 posstack = mcP->GetEntries();
831 else
832 posstack = aLabel;
833
834 aodP = (AliAODMCParticle *) mcP->At(posstack);
835 if (aodP->GetLabel() > posstack) {
836 do {
837 aodP = (AliAODMCParticle *) mcP->At(posstack);
838 if (aodP->GetLabel() == aLabel) return aodP;
839 posstack--;
840 }
841 while (posstack > 0);
842 }
843 else {
844 do {
845 aodP = (AliAODMCParticle *) mcP->At(posstack);
846 if (aodP->GetLabel() == aLabel) return aodP;
847 posstack++;
848 }
849 while (posstack < mcP->GetEntries());
850 }
851
852 return 0;
853}