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