Correct DCA propagation
[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemto / AliFemtoEventReaderAOD.cxx
CommitLineData
de568db1 1////////////////////////////////////////////////////////////////////////////////
2// //
3// AliFemtoEventReaderAOD - the reader class for the Alice AOD //
4// Reads in AOD information and converts it into internal AliFemtoEvent //
5// Authors: Marek Chojnacki mchojnacki@knf.pw.edu.pl //
6// Adam Kisiel kisiel@mps.ohio-state.edu //
7// //
8////////////////////////////////////////////////////////////////////////////////
9
10#include "AliFemtoEventReaderAOD.h"
11
12#include "TFile.h"
13#include "TTree.h"
14#include "AliAODEvent.h"
15#include "AliAODTrack.h"
16#include "AliAODVertex.h"
683877d2 17#include "AliAODMCHeader.h"
4637c3d5 18#include "AliESDtrack.h"
de568db1 19
20#include "AliFmPhysicalHelixD.h"
21#include "AliFmThreeVectorF.h"
22
23#include "SystemOfUnits.h"
24
25#include "AliFemtoEvent.h"
26#include "AliFemtoModelHiddenInfo.h"
683877d2 27#include "AliFemtoModelGlobalHiddenInfo.h"
1309aa0b 28#include "AliPID.h"
29
30#include "AliAODpidUtil.h"
de568db1 31
32ClassImp(AliFemtoEventReaderAOD)
33
34#if !(ST_NO_NAMESPACES)
35 using namespace units;
36#endif
37
38using namespace std;
39//____________________________
40//constructor with 0 parameters , look at default settings
41AliFemtoEventReaderAOD::AliFemtoEventReaderAOD():
de568db1 42 fNumberofEvent(0),
43 fCurEvent(0),
de568db1 44 fEvent(0x0),
45 fAllTrue(160),
46 fAllFalse(160),
47 fFilterBit(0),
fcda1d4e 48 fPWG2AODTracks(0x0),
683877d2 49 fReadMC(0),
948862a7 50 fUsePreCent(0),
4637c3d5 51 fAODpidUtil(0),
fcda1d4e 52 fInputFile(" "),
53 fFileName(" "),
54 fTree(0x0),
4637c3d5 55 fAodFile(0x0)
de568db1 56{
57 // default constructor
58 fAllTrue.ResetAllBits(kTRUE);
59 fAllFalse.ResetAllBits(kFALSE);
948862a7 60 fCentRange[0] = 0;
61 fCentRange[1] = 1000;
de568db1 62}
63
64AliFemtoEventReaderAOD::AliFemtoEventReaderAOD(const AliFemtoEventReaderAOD &aReader) :
fcda1d4e 65 AliFemtoEventReader(),
de568db1 66 fNumberofEvent(0),
67 fCurEvent(0),
de568db1 68 fEvent(0x0),
69 fAllTrue(160),
70 fAllFalse(160),
71 fFilterBit(0),
fcda1d4e 72 fPWG2AODTracks(0x0),
683877d2 73 fReadMC(0),
948862a7 74 fUsePreCent(0),
4637c3d5 75 fAODpidUtil(0),
fcda1d4e 76 fInputFile(" "),
77 fFileName(" "),
78 fTree(0x0),
4637c3d5 79 fAodFile(0x0)
de568db1 80{
81 // copy constructor
82 fInputFile = aReader.fInputFile;
83 fFileName = aReader.fFileName;
84 fNumberofEvent = aReader.fNumberofEvent;
85 fCurEvent = aReader.fCurEvent;
86 fEvent = new AliAODEvent();
87 fAodFile = new TFile(aReader.fAodFile->GetName());
88 fAllTrue.ResetAllBits(kTRUE);
89 fAllFalse.ResetAllBits(kFALSE);
90 fFilterBit = aReader.fFilterBit;
91 fPWG2AODTracks = aReader.fPWG2AODTracks;
1309aa0b 92 fAODpidUtil = aReader.fAODpidUtil;
948862a7 93 fCentRange[0] = aReader.fCentRange[0];
94 fCentRange[1] = aReader.fCentRange[1];
de568db1 95}
96//__________________
97//Destructor
98AliFemtoEventReaderAOD::~AliFemtoEventReaderAOD()
99{
100 // destructor
101 delete fTree;
102 delete fEvent;
103 delete fAodFile;
104 if (fPWG2AODTracks) {
105 fPWG2AODTracks->Delete();
106 delete fPWG2AODTracks;
107 }
108}
109
110//__________________
111AliFemtoEventReaderAOD& AliFemtoEventReaderAOD::operator=(const AliFemtoEventReaderAOD& aReader)
112{
113 // assignment operator
114 if (this == &aReader)
115 return *this;
116
117 fInputFile = aReader.fInputFile;
118 fFileName = aReader.fFileName;
119 fNumberofEvent = aReader.fNumberofEvent;
120 fCurEvent = aReader.fCurEvent;
121 if (fTree) delete fTree;
122 if (fEvent) delete fEvent;
123 fEvent = new AliAODEvent();
124 if (fAodFile) delete fAodFile;
125 fAodFile = new TFile(aReader.fAodFile->GetName());
126 fAllTrue.ResetAllBits(kTRUE);
127 fAllFalse.ResetAllBits(kFALSE);
128 fFilterBit = aReader.fFilterBit;
129 fPWG2AODTracks = aReader.fPWG2AODTracks;
1309aa0b 130 fAODpidUtil = aReader.fAODpidUtil;
948862a7 131 fCentRange[0] = aReader.fCentRange[0];
132 fCentRange[1] = aReader.fCentRange[1];
de568db1 133
134 return *this;
135}
136//__________________
137AliFemtoString AliFemtoEventReaderAOD::Report()
138{
139 // create reader report
140 AliFemtoString temp = "\n This is the AliFemtoEventReaderAOD\n";
141 return temp;
142}
143
144//__________________
145void AliFemtoEventReaderAOD::SetInputFile(const char* inputFile)
146{
147 //setting the name of file where names of AOD file are written
148 //it takes only this files which have good trees
149 char buffer[256];
150 fInputFile=string(inputFile);
de568db1 151 ifstream infile(inputFile);
152
153 fTree = new TChain("aodTree");
154
155 if(infile.good()==true)
156 {
157 //checking if all give files have good tree inside
158 while (infile.eof()==false)
159 {
160 infile.getline(buffer,256);
161 TFile *aodFile=TFile::Open(buffer,"READ");
162 if (aodFile!=0x0)
163 {
164 TTree* tree = (TTree*) aodFile->Get("aodTree");
165 if (tree!=0x0)
166 {
fd7bd835 167 // cout<<"putting file "<<string(buffer)<<" into analysis"<<endl;
de568db1 168 fTree->AddFile(buffer);
169 delete tree;
170 }
171 aodFile->Close();
172 }
173 delete aodFile;
174 }
175 }
176}
177
178AliFemtoEvent* AliFemtoEventReaderAOD::ReturnHbtEvent()
179{
180 // read in a next hbt event from the chain
181 // convert it to AliFemtoEvent and return
182 // for further analysis
183 AliFemtoEvent *hbtEvent = 0;
fd7bd835 184 cout<<"reader"<<endl;
de568db1 185 if (fCurEvent==fNumberofEvent)//open next file
186 {
187 if(fNumberofEvent==0)
188 {
189 fEvent=new AliAODEvent();
190 fEvent->ReadFromTree(fTree);
191
192 // Check for the existence of the additional information
193 fPWG2AODTracks = (TClonesArray *) fEvent->GetList()->FindObject("pwg2aodtracks");
194
195 if (fPWG2AODTracks) {
196 cout << "Found additional PWG2 specific information in the AOD!" << endl;
197 cout << "Reading only tracks with the additional information" << endl;
198 }
199
200 fNumberofEvent=fTree->GetEntries();
683877d2 201 // cout<<"Number of Entries in file "<<fNumberofEvent<<endl;
de568db1 202 fCurEvent=0;
203 }
204 else //no more data to read
205 {
206 cout<<"no more files "<<hbtEvent<<endl;
207 fReaderStatus=1;
208 return hbtEvent;
209 }
210 }
211
212 cout<<"starting to read event "<<fCurEvent<<endl;
213 fTree->GetEvent(fCurEvent);//getting next event
683877d2 214 // cout << "Read event " << fEvent << " from file " << fTree << endl;
de568db1 215
216 hbtEvent = new AliFemtoEvent;
217
218 CopyAODtoFemtoEvent(hbtEvent);
fd7bd835 219 fCurEvent++;
de568db1 220
de568db1 221
222 return hbtEvent;
223}
224
225void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
226{
1309aa0b 227
de568db1 228 // A function that reads in the AOD event
229 // and transfers the neccessary information into
230 // the internal AliFemtoEvent
231
232 // setting global event characteristics
233 tEvent->SetRunNumber(fEvent->GetRunNumber());
234 tEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
235 tEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
236 tEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
237 tEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
238 tEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
239 tEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy(0));
b976850d 240 tEvent->SetZDCParticipants(0);
de568db1 241 tEvent->SetTriggerMask(fEvent->GetTriggerMask());
242 tEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
fd7bd835 243
683877d2 244 // Attempt to access MC header
245 AliAODMCHeader *mcH;
d9ae2120 246 TClonesArray *mcP=0;
683877d2 247 if (fReadMC) {
248 mcH = (AliAODMCHeader *) fEvent->FindListObject(AliAODMCHeader::StdBranchName());
249 if (!mcH) {
250 cout << "AOD MC information requested, but no header found!" << endl;
251 }
252
253 mcP = (TClonesArray *) fEvent->FindListObject(AliAODMCParticle::StdBranchName());
254 if (!mcP) {
255 cout << "AOD MC information requested, but no particle array found!" << endl;
256 }
257 }
258
fee52126 259 tEvent->SetReactionPlaneAngle(fEvent->GetHeader()->GetQTheta(0)/2.0);
260
d9ae2120 261 Int_t *motherids=0;
683877d2 262 if (mcP) {
263 motherids = new Int_t[((AliAODMCParticle *) mcP->At(mcP->GetEntries()-1))->GetLabel()];
264 for (int ip=0; ip<mcP->GetEntries(); ip++) motherids[ip] = 0;
265
266 // Read in mother ids
267 AliAODMCParticle *motherpart;
268 for (int ip=0; ip<mcP->GetEntries(); ip++) {
269 motherpart = (AliAODMCParticle *) mcP->At(ip);
270 if (motherpart->GetDaughter(0) > 0)
271 motherids[motherpart->GetDaughter(0)] = ip;
272 if (motherpart->GetDaughter(1) > 0)
273 motherids[motherpart->GetDaughter(1)] = ip;
274 }
275 }
276
de568db1 277 // Primary Vertex position
278 double fV1[3];
279 fEvent->GetPrimaryVertex()->GetPosition(fV1);
280
281 AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
282 tEvent->SetPrimVertPos(vertex);
283
284 //starting to reading tracks
285 int nofTracks=0; //number of reconstructed tracks in event
286
287 // Check to see whether the additional info exists
288 if (fPWG2AODTracks)
289 nofTracks=fPWG2AODTracks->GetEntries();
290 else
291 nofTracks=fEvent->GetNumberOfTracks();
1309aa0b 292 cout<<"nofTracks: "<<nofTracks<<endl;
de568db1 293
948862a7 294 AliCentrality *cent = fEvent->GetCentrality();
fd7bd835 295
948862a7 296 if (cent && fUsePreCent) {
297 if ((cent->GetCentralityPercentile("V0M")*10 < fCentRange[0]) ||
298 (cent->GetCentralityPercentile("V0M")*10 > fCentRange[1]))
299 {
300 cout << "Centrality " << cent->GetCentralityPercentile("V0M") << " outside of preselection range " << fCentRange[0] << " - " << fCentRange[1] << endl;
1309aa0b 301
948862a7 302 return;
303 }
304 }
305
de568db1 306 int realnofTracks=0; // number of track which we use in a analysis
707be29d 307 int tracksPrim=0;
de568db1 308
63073cb3 309 int labels[20000];
310 for (int il=0; il<20000; il++) labels[il] = -1;
311
312 // looking for global tracks and saving their numbers to copy from them PID information to TPC-only tracks in the main loop over tracks
313 for (int i=0;i<nofTracks;i++) {
314 const AliAODTrack *aodtrack=fEvent->GetTrack(i);
315 if (!aodtrack->TestFilterBit(fFilterBit)) {
316 labels[aodtrack->GetID()] = i;
317 }
318 }
319
a22737c8 320 int tNormMult = 0;
de568db1 321 for (int i=0;i<nofTracks;i++)
322 {
323 AliFemtoTrack* trackCopy = new AliFemtoTrack();
324
325 if (fPWG2AODTracks) {
326 // Read tracks from the additional pwg2 specific AOD part
327 // if they exist
328 // Note that in that case all the AOD tracks without the
329 // additional information will be ignored !
330 AliPWG2AODTrack *pwg2aodtrack = (AliPWG2AODTrack *) fPWG2AODTracks->At(i);
331
332 // Getting the AOD track through the ref of the additional info
333 AliAODTrack *aodtrack = pwg2aodtrack->GetRefAODTrack();
c4b81034 334 if (!aodtrack->TestFilterBit(fFilterBit)) {
335 delete trackCopy;
de568db1 336 continue;
c4b81034 337 }
338
1309aa0b 339
a22737c8 340 if (aodtrack->IsOn(AliESDtrack::kTPCrefit))
341 if (aodtrack->Chi2perNDF() < 6.0)
342 if (aodtrack->Eta() < 0.9)
343 tNormMult++;
344
1309aa0b 345
de568db1 346 CopyAODtoFemtoTrack(aodtrack, trackCopy, pwg2aodtrack);
347
683877d2 348 if (mcP) {
349 // Fill the hidden information with the simulated data
d9ae2120 350 // Int_t pLabel = aodtrack->GetLabel();
683877d2 351 AliAODMCParticle *tPart = GetParticleWithLabel(mcP, (TMath::Abs(aodtrack->GetLabel())));
352
353 // Check the mother information
354
355 // Using the new way of storing the freeze-out information
356 // Final state particle is stored twice on the stack
357 // one copy (mother) is stored with original freeze-out information
358 // and is not tracked
359 // the other one (daughter) is stored with primary vertex position
360 // and is tracked
361
362 // Freeze-out coordinates
363 double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
364 fpx = tPart->Xv() - fV1[0];
365 fpy = tPart->Yv() - fV1[1];
366 fpz = tPart->Zv() - fV1[2];
367 fpt = tPart->T();
368
369 AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
370 tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
371
372 fpx *= 1e13;
373 fpy *= 1e13;
374 fpz *= 1e13;
375 fpt *= 1e13;
376
377 // cout << "Looking for mother ids " << endl;
378 if (motherids[TMath::Abs(aodtrack->GetLabel())]>0) {
379 // cout << "Got mother id" << endl;
380 AliAODMCParticle *mother = GetParticleWithLabel(mcP, motherids[TMath::Abs(aodtrack->GetLabel())]);
381 // Check if this is the same particle stored twice on the stack
382 if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
383 // It is the same particle
384 // Read in the original freeze-out information
385 // and convert it from to [fm]
386
387 // EPOS style
388 // fpx = mother->Xv()*1e13*0.197327;
389 // fpy = mother->Yv()*1e13*0.197327;
390 // fpz = mother->Zv()*1e13*0.197327;
391 // fpt = mother->T() *1e13*0.197327*0.5;
392
393
394 // Therminator style
395 fpx = mother->Xv()*1e13;
396 fpy = mother->Yv()*1e13;
397 fpz = mother->Zv()*1e13;
398 fpt = mother->T() *1e13*3e10;
399
400 }
401 }
402
fd7bd835 403 // if (fRotateToEventPlane) {
404 // double tPhi = TMath::ATan2(fpy, fpx);
405 // double tRad = TMath::Hypot(fpx, fpy);
683877d2 406
fd7bd835 407 // fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
408 // fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
409 // }
683877d2 410
411 tInfo->SetPDGPid(tPart->GetPdgCode());
412
fd7bd835 413 // if (fRotateToEventPlane) {
414 // double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
415 // double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
683877d2 416
fd7bd835 417 // tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
418 // tRad*TMath::Sin(tPhi - tReactionPlane),
419 // tPart->Pz());
420 // }
421 // else
683877d2 422 tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
423 Double_t mass2 = (tPart->E() *tPart->E() -
424 tPart->Px()*tPart->Px() -
425 tPart->Py()*tPart->Py() -
426 tPart->Pz()*tPart->Pz());
427 if (mass2>0.0)
428 tInfo->SetMass(TMath::Sqrt(mass2));
429 else
430 tInfo->SetMass(0.0);
431
432 tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
433 trackCopy->SetHiddenInfo(tInfo);
434
435 }
436
de568db1 437 double pxyz[3];
438 aodtrack->PxPyPz(pxyz);//reading noconstarined momentum
439 const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
440 // Check the sanity of the tracks - reject zero momentum tracks
69c1c8ff 441 if (ktP.Mag() == 0) {
de568db1 442 delete trackCopy;
443 continue;
444 }
445 }
446 else {
447 // No additional information exists
448 // Read in the normal AliAODTracks
fd7bd835 449
450 // const AliAODTrack *aodtrack=fEvent->GetTrack(i); // getting the AODtrack directly
451 AliAODTrack *aodtrack=fEvent->GetTrack(i); // getting the AODtrack directly
707be29d 452
453 if (aodtrack->IsPrimaryCandidate()) tracksPrim++;
de568db1 454
63073cb3 455 if (!aodtrack->TestFilterBit(fFilterBit)) {
948862a7 456 delete trackCopy;
63073cb3 457 continue;
948862a7 458 }
1309aa0b 459
de568db1 460 CopyAODtoFemtoTrack(aodtrack, trackCopy, 0);
63073cb3 461
462 // copying PID information from the correspondent track
fd7bd835 463 // const AliAODTrack *aodtrackpid = fEvent->GetTrack(labels[-1-fEvent->GetTrack(i)->GetID()]);
464 AliAODTrack *aodtrackpid = fEvent->GetTrack(labels[-1-fEvent->GetTrack(i)->GetID()]);
6cc142b2 465 CopyPIDtoFemtoTrack(aodtrackpid, trackCopy);
63073cb3 466
683877d2 467 if (mcP) {
468 // Fill the hidden information with the simulated data
d9ae2120 469 // Int_t pLabel = aodtrack->GetLabel();
683877d2 470 AliAODMCParticle *tPart = GetParticleWithLabel(mcP, (TMath::Abs(aodtrack->GetLabel())));
471
472 AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
473 double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
474 if (!tPart) {
475 fpx = fV1[0];
476 fpy = fV1[1];
477 fpz = fV1[2];
478 tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
479 tInfo->SetPDGPid(0);
480 tInfo->SetTrueMomentum(0.0, 0.0, 0.0);
481 tInfo->SetEmissionPoint(0.0, 0.0, 0.0, 0.0);
482 tInfo->SetMass(0);
483 }
484 else {
485 // Check the mother information
486
487 // Using the new way of storing the freeze-out information
488 // Final state particle is stored twice on the stack
489 // one copy (mother) is stored with original freeze-out information
490 // and is not tracked
491 // the other one (daughter) is stored with primary vertex position
492 // and is tracked
493
494 // Freeze-out coordinates
495 fpx = tPart->Xv() - fV1[0];
496 fpy = tPart->Yv() - fV1[1];
497 fpz = tPart->Zv() - fV1[2];
498 // fpt = tPart->T();
499
500 tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
501
502 fpx *= 1e13;
503 fpy *= 1e13;
504 fpz *= 1e13;
505 // fpt *= 1e13;
506
507 // cout << "Looking for mother ids " << endl;
508 if (motherids[TMath::Abs(aodtrack->GetLabel())]>0) {
509 // cout << "Got mother id" << endl;
510 AliAODMCParticle *mother = GetParticleWithLabel(mcP, motherids[TMath::Abs(aodtrack->GetLabel())]);
511 // Check if this is the same particle stored twice on the stack
512 if (mother) {
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->Xv()*1e13*0.197327;
520 // fpy = mother->Yv()*1e13*0.197327;
521 // fpz = mother->Zv()*1e13*0.197327;
522 // fpt = mother->T() *1e13*0.197327*0.5;
523
524
525 // Therminator style
526 fpx = mother->Xv()*1e13;
527 fpy = mother->Yv()*1e13;
528 fpz = mother->Zv()*1e13;
529 // fpt = mother->T() *1e13*3e10;
530
531 }
532 }
533 }
534
535 // if (fRotateToEventPlane) {
536 // double tPhi = TMath::ATan2(fpy, fpx);
537 // double tRad = TMath::Hypot(fpx, fpy);
538
539 // fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
540 // fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
541 // }
542
543 tInfo->SetPDGPid(tPart->GetPdgCode());
544
545 // if (fRotateToEventPlane) {
546 // double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
547 // double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
548
549 // tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
550 // tRad*TMath::Sin(tPhi - tReactionPlane),
551 // tPart->Pz());
552 // }
553 // else
554 tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
555 Double_t mass2 = (tPart->E() *tPart->E() -
556 tPart->Px()*tPart->Px() -
557 tPart->Py()*tPart->Py() -
558 tPart->Pz()*tPart->Pz());
559 if (mass2>0.0)
560 tInfo->SetMass(TMath::Sqrt(mass2));
561 else
562 tInfo->SetMass(0.0);
563
564 tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
565 }
566 trackCopy->SetHiddenInfo(tInfo);
567 }
568
de568db1 569 double pxyz[3];
570 aodtrack->PxPyPz(pxyz);//reading noconstarined momentum
571 const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
572 // Check the sanity of the tracks - reject zero momentum tracks
69c1c8ff 573 if (ktP.Mag() == 0) {
de568db1 574 delete trackCopy;
575 continue;
576 }
577 }
578
579
580 tEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
581 realnofTracks++;//real number of tracks
582 }
583
584 tEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
707be29d 585 tEvent->SetNormalizedMult(tracksPrim);
586
4637c3d5 587 // AliCentrality *cent = fEvent->GetCentrality();
948862a7 588 if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("V0M")));
63073cb3 589 // if (cent) tEvent->SetNormalizedMult((int) cent->GetCentralityPercentile("V0M"));
de568db1 590
fe77f653 591 if (cent) {
592 tEvent->SetCentralityV0(cent->GetCentralityPercentile("V0M"));
593 // tEvent->SetCentralityFMD(cent->GetCentralityPercentile("FMD"));
594 tEvent->SetCentralitySPD1(cent->GetCentralityPercentile("CL1"));
595 // tEvent->SetCentralityTrk(cent->GetCentralityPercentile("TRK"));
596 }
597
598
4ce766eb 599 if (mcP) delete [] motherids;
683877d2 600
de568db1 601 cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
602}
603
fd7bd835 604void AliFemtoEventReaderAOD::CopyAODtoFemtoTrack(AliAODTrack *tAodTrack,
de568db1 605 AliFemtoTrack *tFemtoTrack,
606 AliPWG2AODTrack *tPWG2AODTrack)
607{
608 // Copy the track information from the AOD into the internal AliFemtoTrack
609 // If it exists, use the additional information from the PWG2 AOD
610
611 // Primary Vertex position
612 double fV1[3];
613 fEvent->GetPrimaryVertex()->GetPosition(fV1);
fd7bd835 614 // fEvent->GetPrimaryVertex()->GetXYZ(fV1);
de568db1 615
616 tFemtoTrack->SetCharge(tAodTrack->Charge());
617
de568db1 618 double pxyz[3];
619 tAodTrack->PxPyPz(pxyz);//reading noconstrained momentum
620 AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
621 tFemtoTrack->SetP(v);//setting momentum
622 tFemtoTrack->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
623 const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
624 //setting track helix
625 const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
626 AliFmPhysicalHelixD helix(ktP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(tFemtoTrack->Charge()));
627 tFemtoTrack->SetHelix(helix);
628
629 // Flags
630 tFemtoTrack->SetTrackId(tAodTrack->GetID());
707be29d 631 tFemtoTrack->SetFlags(tAodTrack->GetFlags());
de568db1 632 tFemtoTrack->SetLabel(tAodTrack->GetLabel());
633
634 // Track quality information
635 float covmat[6];
707be29d 636 tAodTrack->GetCovMatrix(covmat);
637
fd7bd835 638 double impact[2];
639 double covimpact[3];
640 tAodTrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact);
641
63073cb3 642 // if (TMath::Abs(tAodTrack->Xv()) > 0.00000000001)
643 // tFemtoTrack->SetImpactD(TMath::Hypot(tAodTrack->Xv(), tAodTrack->Yv())*(tAodTrack->Xv()/TMath::Abs(tAodTrack->Xv())));
644 // else
645 // tFemtoTrack->SetImpactD(0.0);
646 // tFemtoTrack->SetImpactD(tAodTrack->DCA());
707be29d 647
63073cb3 648 // tFemtoTrack->SetImpactZ(tAodTrack->ZAtDCA());
fd7bd835 649
650
651 // tFemtoTrack->SetImpactD(TMath::Hypot(tAodTrack->Xv() - fV1[0], tAodTrack->Yv() - fV1[1]));
652 // tFemtoTrack->SetImpactZ(tAodTrack->Zv() - fV1[2]);
653 tFemtoTrack->SetImpactD(impact[0]);
654 tFemtoTrack->SetImpactZ(impact[1]);
655
656// cout
657// // << "dca" << TMath::Hypot(tAodTrack->Xv() - fV1[0], tAodTrack->Yv() - fV1[1])
658// // << "xv - fv10 = "<< tAodTrack->Xv() - fV1[0]
659// // << tAodTrack->Yv() - fV1[1]
660// << "xv = " << tAodTrack->Xv() << endl
661// << "fv1[0] = " << fV1[0] << endl
662// << "yv = " << tAodTrack->Yv() << endl
663// << "fv1[1] = " << fV1[1] << endl
664// << "zv = " << tAodTrack->Zv() << endl
665// << "fv1[2] = " << fV1[2] << endl
666// << "impact[0] = " << impact[0] << endl
667// << "impact[1] = " << impact[1] << endl
668// << endl << endl ;
707be29d 669
707be29d 670 tFemtoTrack->SetCdd(covmat[0]);
671 tFemtoTrack->SetCdz(covmat[1]);
672 tFemtoTrack->SetCzz(covmat[2]);
de568db1 673 tFemtoTrack->SetITSchi2(tAodTrack->Chi2perNDF());
707be29d 674 tFemtoTrack->SetITSncls(tAodTrack->GetITSNcls());
de568db1 675 tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF());
707be29d 676 tFemtoTrack->SetTPCncls(tAodTrack->GetTPCNcls());
677 tFemtoTrack->SetTPCnclsF(tAodTrack->GetTPCNcls());
683877d2 678 tFemtoTrack->SetTPCsignalN(1);
679 tFemtoTrack->SetTPCsignalS(1);
1309aa0b 680 tFemtoTrack->SetTPCsignal(tAodTrack->GetTPCsignal());
de568db1 681
682 if (tPWG2AODTrack) {
683 // Copy the PWG2 specific information if it exists
684 tFemtoTrack->SetTPCClusterMap(tPWG2AODTrack->GetTPCClusterMap());
685 tFemtoTrack->SetTPCSharedMap(tPWG2AODTrack->GetTPCSharedMap());
686
687 double xtpc[3] = {0,0,0};
688 tPWG2AODTrack->GetTPCNominalEntrancePoint(xtpc);
689 tFemtoTrack->SetNominalTPCEntrancePoint(xtpc);
690 tPWG2AODTrack->GetTPCNominalExitPoint(xtpc);
691 tFemtoTrack->SetNominalTPCExitPoint(xtpc);
692 }
693 else {
694 // If not use dummy values
707be29d 695 tFemtoTrack->SetTPCClusterMap(tAodTrack->GetTPCClusterMap());
696 tFemtoTrack->SetTPCSharedMap(tAodTrack->GetTPCSharedMap());
de568db1 697
698 double xtpc[3] = {0,0,0};
699 tFemtoTrack->SetNominalTPCEntrancePoint(xtpc);
700 tFemtoTrack->SetNominalTPCExitPoint(xtpc);
701 }
702
707be29d 703 // cout << "Track has " << TMath::Hypot(tAodTrack->Xv(), tAodTrack->Yv()) << " " << tAodTrack->Zv() << " " << tAodTrack->GetTPCNcls() << endl;
704
705
de568db1 706 int indexes[3];
707 for (int ik=0; ik<3; ik++) {
708 indexes[ik] = 0;
709 }
710 tFemtoTrack->SetKinkIndexes(indexes);
711}
712
713void AliFemtoEventReaderAOD::SetFilterBit(UInt_t ibit)
714{
715 fFilterBit = (1 << (ibit));
716}
717
683877d2 718void AliFemtoEventReaderAOD::SetReadMC(unsigned char a)
719{
720 fReadMC = a;
721}
722
723AliAODMCParticle* AliFemtoEventReaderAOD::GetParticleWithLabel(TClonesArray *mcP, Int_t aLabel)
724{
725 if (aLabel < 0) return 0;
726 AliAODMCParticle *aodP;
727 Int_t posstack = 0;
728 if (aLabel > mcP->GetEntries())
729 posstack = mcP->GetEntries();
730 else
731 posstack = aLabel;
732
733 aodP = (AliAODMCParticle *) mcP->At(posstack);
734 if (aodP->GetLabel() > posstack) {
735 do {
736 aodP = (AliAODMCParticle *) mcP->At(posstack);
737 if (aodP->GetLabel() == aLabel) return aodP;
738 posstack--;
739 }
740 while (posstack > 0);
741 }
742 else {
743 do {
744 aodP = (AliAODMCParticle *) mcP->At(posstack);
745 if (aodP->GetLabel() == aLabel) return aodP;
746 posstack++;
747 }
748 while (posstack < mcP->GetEntries());
749 }
750
751 return 0;
752}
de568db1 753
fd7bd835 754void AliFemtoEventReaderAOD::CopyPIDtoFemtoTrack(AliAODTrack *tAodTrack,
6cc142b2 755 AliFemtoTrack *tFemtoTrack)
63073cb3 756{
757 double aodpid[10];
758 tAodTrack->GetPID(aodpid);
759 tFemtoTrack->SetPidProbElectron(aodpid[0]);
760 tFemtoTrack->SetPidProbMuon(aodpid[1]);
761 tFemtoTrack->SetPidProbPion(aodpid[2]);
762 tFemtoTrack->SetPidProbKaon(aodpid[3]);
763 tFemtoTrack->SetPidProbProton(aodpid[4]);
764
765 aodpid[0] = -100000.0;
766 aodpid[1] = -100000.0;
767 aodpid[2] = -100000.0;
768 aodpid[3] = -100000.0;
769 aodpid[4] = -100000.0;
770
771 double tTOF = 0.0;
772
773 if (tAodTrack->GetStatus() & AliESDtrack::kTOFpid) { //AliESDtrack::kTOFpid=0x8000
774 tTOF = tAodTrack->GetTOFsignal();
775 tAodTrack->GetIntegratedTimes(aodpid);
776 }
777
778 tFemtoTrack->SetTofExpectedTimes(tTOF-aodpid[2], tTOF-aodpid[3], tTOF-aodpid[4]);
779
780 ////// TPC ////////////////////////////////////////////
781
782 float nsigmaTPCK=-1000.;
783 float nsigmaTPCPi=-1000.;
784 float nsigmaTPCP=-1000.;
785
786 // cout<<"in reader fESDpid"<<fESDpid<<endl;
787
788 if (tAodTrack->IsOn(AliESDtrack::kTPCpid)){ //AliESDtrack::kTPCpid=0x0080
789 nsigmaTPCK = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kKaon);
790 nsigmaTPCPi = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kPion);
791 nsigmaTPCP = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kProton);
792 }
793
794 tFemtoTrack->SetNSigmaTPCPi(nsigmaTPCPi);
795 tFemtoTrack->SetNSigmaTPCK(nsigmaTPCK);
796 tFemtoTrack->SetNSigmaTPCP(nsigmaTPCP);
797
798 tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF());
799 tFemtoTrack->SetTPCncls(tAodTrack->GetTPCNcls());
800 tFemtoTrack->SetTPCnclsF(tAodTrack->GetTPCNcls());
801
802 tFemtoTrack->SetTPCsignalN(1);
803 tFemtoTrack->SetTPCsignalS(1);
804 tFemtoTrack->SetTPCsignal(tAodTrack->GetTPCsignal());
805
806 ///////TOF//////////////////////
807
808 float vp=-1000.;
809 float nsigmaTOFPi=-1000.;
810 float nsigmaTOFK=-1000.;
811 float nsigmaTOFP=-1000.;
812
813 if ((tAodTrack->GetStatus() & AliESDtrack::kTOFpid) && //AliESDtrack::kTOFpid=0x8000
814 (tAodTrack->GetStatus() & AliESDtrack::kTOFout) && //AliESDtrack::kTOFout=0x2000
815 (tAodTrack->GetStatus() & AliESDtrack::kTIME) && //AliESDtrack::kTIME=0x80000000
816 !(tAodTrack->GetStatus() & AliESDtrack::kTOFmismatch)) //AliESDtrack::kTOFmismatch=0x100000
817 {
818 if(tAodTrack->IsOn(AliESDtrack::kTOFpid)) //AliESDtrack::kTOFpid=0x8000
819 {
820
821 nsigmaTOFPi = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kPion);
822 nsigmaTOFK = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kKaon);
823 nsigmaTOFP = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kProton);
824
825 Double_t len=200;// esdtrack->GetIntegratedLength(); !!!!!
826 Double_t tof=tAodTrack->GetTOFsignal();
827 if(tof > 0.) vp=len/tof/0.03;
828 }
829 }
830 tFemtoTrack->SetVTOF(vp);
831 tFemtoTrack->SetNSigmaTOFPi(nsigmaTOFPi);
832 tFemtoTrack->SetNSigmaTOFK(nsigmaTOFK);
833 tFemtoTrack->SetNSigmaTOFP(nsigmaTOFP);
834
835
836 //////////////////////////////////////
837
838}
839
948862a7 840void AliFemtoEventReaderAOD::SetCentralityPreSelection(double min, double max)
841{
842 fCentRange[0] = min; fCentRange[1] = max;
63073cb3 843 fUsePreCent = 1;
948862a7 844}
de568db1 845
846
1309aa0b 847void AliFemtoEventReaderAOD::SetAODpidUtil(AliAODpidUtil *aAODpidUtil)
848{
849 fAODpidUtil = aAODpidUtil;
4637c3d5 850 // printf("fAODpidUtil: %x\n",fAODpidUtil);
1309aa0b 851}
852
853
854
de568db1 855
856