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