]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/FEMTOSCOPY/AliFemto/AliFemtoEventReaderAOD.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
This patch covers
[u/mrichter/AliRoot.git] / PWGCF / FEMTOSCOPY / AliFemto / AliFemtoEventReaderAOD.cxx
CommitLineData
76ce4b5b 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"
17#include "AliAODMCHeader.h"
18#include "AliESDtrack.h"
19
20#include "AliFmPhysicalHelixD.h"
21#include "AliFmThreeVectorF.h"
22
23#include "SystemOfUnits.h"
24
25#include "AliFemtoEvent.h"
26#include "AliFemtoModelHiddenInfo.h"
27#include "AliFemtoModelGlobalHiddenInfo.h"
28#include "AliPID.h"
29
30#include "AliAODpidUtil.h"
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():
42 fNumberofEvent(0),
43 fCurEvent(0),
44 fEvent(0x0),
45 fAllTrue(160),
46 fAllFalse(160),
47 fFilterBit(0),
48 // fPWG2AODTracks(0x0),
49 fReadMC(0),
50 fUsePreCent(0),
51 fAODpidUtil(0),
52 fInputFile(" "),
53 fFileName(" "),
54 fTree(0x0),
55 fAodFile(0x0)
56{
57 // default constructor
58 fAllTrue.ResetAllBits(kTRUE);
59 fAllFalse.ResetAllBits(kFALSE);
60 fCentRange[0] = 0;
61 fCentRange[1] = 1000;
62}
63
64AliFemtoEventReaderAOD::AliFemtoEventReaderAOD(const AliFemtoEventReaderAOD &aReader) :
65 AliFemtoEventReader(),
66 fNumberofEvent(0),
67 fCurEvent(0),
68 fEvent(0x0),
69 fAllTrue(160),
70 fAllFalse(160),
71 fFilterBit(0),
72 // fPWG2AODTracks(0x0),
73 fReadMC(0),
74 fUsePreCent(0),
75 fAODpidUtil(0),
76 fInputFile(" "),
77 fFileName(" "),
78 fTree(0x0),
79 fAodFile(0x0)
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;
92 fAODpidUtil = aReader.fAODpidUtil;
93 fCentRange[0] = aReader.fCentRange[0];
94 fCentRange[1] = aReader.fCentRange[1];
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;
130 fAODpidUtil = aReader.fAODpidUtil;
131 fCentRange[0] = aReader.fCentRange[0];
132 fCentRange[1] = aReader.fCentRange[1];
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);
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 {
167 // cout<<"putting file "<<string(buffer)<<" into analysis"<<endl;
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;
184 cout<<"reader"<<endl;
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();
201 // cout<<"Number of Entries in file "<<fNumberofEvent<<endl;
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
214 // cout << "Read event " << fEvent << " from file " << fTree << endl;
215
216 hbtEvent = new AliFemtoEvent;
217
218 CopyAODtoFemtoEvent(hbtEvent);
219 fCurEvent++;
220
221
222 return hbtEvent;
223}
224
225void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
226{
227
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));
240 tEvent->SetZDCParticipants(0);
241 tEvent->SetTriggerMask(fEvent->GetTriggerMask());
242 tEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
243
244 // Attempt to access MC header
245 AliAODMCHeader *mcH;
246 TClonesArray *mcP=0;
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
259 tEvent->SetReactionPlaneAngle(fEvent->GetHeader()->GetQTheta(0)/2.0);
260
261 Int_t *motherids=0;
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
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();
292 cout<<"nofTracks: "<<nofTracks<<endl;
293
294 AliCentrality *cent = fEvent->GetCentrality();
295
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;
301
302 return;
303 }
304 }
305
306 int realnofTracks=0; // number of track which we use in a analysis
307 int tracksPrim=0;
308
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
320 // int tNormMult = 0;
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();
334// if (!aodtrack->TestFilterBit(fFilterBit)) {
335// delete trackCopy;
336// continue;
337// }
338
339
340// if (aodtrack->IsOn(AliESDtrack::kTPCrefit))
341// if (aodtrack->Chi2perNDF() < 6.0)
342// if (aodtrack->Eta() < 0.9)
343// tNormMult++;
344
345
346// CopyAODtoFemtoTrack(aodtrack, trackCopy, pwg2aodtrack);
347
348// if (mcP) {
349// // Fill the hidden information with the simulated data
350// // Int_t pLabel = aodtrack->GetLabel();
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
403// // if (fRotateToEventPlane) {
404// // double tPhi = TMath::ATan2(fpy, fpx);
405// // double tRad = TMath::Hypot(fpx, fpy);
406
407// // fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
408// // fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
409// // }
410
411// tInfo->SetPDGPid(tPart->GetPdgCode());
412
413// // if (fRotateToEventPlane) {
414// // double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
415// // double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
416
417// // tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
418// // tRad*TMath::Sin(tPhi - tReactionPlane),
419// // tPart->Pz());
420// // }
421// // else
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
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
441// if (ktP.Mag() == 0) {
442// delete trackCopy;
443// continue;
444// }
445// }
446// else {
447 // No additional information exists
448 // Read in the normal AliAODTracks
449
450 // const AliAODTrack *aodtrack=fEvent->GetTrack(i); // getting the AODtrack directly
451 AliAODTrack *aodtrack=fEvent->GetTrack(i); // getting the AODtrack directly
452
453 if (aodtrack->IsPrimaryCandidate()) tracksPrim++;
454
455 if (!aodtrack->TestFilterBit(fFilterBit)) {
456 delete trackCopy;
457 continue;
458 }
459
460 CopyAODtoFemtoTrack(aodtrack, trackCopy);
461
462 // copying PID information from the correspondent track
463 // const AliAODTrack *aodtrackpid = fEvent->GetTrack(labels[-1-fEvent->GetTrack(i)->GetID()]);
464 AliAODTrack *aodtrackpid = fEvent->GetTrack(labels[-1-fEvent->GetTrack(i)->GetID()]);
465 CopyPIDtoFemtoTrack(aodtrackpid, trackCopy);
466
467 if (mcP) {
468 // Fill the hidden information with the simulated data
469 // Int_t pLabel = aodtrack->GetLabel();
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
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
573 if (ktP.Mag() == 0) {
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
585 tEvent->SetNormalizedMult(tracksPrim);
586
587 // AliCentrality *cent = fEvent->GetCentrality();
588 if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("V0M")));
589 // if (cent) tEvent->SetNormalizedMult((int) cent->GetCentralityPercentile("V0M"));
590
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
599 if (mcP) delete [] motherids;
600
601 cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
602}
603
604void AliFemtoEventReaderAOD::CopyAODtoFemtoTrack(AliAODTrack *tAodTrack,
605 AliFemtoTrack *tFemtoTrack
606 // AliPWG2AODTrack *tPWG2AODTrack
607 )
608{
609 // Copy the track information from the AOD into the internal AliFemtoTrack
610 // If it exists, use the additional information from the PWG2 AOD
611
612 // Primary Vertex position
613 double fV1[3];
614 fEvent->GetPrimaryVertex()->GetPosition(fV1);
615 // fEvent->GetPrimaryVertex()->GetXYZ(fV1);
616
617 tFemtoTrack->SetCharge(tAodTrack->Charge());
618
619 double pxyz[3];
620 tAodTrack->PxPyPz(pxyz);//reading noconstrained momentum
621 AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
622 tFemtoTrack->SetP(v);//setting momentum
623 tFemtoTrack->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
624 const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
625 //setting track helix
626 const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
627 AliFmPhysicalHelixD helix(ktP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(tFemtoTrack->Charge()));
628 tFemtoTrack->SetHelix(helix);
629
630 // Flags
631 tFemtoTrack->SetTrackId(tAodTrack->GetID());
632 tFemtoTrack->SetFlags(tAodTrack->GetFlags());
633 tFemtoTrack->SetLabel(tAodTrack->GetLabel());
634
635 // Track quality information
636 float covmat[6];
637 tAodTrack->GetCovMatrix(covmat);
638
639 double impact[2];
640 double covimpact[3];
641
642 if (!tAodTrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
643 //cout << "sth went wrong with dca propagation" << endl;
644 tFemtoTrack->SetImpactD(-1000.0);
645 tFemtoTrack->SetImpactZ(-1000.0);
646
647 }
648 else {
649 tFemtoTrack->SetImpactD(impact[0]);
650 tFemtoTrack->SetImpactZ(impact[1]);
651 }
652
653 // if (TMath::Abs(tAodTrack->Xv()) > 0.00000000001)
654 // tFemtoTrack->SetImpactD(TMath::Hypot(tAodTrack->Xv(), tAodTrack->Yv())*(tAodTrack->Xv()/TMath::Abs(tAodTrack->Xv())));
655 // else
656 // tFemtoTrack->SetImpactD(0.0);
657 // tFemtoTrack->SetImpactD(tAodTrack->DCA());
658
659 // tFemtoTrack->SetImpactZ(tAodTrack->ZAtDCA());
660
661
662 // tFemtoTrack->SetImpactD(TMath::Hypot(tAodTrack->Xv() - fV1[0], tAodTrack->Yv() - fV1[1]));
663 // tFemtoTrack->SetImpactZ(tAodTrack->Zv() - fV1[2]);
664
665
666 // cout
667 // << "dca" << TMath::Hypot(tAodTrack->Xv() - fV1[0], tAodTrack->Yv() - fV1[1])
668 // << "xv - fv10 = "<< tAodTrack->Xv() - fV1[0]
669 // << tAodTrack->Yv() - fV1[1]
670// << "xv = " << tAodTrack->Xv() << endl
671// << "fv1[0] = " << fV1[0] << endl
672// << "yv = " << tAodTrack->Yv() << endl
673// << "fv1[1] = " << fV1[1] << endl
674// << "zv = " << tAodTrack->Zv() << endl
675// << "fv1[2] = " << fV1[2] << endl
676// << "impact[0] = " << impact[0] << endl
677// << "impact[1] = " << impact[1] << endl
678// << endl << endl ;
679
680 tFemtoTrack->SetCdd(covmat[0]);
681 tFemtoTrack->SetCdz(covmat[1]);
682 tFemtoTrack->SetCzz(covmat[2]);
683 tFemtoTrack->SetITSchi2(tAodTrack->Chi2perNDF());
684 tFemtoTrack->SetITSncls(tAodTrack->GetITSNcls());
685 tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF());
686 tFemtoTrack->SetTPCncls(tAodTrack->GetTPCNcls());
687 tFemtoTrack->SetTPCnclsF(tAodTrack->GetTPCNcls());
688 tFemtoTrack->SetTPCsignalN(1);
689 tFemtoTrack->SetTPCsignalS(1);
690 tFemtoTrack->SetTPCsignal(tAodTrack->GetTPCsignal());
691
692// if (tPWG2AODTrack) {
693// // Copy the PWG2 specific information if it exists
694// tFemtoTrack->SetTPCClusterMap(tPWG2AODTrack->GetTPCClusterMap());
695// tFemtoTrack->SetTPCSharedMap(tPWG2AODTrack->GetTPCSharedMap());
696
697// double xtpc[3] = {0,0,0};
698// tPWG2AODTrack->GetTPCNominalEntrancePoint(xtpc);
699// tFemtoTrack->SetNominalTPCEntrancePoint(xtpc);
700// tPWG2AODTrack->GetTPCNominalExitPoint(xtpc);
701// tFemtoTrack->SetNominalTPCExitPoint(xtpc);
702// }
703// else {
704 // If not use dummy values
705 tFemtoTrack->SetTPCClusterMap(tAodTrack->GetTPCClusterMap());
706 tFemtoTrack->SetTPCSharedMap(tAodTrack->GetTPCSharedMap());
707
708 double xtpc[3] = {0,0,0};
709 tFemtoTrack->SetNominalTPCEntrancePoint(xtpc);
710 tFemtoTrack->SetNominalTPCExitPoint(xtpc);
711 // }
712
713 // // cout << "Track has " << TMath::Hypot(tAodTrack->Xv(), tAodTrack->Yv()) << " " << tAodTrack->Zv() << " " << tAodTrack->GetTPCNcls() << endl;
714
715
716 int indexes[3];
717 for (int ik=0; ik<3; ik++) {
718 indexes[ik] = 0;
719 }
720 tFemtoTrack->SetKinkIndexes(indexes);
721}
722
723void AliFemtoEventReaderAOD::SetFilterBit(UInt_t ibit)
724{
725 fFilterBit = (1 << (ibit));
726}
727
728void AliFemtoEventReaderAOD::SetReadMC(unsigned char a)
729{
730 fReadMC = a;
731}
732
733AliAODMCParticle* AliFemtoEventReaderAOD::GetParticleWithLabel(TClonesArray *mcP, Int_t aLabel)
734{
735 if (aLabel < 0) return 0;
736 AliAODMCParticle *aodP;
737 Int_t posstack = 0;
738 if (aLabel > mcP->GetEntries())
739 posstack = mcP->GetEntries();
740 else
741 posstack = aLabel;
742
743 aodP = (AliAODMCParticle *) mcP->At(posstack);
744 if (aodP->GetLabel() > posstack) {
745 do {
746 aodP = (AliAODMCParticle *) mcP->At(posstack);
747 if (aodP->GetLabel() == aLabel) return aodP;
748 posstack--;
749 }
750 while (posstack > 0);
751 }
752 else {
753 do {
754 aodP = (AliAODMCParticle *) mcP->At(posstack);
755 if (aodP->GetLabel() == aLabel) return aodP;
756 posstack++;
757 }
758 while (posstack < mcP->GetEntries());
759 }
760
761 return 0;
762}
763
764void AliFemtoEventReaderAOD::CopyPIDtoFemtoTrack(AliAODTrack *tAodTrack,
765 AliFemtoTrack *tFemtoTrack)
766{
767 double aodpid[10];
768 tAodTrack->GetPID(aodpid);
769 tFemtoTrack->SetPidProbElectron(aodpid[0]);
770 tFemtoTrack->SetPidProbMuon(aodpid[1]);
771 tFemtoTrack->SetPidProbPion(aodpid[2]);
772 tFemtoTrack->SetPidProbKaon(aodpid[3]);
773 tFemtoTrack->SetPidProbProton(aodpid[4]);
774
775 aodpid[0] = -100000.0;
776 aodpid[1] = -100000.0;
777 aodpid[2] = -100000.0;
778 aodpid[3] = -100000.0;
779 aodpid[4] = -100000.0;
780
781 double tTOF = 0.0;
782
783 if (tAodTrack->GetStatus() & AliESDtrack::kTOFpid) { //AliESDtrack::kTOFpid=0x8000
784 tTOF = tAodTrack->GetTOFsignal();
785 tAodTrack->GetIntegratedTimes(aodpid);
786 }
787
788 tFemtoTrack->SetTofExpectedTimes(tTOF-aodpid[2], tTOF-aodpid[3], tTOF-aodpid[4]);
789
790 ////// TPC ////////////////////////////////////////////
791
792 float nsigmaTPCK=-1000.;
793 float nsigmaTPCPi=-1000.;
794 float nsigmaTPCP=-1000.;
795
796 // cout<<"in reader fESDpid"<<fESDpid<<endl;
797
798 if (tAodTrack->IsOn(AliESDtrack::kTPCpid)){ //AliESDtrack::kTPCpid=0x0080
799 nsigmaTPCK = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kKaon);
800 nsigmaTPCPi = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kPion);
801 nsigmaTPCP = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kProton);
802 }
803
804 tFemtoTrack->SetNSigmaTPCPi(nsigmaTPCPi);
805 tFemtoTrack->SetNSigmaTPCK(nsigmaTPCK);
806 tFemtoTrack->SetNSigmaTPCP(nsigmaTPCP);
807
808 tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF());
809 tFemtoTrack->SetTPCncls(tAodTrack->GetTPCNcls());
810 tFemtoTrack->SetTPCnclsF(tAodTrack->GetTPCNcls());
811
812 tFemtoTrack->SetTPCsignalN(1);
813 tFemtoTrack->SetTPCsignalS(1);
814 tFemtoTrack->SetTPCsignal(tAodTrack->GetTPCsignal());
815
816 ///////TOF//////////////////////
817
818 float vp=-1000.;
819 float nsigmaTOFPi=-1000.;
820 float nsigmaTOFK=-1000.;
821 float nsigmaTOFP=-1000.;
822
823 if ((tAodTrack->GetStatus() & AliESDtrack::kTOFpid) && //AliESDtrack::kTOFpid=0x8000
824 (tAodTrack->GetStatus() & AliESDtrack::kTOFout) && //AliESDtrack::kTOFout=0x2000
825 (tAodTrack->GetStatus() & AliESDtrack::kTIME) && //AliESDtrack::kTIME=0x80000000
826 !(tAodTrack->GetStatus() & AliESDtrack::kTOFmismatch)) //AliESDtrack::kTOFmismatch=0x100000
827 {
828 if(tAodTrack->IsOn(AliESDtrack::kTOFpid)) //AliESDtrack::kTOFpid=0x8000
829 {
830
831 nsigmaTOFPi = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kPion);
832 nsigmaTOFK = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kKaon);
833 nsigmaTOFP = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kProton);
834
835 Double_t len=200;// esdtrack->GetIntegratedLength(); !!!!!
836 Double_t tof=tAodTrack->GetTOFsignal();
837 if(tof > 0.) vp=len/tof/0.03;
838 }
839 }
840 tFemtoTrack->SetVTOF(vp);
841 tFemtoTrack->SetNSigmaTOFPi(nsigmaTOFPi);
842 tFemtoTrack->SetNSigmaTOFK(nsigmaTOFK);
843 tFemtoTrack->SetNSigmaTOFP(nsigmaTOFP);
844
845
846 //////////////////////////////////////
847
848}
849
850void AliFemtoEventReaderAOD::SetCentralityPreSelection(double min, double max)
851{
852 fCentRange[0] = min; fCentRange[1] = max;
853 fUsePreCent = 1;
854}
855
856
857void AliFemtoEventReaderAOD::SetAODpidUtil(AliAODpidUtil *aAODpidUtil)
858{
859 fAODpidUtil = aAODpidUtil;
860 // printf("fAODpidUtil: %x\n",fAODpidUtil);
861}
862
863
864
865
866