]> 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
updated
[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;
973a91f8 39
40double fV1[3];
41
76ce4b5b 42//____________________________
43//constructor with 0 parameters , look at default settings
44AliFemtoEventReaderAOD::AliFemtoEventReaderAOD():
45 fNumberofEvent(0),
46 fCurEvent(0),
47 fEvent(0x0),
48 fAllTrue(160),
49 fAllFalse(160),
50 fFilterBit(0),
64536eaf 51 fFilterMask(0),
76ce4b5b 52 // fPWG2AODTracks(0x0),
53 fReadMC(0),
973a91f8 54 fReadV0(0),
76ce4b5b 55 fUsePreCent(0),
1777446b 56 fEstEventMult(kCentrality),
76ce4b5b 57 fAODpidUtil(0),
1777446b 58 fAODheader(0),
76ce4b5b 59 fInputFile(" "),
60 fFileName(" "),
61 fTree(0x0),
ba3c23a4 62 fAodFile(0x0),
5e2038a9 63 fMagFieldSign(1),
64 fisEPVZ(kTRUE)
76ce4b5b 65{
66 // default constructor
67 fAllTrue.ResetAllBits(kTRUE);
68 fAllFalse.ResetAllBits(kFALSE);
69 fCentRange[0] = 0;
70 fCentRange[1] = 1000;
71}
72
73AliFemtoEventReaderAOD::AliFemtoEventReaderAOD(const AliFemtoEventReaderAOD &aReader) :
74 AliFemtoEventReader(),
75 fNumberofEvent(0),
76 fCurEvent(0),
77 fEvent(0x0),
78 fAllTrue(160),
79 fAllFalse(160),
80 fFilterBit(0),
64536eaf 81 fFilterMask(0),
76ce4b5b 82 // fPWG2AODTracks(0x0),
83 fReadMC(0),
973a91f8 84 fReadV0(0),
76ce4b5b 85 fUsePreCent(0),
1777446b 86 fEstEventMult(kCentrality),
76ce4b5b 87 fAODpidUtil(0),
1777446b 88 fAODheader(0),
76ce4b5b 89 fInputFile(" "),
90 fFileName(" "),
91 fTree(0x0),
ba3c23a4 92 fAodFile(0x0),
5e2038a9 93 fMagFieldSign(1),
94 fisEPVZ(kTRUE)
76ce4b5b 95{
96 // copy constructor
973a91f8 97 fReadMC = aReader.fReadMC;
98 fReadV0 = aReader.fReadV0;
76ce4b5b 99 fInputFile = aReader.fInputFile;
100 fFileName = aReader.fFileName;
101 fNumberofEvent = aReader.fNumberofEvent;
102 fCurEvent = aReader.fCurEvent;
103 fEvent = new AliAODEvent();
104 fAodFile = new TFile(aReader.fAodFile->GetName());
105 fAllTrue.ResetAllBits(kTRUE);
106 fAllFalse.ResetAllBits(kFALSE);
107 fFilterBit = aReader.fFilterBit;
108 // fPWG2AODTracks = aReader.fPWG2AODTracks;
109 fAODpidUtil = aReader.fAODpidUtil;
1777446b 110 fAODheader = aReader.fAODheader;
76ce4b5b 111 fCentRange[0] = aReader.fCentRange[0];
112 fCentRange[1] = aReader.fCentRange[1];
1777446b 113 fEstEventMult = aReader.fEstEventMult;
973a91f8 114 fUsePreCent = aReader.fUsePreCent;
76ce4b5b 115}
116//__________________
117//Destructor
118AliFemtoEventReaderAOD::~AliFemtoEventReaderAOD()
119{
120 // destructor
121 delete fTree;
122 delete fEvent;
123 delete fAodFile;
124// if (fPWG2AODTracks) {
125// fPWG2AODTracks->Delete();
126// delete fPWG2AODTracks;
127// }
128}
129
130//__________________
131AliFemtoEventReaderAOD& AliFemtoEventReaderAOD::operator=(const AliFemtoEventReaderAOD& aReader)
132{
133 // assignment operator
134 if (this == &aReader)
1777446b 135 return *this;
76ce4b5b 136
137 fInputFile = aReader.fInputFile;
138 fFileName = aReader.fFileName;
139 fNumberofEvent = aReader.fNumberofEvent;
140 fCurEvent = aReader.fCurEvent;
141 if (fTree) delete fTree;
142 if (fEvent) delete fEvent;
143 fEvent = new AliAODEvent();
144 if (fAodFile) delete fAodFile;
145 fAodFile = new TFile(aReader.fAodFile->GetName());
146 fAllTrue.ResetAllBits(kTRUE);
147 fAllFalse.ResetAllBits(kFALSE);
148 fFilterBit = aReader.fFilterBit;
64536eaf 149 fFilterMask = aReader.fFilterMask;
76ce4b5b 150 // fPWG2AODTracks = aReader.fPWG2AODTracks;
151 fAODpidUtil = aReader.fAODpidUtil;
1777446b 152 fAODheader = aReader.fAODheader;
76ce4b5b 153 fCentRange[0] = aReader.fCentRange[0];
154 fCentRange[1] = aReader.fCentRange[1];
973a91f8 155 fUsePreCent = aReader.fUsePreCent;
1777446b 156 fEstEventMult = aReader.fEstEventMult;
76ce4b5b 157
158 return *this;
159}
160//__________________
161AliFemtoString AliFemtoEventReaderAOD::Report()
162{
163 // create reader report
164 AliFemtoString temp = "\n This is the AliFemtoEventReaderAOD\n";
165 return temp;
166}
167
168//__________________
169void AliFemtoEventReaderAOD::SetInputFile(const char* inputFile)
170{
171 //setting the name of file where names of AOD file are written
172 //it takes only this files which have good trees
173 char buffer[256];
174 fInputFile=string(inputFile);
175 ifstream infile(inputFile);
176
177 fTree = new TChain("aodTree");
178
179 if(infile.good()==true)
180 {
181 //checking if all give files have good tree inside
182 while (infile.eof()==false)
183 {
184 infile.getline(buffer,256);
185 TFile *aodFile=TFile::Open(buffer,"READ");
186 if (aodFile!=0x0)
187 {
188 TTree* tree = (TTree*) aodFile->Get("aodTree");
189 if (tree!=0x0)
190 {
191 // cout<<"putting file "<<string(buffer)<<" into analysis"<<endl;
192 fTree->AddFile(buffer);
193 delete tree;
194 }
195 aodFile->Close();
196 }
197 delete aodFile;
198 }
199 }
200}
201
202AliFemtoEvent* AliFemtoEventReaderAOD::ReturnHbtEvent()
203{
204 // read in a next hbt event from the chain
205 // convert it to AliFemtoEvent and return
206 // for further analysis
207 AliFemtoEvent *hbtEvent = 0;
5e2038a9 208 // cout<<"reader"<<endl;
76ce4b5b 209 if (fCurEvent==fNumberofEvent)//open next file
210 {
211 if(fNumberofEvent==0)
212 {
213 fEvent=new AliAODEvent();
214 fEvent->ReadFromTree(fTree);
215
216 // Check for the existence of the additional information
217// fPWG2AODTracks = (TClonesArray *) fEvent->GetList()->FindObject("pwg2aodtracks");
218
219// if (fPWG2AODTracks) {
220// cout << "Found additional PWG2 specific information in the AOD!" << endl;
221// cout << "Reading only tracks with the additional information" << endl;
222// }
223
224 fNumberofEvent=fTree->GetEntries();
225 // cout<<"Number of Entries in file "<<fNumberofEvent<<endl;
226 fCurEvent=0;
227 }
228 else //no more data to read
229 {
5e2038a9 230 // cout<<"no more files "<<hbtEvent<<endl;
76ce4b5b 231 fReaderStatus=1;
232 return hbtEvent;
233 }
234 }
235
5e2038a9 236 // cout<<"starting to read event "<<fCurEvent<<endl;
76ce4b5b 237 fTree->GetEvent(fCurEvent);//getting next event
238 // cout << "Read event " << fEvent << " from file " << fTree << endl;
239
240 hbtEvent = new AliFemtoEvent;
241
242 CopyAODtoFemtoEvent(hbtEvent);
243 fCurEvent++;
244
245
246 return hbtEvent;
247}
248
249void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
250{
251
252 // A function that reads in the AOD event
253 // and transfers the neccessary information into
254 // the internal AliFemtoEvent
255
256 // setting global event characteristics
257 tEvent->SetRunNumber(fEvent->GetRunNumber());
258 tEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
259 tEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
260 tEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
261 tEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
262 tEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
263 tEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy(0));
264 tEvent->SetZDCParticipants(0);
265 tEvent->SetTriggerMask(fEvent->GetTriggerMask());
266 tEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
267
268 // Attempt to access MC header
269 AliAODMCHeader *mcH;
270 TClonesArray *mcP=0;
271 if (fReadMC) {
272 mcH = (AliAODMCHeader *) fEvent->FindListObject(AliAODMCHeader::StdBranchName());
273 if (!mcH) {
274 cout << "AOD MC information requested, but no header found!" << endl;
275 }
276
277 mcP = (TClonesArray *) fEvent->FindListObject(AliAODMCParticle::StdBranchName());
278 if (!mcP) {
279 cout << "AOD MC information requested, but no particle array found!" << endl;
280 }
281 }
282
283 tEvent->SetReactionPlaneAngle(fEvent->GetHeader()->GetQTheta(0)/2.0);
284
285 Int_t *motherids=0;
286 if (mcP) {
287 motherids = new Int_t[((AliAODMCParticle *) mcP->At(mcP->GetEntries()-1))->GetLabel()];
288 for (int ip=0; ip<mcP->GetEntries(); ip++) motherids[ip] = 0;
289
290 // Read in mother ids
291 AliAODMCParticle *motherpart;
292 for (int ip=0; ip<mcP->GetEntries(); ip++) {
293 motherpart = (AliAODMCParticle *) mcP->At(ip);
294 if (motherpart->GetDaughter(0) > 0)
295 motherids[motherpart->GetDaughter(0)] = ip;
296 if (motherpart->GetDaughter(1) > 0)
297 motherids[motherpart->GetDaughter(1)] = ip;
298 }
299 }
300
301 // Primary Vertex position
973a91f8 302 // double fV1[3];
76ce4b5b 303 fEvent->GetPrimaryVertex()->GetPosition(fV1);
304
305 AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
306 tEvent->SetPrimVertPos(vertex);
307
308 //starting to reading tracks
309 int nofTracks=0; //number of reconstructed tracks in event
310
311 // Check to see whether the additional info exists
1777446b 312 // if (fPWG2AODTracks)
313 // nofTracks=fPWG2AODTracks->GetEntries();
314 // else
315 nofTracks=fEvent->GetNumberOfTracks();
316
317 AliEventplane *ep = fEvent->GetEventplane();
318 if (ep) {
319 tEvent->SetEP(ep);
5e2038a9 320 if (fisEPVZ)
321 tEvent->SetReactionPlaneAngle(ep->GetEventplane("V0",fEvent,2));
322 else
1777446b 323 tEvent->SetReactionPlaneAngle(ep->GetEventplane("Q"));
324 }
76ce4b5b 325
36cfef5c 326 AliCentrality *cent = fEvent->GetCentrality();
4eac0b05 327
1777446b 328 if (!fEstEventMult && cent && fUsePreCent) {
76ce4b5b 329 if ((cent->GetCentralityPercentile("V0M")*10 < fCentRange[0]) ||
330 (cent->GetCentralityPercentile("V0M")*10 > fCentRange[1]))
331 {
5e2038a9 332 // cout << "Centrality " << cent->GetCentralityPercentile("V0M") << " outside of preselection range " << fCentRange[0] << " - " << fCentRange[1] << endl;
76ce4b5b 333
334 return;
335 }
336 }
337
338 int realnofTracks=0; // number of track which we use in a analysis
339 int tracksPrim=0;
340
341 int labels[20000];
342 for (int il=0; il<20000; il++) labels[il] = -1;
343
344 // looking for global tracks and saving their numbers to copy from them PID information to TPC-only tracks in the main loop over tracks
345 for (int i=0;i<nofTracks;i++) {
346 const AliAODTrack *aodtrack=fEvent->GetTrack(i);
347 if (!aodtrack->TestFilterBit(fFilterBit)) {
973a91f8 348 if(aodtrack->GetID() < 0) continue;
76ce4b5b 349 labels[aodtrack->GetID()] = i;
350 }
351 }
352
ba3c23a4 353 int tNormMult = 0;
76ce4b5b 354 for (int i=0;i<nofTracks;i++)
355 {
356 AliFemtoTrack* trackCopy = new AliFemtoTrack();
357
358// if (fPWG2AODTracks) {
359// // Read tracks from the additional pwg2 specific AOD part
360// // if they exist
361// // Note that in that case all the AOD tracks without the
362// // additional information will be ignored !
363// AliPWG2AODTrack *pwg2aodtrack = (AliPWG2AODTrack *) fPWG2AODTracks->At(i);
364
365// // Getting the AOD track through the ref of the additional info
366// AliAODTrack *aodtrack = pwg2aodtrack->GetRefAODTrack();
367// if (!aodtrack->TestFilterBit(fFilterBit)) {
368// delete trackCopy;
369// continue;
370// }
371
372
373// if (aodtrack->IsOn(AliESDtrack::kTPCrefit))
374// if (aodtrack->Chi2perNDF() < 6.0)
375// if (aodtrack->Eta() < 0.9)
376// tNormMult++;
377
378
379// CopyAODtoFemtoTrack(aodtrack, trackCopy, pwg2aodtrack);
380
381// if (mcP) {
382// // Fill the hidden information with the simulated data
383// // Int_t pLabel = aodtrack->GetLabel();
384// AliAODMCParticle *tPart = GetParticleWithLabel(mcP, (TMath::Abs(aodtrack->GetLabel())));
385
386// // Check the mother information
387
388// // Using the new way of storing the freeze-out information
389// // Final state particle is stored twice on the stack
390// // one copy (mother) is stored with original freeze-out information
391// // and is not tracked
392// // the other one (daughter) is stored with primary vertex position
393// // and is tracked
394
395// // Freeze-out coordinates
396// double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
397// fpx = tPart->Xv() - fV1[0];
398// fpy = tPart->Yv() - fV1[1];
399// fpz = tPart->Zv() - fV1[2];
400// fpt = tPart->T();
401
402// AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
403// tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
404
405// fpx *= 1e13;
406// fpy *= 1e13;
407// fpz *= 1e13;
408// fpt *= 1e13;
409
410// // cout << "Looking for mother ids " << endl;
411// if (motherids[TMath::Abs(aodtrack->GetLabel())]>0) {
412// // cout << "Got mother id" << endl;
413// AliAODMCParticle *mother = GetParticleWithLabel(mcP, motherids[TMath::Abs(aodtrack->GetLabel())]);
414// // Check if this is the same particle stored twice on the stack
415// if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
416// // It is the same particle
417// // Read in the original freeze-out information
418// // and convert it from to [fm]
419
420// // EPOS style
421// // fpx = mother->Xv()*1e13*0.197327;
422// // fpy = mother->Yv()*1e13*0.197327;
423// // fpz = mother->Zv()*1e13*0.197327;
424// // fpt = mother->T() *1e13*0.197327*0.5;
425
426
427// // Therminator style
428// fpx = mother->Xv()*1e13;
429// fpy = mother->Yv()*1e13;
430// fpz = mother->Zv()*1e13;
431// fpt = mother->T() *1e13*3e10;
432
433// }
434// }
435
436// // if (fRotateToEventPlane) {
437// // double tPhi = TMath::ATan2(fpy, fpx);
438// // double tRad = TMath::Hypot(fpx, fpy);
439
440// // fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
441// // fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
442// // }
443
444// tInfo->SetPDGPid(tPart->GetPdgCode());
445
446// // if (fRotateToEventPlane) {
447// // double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
448// // double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
449
450// // tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
451// // tRad*TMath::Sin(tPhi - tReactionPlane),
452// // tPart->Pz());
453// // }
454// // else
455// tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
456// Double_t mass2 = (tPart->E() *tPart->E() -
457// tPart->Px()*tPart->Px() -
458// tPart->Py()*tPart->Py() -
459// tPart->Pz()*tPart->Pz());
460// if (mass2>0.0)
461// tInfo->SetMass(TMath::Sqrt(mass2));
462// else
463// tInfo->SetMass(0.0);
464
465// tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
466// trackCopy->SetHiddenInfo(tInfo);
467
468// }
469
470// double pxyz[3];
471// aodtrack->PxPyPz(pxyz);//reading noconstarined momentum
472// const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
473// // Check the sanity of the tracks - reject zero momentum tracks
474// if (ktP.Mag() == 0) {
475// delete trackCopy;
476// continue;
477// }
478// }
479// else {
480 // No additional information exists
481 // Read in the normal AliAODTracks
482
483 // const AliAODTrack *aodtrack=fEvent->GetTrack(i); // getting the AODtrack directly
484 AliAODTrack *aodtrack=fEvent->GetTrack(i); // getting the AODtrack directly
ba3c23a4 485
4eac0b05 486
487
76ce4b5b 488 if (aodtrack->IsPrimaryCandidate()) tracksPrim++;
489
64536eaf 490 if (fFilterBit && !aodtrack->TestFilterBit(fFilterBit)) {
76ce4b5b 491 delete trackCopy;
492 continue;
64536eaf 493 }
494
495 if (fFilterMask && !aodtrack->TestFilterBit(fFilterMask)) {
496 delete trackCopy;
497 continue;
498 }
499
ba3c23a4 500 //counting particles to set multiplicity
501 double impact[2];
502 double covimpact[3];
503 if (aodtrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
1777446b 504 if(impact[0]<0.2 && TMath::Abs(impact[1]+fV1[2])<2.0)
505 //if (aodtrack->IsPrimaryCandidate()) //? instead of kinks?
ba3c23a4 506 if (aodtrack->Chi2perNDF() < 4.0)
507 if (aodtrack->Pt() > 0.15 && aodtrack->Pt() < 20)
1777446b 508 if (aodtrack->GetTPCNcls() > 70)
ba3c23a4 509 if (aodtrack->Eta() < 0.8)
510 tNormMult++;
ba3c23a4 511 }
512
76ce4b5b 513 CopyAODtoFemtoTrack(aodtrack, trackCopy);
514
515 // copying PID information from the correspondent track
516 // const AliAODTrack *aodtrackpid = fEvent->GetTrack(labels[-1-fEvent->GetTrack(i)->GetID()]);
ba3c23a4 517
518
519 AliAODTrack *aodtrackpid;
64536eaf 520 if((fFilterBit == (1 << (7))) || fFilterMask==128) //for TPC Only tracks we have to copy PID information from corresponding global tracks
ba3c23a4 521 aodtrackpid = fEvent->GetTrack(labels[-1-fEvent->GetTrack(i)->GetID()]);
522 else
523 aodtrackpid = fEvent->GetTrack(i);
63e066f7 524 CopyPIDtoFemtoTrack(aodtrackpid, trackCopy);
973a91f8 525
76ce4b5b 526 if (mcP) {
527 // Fill the hidden information with the simulated data
528 // Int_t pLabel = aodtrack->GetLabel();
529 AliAODMCParticle *tPart = GetParticleWithLabel(mcP, (TMath::Abs(aodtrack->GetLabel())));
530
531 AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
532 double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
533 if (!tPart) {
534 fpx = fV1[0];
535 fpy = fV1[1];
536 fpz = fV1[2];
537 tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
538 tInfo->SetPDGPid(0);
539 tInfo->SetTrueMomentum(0.0, 0.0, 0.0);
540 tInfo->SetEmissionPoint(0.0, 0.0, 0.0, 0.0);
541 tInfo->SetMass(0);
542 }
543 else {
544 // Check the mother information
545
546 // Using the new way of storing the freeze-out information
547 // Final state particle is stored twice on the stack
548 // one copy (mother) is stored with original freeze-out information
549 // and is not tracked
550 // the other one (daughter) is stored with primary vertex position
551 // and is tracked
552
553 // Freeze-out coordinates
554 fpx = tPart->Xv() - fV1[0];
555 fpy = tPart->Yv() - fV1[1];
556 fpz = tPart->Zv() - fV1[2];
557 // fpt = tPart->T();
558
559 tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
560
561 fpx *= 1e13;
562 fpy *= 1e13;
563 fpz *= 1e13;
564 // fpt *= 1e13;
565
566 // cout << "Looking for mother ids " << endl;
567 if (motherids[TMath::Abs(aodtrack->GetLabel())]>0) {
568 // cout << "Got mother id" << endl;
569 AliAODMCParticle *mother = GetParticleWithLabel(mcP, motherids[TMath::Abs(aodtrack->GetLabel())]);
570 // Check if this is the same particle stored twice on the stack
571 if (mother) {
572 if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
573 // It is the same particle
574 // Read in the original freeze-out information
575 // and convert it from to [fm]
576
577 // EPOS style
578 // fpx = mother->Xv()*1e13*0.197327;
579 // fpy = mother->Yv()*1e13*0.197327;
580 // fpz = mother->Zv()*1e13*0.197327;
581 // fpt = mother->T() *1e13*0.197327*0.5;
582
583
584 // Therminator style
585 fpx = mother->Xv()*1e13;
586 fpy = mother->Yv()*1e13;
587 fpz = mother->Zv()*1e13;
588 // fpt = mother->T() *1e13*3e10;
589
590 }
591 }
592 }
593
594 // if (fRotateToEventPlane) {
595 // double tPhi = TMath::ATan2(fpy, fpx);
596 // double tRad = TMath::Hypot(fpx, fpy);
597
598 // fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
599 // fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
600 // }
601
602 tInfo->SetPDGPid(tPart->GetPdgCode());
603
604 // if (fRotateToEventPlane) {
605 // double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
606 // double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
607
608 // tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
609 // tRad*TMath::Sin(tPhi - tReactionPlane),
610 // tPart->Pz());
611 // }
612 // else
613 tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
614 Double_t mass2 = (tPart->E() *tPart->E() -
615 tPart->Px()*tPart->Px() -
616 tPart->Py()*tPart->Py() -
617 tPart->Pz()*tPart->Pz());
618 if (mass2>0.0)
619 tInfo->SetMass(TMath::Sqrt(mass2));
620 else
621 tInfo->SetMass(0.0);
622
623 tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
624 }
625 trackCopy->SetHiddenInfo(tInfo);
626 }
627
628 double pxyz[3];
7eb9f5d0 629
630 //AliExternalTrackParam *param = new AliExternalTrackParam(*aodtrack->GetInnerParam());
631 trackCopy->SetInnerMomentum(aodtrack->GetTPCmomentum());
632
76ce4b5b 633 aodtrack->PxPyPz(pxyz);//reading noconstarined momentum
634 const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
635 // Check the sanity of the tracks - reject zero momentum tracks
636 if (ktP.Mag() == 0) {
637 delete trackCopy;
638 continue;
639 }
640 // }
641
4eac0b05 642
76ce4b5b 643 tEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
644 realnofTracks++;//real number of tracks
645 }
646
647 tEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
648 tEvent->SetNormalizedMult(tracksPrim);
649
1777446b 650 if (fEstEventMult==kCentrality) {
973a91f8 651 // AliCentrality *cent = fEvent->GetCentrality();
4eac0b05 652 //cout<<"AliFemtoEventReaderAOD:"<<lrint(10*cent->GetCentralityPercentile("V0M"))<<endl;
973a91f8 653 if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("V0M")));
654 // if (cent) tEvent->SetNormalizedMult((int) cent->GetCentralityPercentile("V0M"));
655
656 if (cent) {
657 tEvent->SetCentralityV0(cent->GetCentralityPercentile("V0M"));
658 // tEvent->SetCentralityFMD(cent->GetCentralityPercentile("FMD"));
659 tEvent->SetCentralitySPD1(cent->GetCentralityPercentile("CL1"));
660 // tEvent->SetCentralityTrk(cent->GetCentralityPercentile("TRK"));
661 }
76ce4b5b 662 }
1777446b 663 else if(fEstEventMult==kGlobalCount){
ba3c23a4 664 tEvent->SetNormalizedMult(tNormMult); //particles counted in the loop, trying to reproduce GetReferenceMultiplicity. If better (default) method appears it should be changed
665 }
1777446b 666 else if(fEstEventMult==kReference)
667 {
668 tEvent->SetNormalizedMult(fAODheader->GetRefMultiplicity());
669 }
670 else if(fEstEventMult==kTPCOnlyRef)
671 {
672 tEvent->SetNormalizedMult(fAODheader->GetTPConlyRefMultiplicity());
673 }
1938eadf 674 else if(fEstEventMult == kVZERO)
675 {
676 Float_t multV0 = 0;
677 for (Int_t i=0; i<64; i++)
678 multV0 += fEvent->GetVZEROData()->GetMultiplicity(i);
679 tEvent->SetNormalizedMult(multV0);
680 }
76ce4b5b 681
682 if (mcP) delete [] motherids;
683
5e2038a9 684 // cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
973a91f8 685
686 if(fReadV0)
687 {
ce7b3d98 688 int count_pass = 0;
973a91f8 689 for (Int_t i = 0; i < fEvent->GetNumberOfV0s(); i++) {
690 AliAODv0* aodv0 = fEvent->GetV0(i);
691 if (!aodv0) continue;
692 if(aodv0->GetNDaughters()>2) continue;
693 if(aodv0->GetNProngs()>2) continue;
694 if(aodv0->GetCharge()!=0) continue;
695 if(aodv0->ChargeProng(0)==aodv0->ChargeProng(1)) continue;
ce7b3d98 696 if(aodv0->CosPointingAngle(fV1)<0.998) continue;
973a91f8 697 AliFemtoV0* trackCopyV0 = new AliFemtoV0();
ce7b3d98 698 count_pass++;
973a91f8 699 CopyAODtoFemtoV0(aodv0, trackCopyV0);
700 tEvent->V0Collection()->push_back(trackCopyV0);
701 //cout<<"Pushback v0 to v0collection"<<endl;
702 }
703 }
704
76ce4b5b 705}
706
707void AliFemtoEventReaderAOD::CopyAODtoFemtoTrack(AliAODTrack *tAodTrack,
708 AliFemtoTrack *tFemtoTrack
709 // AliPWG2AODTrack *tPWG2AODTrack
710 )
711{
712 // Copy the track information from the AOD into the internal AliFemtoTrack
713 // If it exists, use the additional information from the PWG2 AOD
714
715 // Primary Vertex position
973a91f8 716
76ce4b5b 717 fEvent->GetPrimaryVertex()->GetPosition(fV1);
718 // fEvent->GetPrimaryVertex()->GetXYZ(fV1);
719
720 tFemtoTrack->SetCharge(tAodTrack->Charge());
721
722 double pxyz[3];
723 tAodTrack->PxPyPz(pxyz);//reading noconstrained momentum
724 AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
725 tFemtoTrack->SetP(v);//setting momentum
726 tFemtoTrack->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
727 const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
728 //setting track helix
729 const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
730 AliFmPhysicalHelixD helix(ktP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(tFemtoTrack->Charge()));
731 tFemtoTrack->SetHelix(helix);
732
733 // Flags
734 tFemtoTrack->SetTrackId(tAodTrack->GetID());
735 tFemtoTrack->SetFlags(tAodTrack->GetFlags());
736 tFemtoTrack->SetLabel(tAodTrack->GetLabel());
737
738 // Track quality information
739 float covmat[6];
740 tAodTrack->GetCovMatrix(covmat);
741
41cb6c1e 742 // ! DCA information is done in CopyPIDtoFemtoTrack()
76ce4b5b 743
41cb6c1e 744 // double impact[2];
745 // double covimpact[3];
76ce4b5b 746
41cb6c1e 747 // if (!tAodTrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
748 // //cout << "sth went wrong with dca propagation" << endl;
749 // tFemtoTrack->SetImpactD(-1000.0);
750 // tFemtoTrack->SetImpactZ(-1000.0);
751
752 // }
753 // else {
754 // tFemtoTrack->SetImpactD(impact[0]);
755 // tFemtoTrack->SetImpactZ(impact[1]+fV1[2]);
756 // }
76ce4b5b 757
758 // if (TMath::Abs(tAodTrack->Xv()) > 0.00000000001)
759 // tFemtoTrack->SetImpactD(TMath::Hypot(tAodTrack->Xv(), tAodTrack->Yv())*(tAodTrack->Xv()/TMath::Abs(tAodTrack->Xv())));
760 // else
761 // tFemtoTrack->SetImpactD(0.0);
762 // tFemtoTrack->SetImpactD(tAodTrack->DCA());
763
764 // tFemtoTrack->SetImpactZ(tAodTrack->ZAtDCA());
765
766
767 // tFemtoTrack->SetImpactD(TMath::Hypot(tAodTrack->Xv() - fV1[0], tAodTrack->Yv() - fV1[1]));
768 // tFemtoTrack->SetImpactZ(tAodTrack->Zv() - fV1[2]);
769
770
771 // cout
772 // << "dca" << TMath::Hypot(tAodTrack->Xv() - fV1[0], tAodTrack->Yv() - fV1[1])
773 // << "xv - fv10 = "<< tAodTrack->Xv() - fV1[0]
774 // << tAodTrack->Yv() - fV1[1]
775// << "xv = " << tAodTrack->Xv() << endl
776// << "fv1[0] = " << fV1[0] << endl
777// << "yv = " << tAodTrack->Yv() << endl
778// << "fv1[1] = " << fV1[1] << endl
779// << "zv = " << tAodTrack->Zv() << endl
780// << "fv1[2] = " << fV1[2] << endl
781// << "impact[0] = " << impact[0] << endl
782// << "impact[1] = " << impact[1] << endl
783// << endl << endl ;
784
785 tFemtoTrack->SetCdd(covmat[0]);
786 tFemtoTrack->SetCdz(covmat[1]);
787 tFemtoTrack->SetCzz(covmat[2]);
973a91f8 788 tFemtoTrack->SetITSchi2(tAodTrack->Chi2perNDF());
789 tFemtoTrack->SetITSncls(tAodTrack->GetITSNcls());
790 tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF());
791 tFemtoTrack->SetTPCncls(tAodTrack->GetTPCNcls());
792 tFemtoTrack->SetTPCnclsF(tAodTrack->GetTPCNcls());
76ce4b5b 793 tFemtoTrack->SetTPCsignalN(1);
794 tFemtoTrack->SetTPCsignalS(1);
795 tFemtoTrack->SetTPCsignal(tAodTrack->GetTPCsignal());
796
797// if (tPWG2AODTrack) {
798// // Copy the PWG2 specific information if it exists
799// tFemtoTrack->SetTPCClusterMap(tPWG2AODTrack->GetTPCClusterMap());
800// tFemtoTrack->SetTPCSharedMap(tPWG2AODTrack->GetTPCSharedMap());
801
802// double xtpc[3] = {0,0,0};
803// tPWG2AODTrack->GetTPCNominalEntrancePoint(xtpc);
804// tFemtoTrack->SetNominalTPCEntrancePoint(xtpc);
805// tPWG2AODTrack->GetTPCNominalExitPoint(xtpc);
806// tFemtoTrack->SetNominalTPCExitPoint(xtpc);
807// }
808// else {
809 // If not use dummy values
810 tFemtoTrack->SetTPCClusterMap(tAodTrack->GetTPCClusterMap());
811 tFemtoTrack->SetTPCSharedMap(tAodTrack->GetTPCSharedMap());
812
ce7b3d98 813
814 float globalPositionsAtRadii[9][3];
815 float bfield = 5*fMagFieldSign;
816 GetGlobalPositionAtGlobalRadiiThroughTPC(tAodTrack,bfield,globalPositionsAtRadii);
817 double tpcEntrance[3]={globalPositionsAtRadii[0][0],globalPositionsAtRadii[0][1],globalPositionsAtRadii[0][2]};
818 double **tpcPositions;
819 tpcPositions = new double*[9];
820 for(int i=0;i<9;i++)
821 tpcPositions[i] = new double[3];
822 double tpcExit[3]={globalPositionsAtRadii[8][0],globalPositionsAtRadii[8][1],globalPositionsAtRadii[8][2]};
823 for(int i=0;i<9;i++)
824 {
825 tpcPositions[i][0] = globalPositionsAtRadii[i][0];
826 tpcPositions[i][1] = globalPositionsAtRadii[i][1];
827 tpcPositions[i][2] = globalPositionsAtRadii[i][2];
828 }
829 tFemtoTrack->SetNominalTPCEntrancePoint(tpcEntrance);
830 tFemtoTrack->SetNominalTPCPoints(tpcPositions);
831 tFemtoTrack->SetNominalTPCExitPoint(tpcExit);
ba3c23a4 832
76ce4b5b 833 // }
834
835 // // cout << "Track has " << TMath::Hypot(tAodTrack->Xv(), tAodTrack->Yv()) << " " << tAodTrack->Zv() << " " << tAodTrack->GetTPCNcls() << endl;
836
837
838 int indexes[3];
839 for (int ik=0; ik<3; ik++) {
840 indexes[ik] = 0;
841 }
842 tFemtoTrack->SetKinkIndexes(indexes);
4eac0b05 843
844
845 for (int ii=0; ii<6; ii++){
846 tFemtoTrack->SetITSHitOnLayer(ii,tAodTrack->HasPointOnITSLayer(ii));
847 }
848
849
76ce4b5b 850}
851
973a91f8 852void AliFemtoEventReaderAOD::CopyAODtoFemtoV0(AliAODv0 *tAODv0, AliFemtoV0 *tFemtoV0)
853{
ce7b3d98 854 tFemtoV0->SetdecayLengthV0(tAODv0->DecayLength(fV1));
973a91f8 855 tFemtoV0->SetdecayVertexV0X(tAODv0->DecayVertexV0X());
856 tFemtoV0->SetdecayVertexV0Y(tAODv0->DecayVertexV0Y());
857 tFemtoV0->SetdecayVertexV0Z(tAODv0->DecayVertexV0Z());
858 AliFemtoThreeVector decayvertex(tAODv0->DecayVertexV0X(),tAODv0->DecayVertexV0Y(),tAODv0->DecayVertexV0Z());
859 tFemtoV0->SetdecayVertexV0(decayvertex);
860 tFemtoV0->SetdcaV0Daughters(tAODv0->DcaV0Daughters());
861 tFemtoV0->SetdcaV0ToPrimVertex(tAODv0->DcaV0ToPrimVertex());
862 tFemtoV0->SetdcaPosToPrimVertex(tAODv0->DcaPosToPrimVertex());
863 tFemtoV0->SetdcaNegToPrimVertex(tAODv0->DcaNegToPrimVertex());
864 tFemtoV0->SetmomPosX(tAODv0->MomPosX());
865 tFemtoV0->SetmomPosY(tAODv0->MomPosY());
866 tFemtoV0->SetmomPosZ(tAODv0->MomPosZ());
867 AliFemtoThreeVector mompos(tAODv0->MomPosX(),tAODv0->MomPosY(),tAODv0->MomPosZ());
868 tFemtoV0->SetmomPos(mompos);
869 tFemtoV0->SetmomNegX(tAODv0->MomNegX());
870 tFemtoV0->SetmomNegY(tAODv0->MomNegY());
871 tFemtoV0->SetmomNegZ(tAODv0->MomNegZ());
872 AliFemtoThreeVector momneg(tAODv0->MomNegX(),tAODv0->MomNegY(),tAODv0->MomNegZ());
873 tFemtoV0->SetmomNeg(momneg);
874
875 //jest cos takiego w AliFemtoV0.h czego nie ma w AliAODv0.h
876 //void SettpcHitsPos(const int& i);
877 //void SettpcHitsNeg(const int& i);
878
879 //void SetTrackTopologyMapPos(unsigned int word, const unsigned long& m);
880 //void SetTrackTopologyMapNeg(unsigned int word, const unsigned long& m);
881
882 tFemtoV0->SetmomV0X(tAODv0->MomV0X());
883 tFemtoV0->SetmomV0Y(tAODv0->MomV0Y());
884 tFemtoV0->SetmomV0Z(tAODv0->MomV0Z());
885 AliFemtoThreeVector momv0(tAODv0->MomV0X(),tAODv0->MomV0Y(),tAODv0->MomV0Z());
886 tFemtoV0->SetmomV0(momv0);
887 tFemtoV0->SetalphaV0(tAODv0->AlphaV0());
888 tFemtoV0->SetptArmV0(tAODv0->PtArmV0());
889 tFemtoV0->SeteLambda(tAODv0->ELambda());
890 tFemtoV0->SeteK0Short(tAODv0->EK0Short());
891 tFemtoV0->SetePosProton(tAODv0->EPosProton());
892 tFemtoV0->SeteNegProton(tAODv0->ENegProton());
893 tFemtoV0->SetmassLambda(tAODv0->MassLambda());
894 tFemtoV0->SetmassAntiLambda(tAODv0->MassAntiLambda());
895 tFemtoV0->SetmassK0Short(tAODv0->MassK0Short());
896 tFemtoV0->SetrapLambda(tAODv0->RapLambda());
897 tFemtoV0->SetrapK0Short(tAODv0->RapK0Short());
898
899 //void SetcTauLambda( float x);
900 //void SetcTauK0Short( float x);
901
ce7b3d98 902 //tFemtoV0->SetptV0(::sqrt(tAODv0->Pt2V0())); //!
903 tFemtoV0->SetptV0(tAODv0->Pt());
973a91f8 904 tFemtoV0->SetptotV0(::sqrt(tAODv0->Ptot2V0()));
ce7b3d98 905 //tFemtoV0->SetptPos(::sqrt(tAODv0->MomPosX()*tAODv0->MomPosX()+tAODv0->MomPosY()*tAODv0->MomPosY()));
906 //tFemtoV0->SetptotPos(::sqrt(tAODv0->Ptot2Pos()));
907 //tFemtoV0->SetptNeg(::sqrt(tAODv0->MomNegX()*tAODv0->MomNegX()+tAODv0->MomNegY()*tAODv0->MomNegY()));
908 //tFemtoV0->SetptotNeg(::sqrt(tAODv0->Ptot2Neg()));
973a91f8 909
910 tFemtoV0->SetidNeg(tAODv0->GetNegID());
911 //cout<<"tAODv0->GetNegID(): "<<tAODv0->GetNegID()<<endl;
912 //cout<<"tFemtoV0->IdNeg(): "<<tFemtoV0->IdNeg()<<endl;
913 tFemtoV0->SetidPos(tAODv0->GetPosID());
914
915 tFemtoV0->SetEtaV0(tAODv0->Eta());
916 tFemtoV0->SetPhiV0(tAODv0->Phi());
917 tFemtoV0->SetCosPointingAngle(tAODv0->CosPointingAngle(fV1));
918 //tFemtoV0->SetYV0(tAODv0->Y());
919
920 //void SetdedxNeg(float x);
921 //void SeterrdedxNeg(float x);//Gael 04Fev2002
922 //void SetlendedxNeg(float x);//Gael 04Fev2002
923 //void SetdedxPos(float x);
924 //void SeterrdedxPos(float x);//Gael 04Fev2002
925 //void SetlendedxPos(float x);//Gael 04Fev2002
926
ce7b3d98 927 //tFemtoV0->SetEtaPos(tAODv0->PseudoRapPos());
928 //tFemtoV0->SetEtaNeg(tAODv0->PseudoRapNeg());
973a91f8 929
930 AliAODTrack *trackpos = (AliAODTrack*)tAODv0->GetDaughter(0);
931 AliAODTrack *trackneg = (AliAODTrack*)tAODv0->GetDaughter(1);
932
933 if(trackpos && trackneg)
934 {
ce7b3d98 935 tFemtoV0->SetEtaPos(trackpos->Eta());
936 tFemtoV0->SetEtaNeg(trackneg->Eta());
937 tFemtoV0->SetptotPos(tAODv0->PProng(0));
938 tFemtoV0->SetptotNeg(tAODv0->PProng(1));
939 tFemtoV0->SetptPos(trackpos->Pt());
940 tFemtoV0->SetptNeg(trackneg->Pt());
941
973a91f8 942 //tFemtoV0->SetEtaPos(trackpos->Eta()); //tAODv0->PseudoRapPos()
943 //tFemtoV0->SetEtaNeg(trackneg->Eta()); //tAODv0->PseudoRapNeg()
944 tFemtoV0->SetTPCNclsPos(trackpos->GetTPCNcls());
945 tFemtoV0->SetTPCNclsNeg(trackneg->GetTPCNcls());
946 tFemtoV0->SetTPCclustersPos(trackpos->GetTPCClusterMap());
947 tFemtoV0->SetTPCclustersNeg(trackneg->GetTPCClusterMap());
948 tFemtoV0->SetTPCsharingPos(trackpos->GetTPCSharedMap());
949 tFemtoV0->SetTPCsharingNeg(trackneg->GetTPCSharedMap());
950 tFemtoV0->SetNdofPos(trackpos->Chi2perNDF());
951 tFemtoV0->SetNdofNeg(trackneg->Chi2perNDF());
952 tFemtoV0->SetStatusPos(trackpos->GetStatus());
953 tFemtoV0->SetStatusNeg(trackneg->GetStatus());
954
955 tFemtoV0->SetPosNSigmaTPCK(fAODpidUtil->NumberOfSigmasTPC(trackpos,AliPID::kKaon));
956 tFemtoV0->SetNegNSigmaTPCK(fAODpidUtil->NumberOfSigmasTPC(trackneg,AliPID::kKaon));
957 tFemtoV0->SetPosNSigmaTPCP(fAODpidUtil->NumberOfSigmasTPC(trackpos,AliPID::kProton));
958 tFemtoV0->SetNegNSigmaTPCP(fAODpidUtil->NumberOfSigmasTPC(trackneg,AliPID::kProton));
959 tFemtoV0->SetPosNSigmaTPCPi(fAODpidUtil->NumberOfSigmasTPC(trackpos,AliPID::kPion));
960 tFemtoV0->SetNegNSigmaTPCPi(fAODpidUtil->NumberOfSigmasTPC(trackneg,AliPID::kPion));
961
ce7b3d98 962
963 float bfield = 5*fMagFieldSign;
964 float globalPositionsAtRadiiPos[9][3];
965 GetGlobalPositionAtGlobalRadiiThroughTPC(trackpos,bfield,globalPositionsAtRadiiPos);
966 double tpcEntrancePos[3]={globalPositionsAtRadiiPos[0][0],globalPositionsAtRadiiPos[0][1],globalPositionsAtRadiiPos[0][2]};
967 double tpcExitPos[3]={globalPositionsAtRadiiPos[8][0],globalPositionsAtRadiiPos[8][1],globalPositionsAtRadiiPos[8][2]};
968
969 float globalPositionsAtRadiiNeg[9][3];
970 GetGlobalPositionAtGlobalRadiiThroughTPC(trackneg,bfield,globalPositionsAtRadiiNeg);
971 double tpcEntranceNeg[3]={globalPositionsAtRadiiNeg[0][0],globalPositionsAtRadiiNeg[0][1],globalPositionsAtRadiiNeg[0][2]};
972 double tpcExitNeg[3]={globalPositionsAtRadiiNeg[8][0],globalPositionsAtRadiiNeg[8][1],globalPositionsAtRadiiNeg[8][2]};
973
974 AliFemtoThreeVector tmpVec;
975 tmpVec.SetX(tpcEntrancePos[0]); tmpVec.SetX(tpcEntrancePos[1]); tmpVec.SetX(tpcEntrancePos[2]);
ba3c23a4 976 tFemtoV0->SetNominalTpcEntrancePointPos(tmpVec);
ce7b3d98 977
978 tmpVec.SetX(tpcExitPos[0]); tmpVec.SetX(tpcExitPos[1]); tmpVec.SetX(tpcExitPos[2]);
ba3c23a4 979 tFemtoV0->SetNominalTpcExitPointPos(tmpVec);
ce7b3d98 980
981 tmpVec.SetX(tpcEntranceNeg[0]); tmpVec.SetX(tpcEntranceNeg[1]); tmpVec.SetX(tpcEntranceNeg[2]);
982 tFemtoV0->SetNominalTpcEntrancePointNeg(tmpVec);
983
984 tmpVec.SetX(tpcExitNeg[0]); tmpVec.SetX(tpcExitNeg[1]); tmpVec.SetX(tpcExitNeg[2]);
ba3c23a4 985 tFemtoV0->SetNominalTpcExitPointNeg(tmpVec);
973a91f8 986
ce7b3d98 987 AliFemtoThreeVector vecTpcPos[9];
988 AliFemtoThreeVector vecTpcNeg[9];
989 for(int i=0;i<9;i++)
990 {
991 vecTpcPos[i].SetX(globalPositionsAtRadiiPos[i][0]); vecTpcPos[i].SetY(globalPositionsAtRadiiPos[i][1]); vecTpcPos[i].SetZ(globalPositionsAtRadiiPos[i][2]);
992 vecTpcNeg[i].SetX(globalPositionsAtRadiiNeg[i][0]); vecTpcNeg[i].SetY(globalPositionsAtRadiiNeg[i][1]); vecTpcNeg[i].SetZ(globalPositionsAtRadiiNeg[i][2]);
993 }
994 tFemtoV0->SetNominalTpcPointPos(vecTpcPos);
995 tFemtoV0->SetNominalTpcPointNeg(vecTpcNeg);
996
997 tFemtoV0->SetTPCMomentumPos(trackpos->GetTPCmomentum());
998 tFemtoV0->SetTPCMomentumNeg(trackneg->GetTPCmomentum());
999
1000 tFemtoV0->SetdedxPos(trackpos->GetTPCsignal());
1001 tFemtoV0->SetdedxNeg(trackneg->GetTPCsignal());
1002
63e066f7 1003 if((tFemtoV0->StatusPos()&AliESDtrack::kTOFpid)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTOFout)==0 )
973a91f8 1004 {
63e066f7 1005 if((tFemtoV0->StatusNeg()&AliESDtrack::kTOFpid)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTOFout)==0 )
973a91f8 1006 {
ce7b3d98 1007 tFemtoV0->SetPosNSigmaTOFK(-1000);
1008 tFemtoV0->SetNegNSigmaTOFK(-1000);
1009 tFemtoV0->SetPosNSigmaTOFP(-1000);
1010 tFemtoV0->SetNegNSigmaTOFP(-1000);
1011 tFemtoV0->SetPosNSigmaTOFPi(-1000);
1012 tFemtoV0->SetNegNSigmaTOFPi(-1000);
1013
1014 tFemtoV0->SetTOFProtonTimePos(-1000);
1015 tFemtoV0->SetTOFPionTimePos(-1000);
1016 tFemtoV0->SetTOFKaonTimePos(-1000);
1017 tFemtoV0->SetTOFProtonTimeNeg(-1000);
1018 tFemtoV0->SetTOFPionTimeNeg(-1000);
1019 tFemtoV0->SetTOFKaonTimeNeg(-1000);
973a91f8 1020 }
1021 }
1022 else
1023 {
1024 tFemtoV0->SetPosNSigmaTOFK(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kKaon));
1025 tFemtoV0->SetNegNSigmaTOFK(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kKaon));
1026 tFemtoV0->SetPosNSigmaTOFP(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kProton));
1027 tFemtoV0->SetNegNSigmaTOFP(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kProton));
1028 tFemtoV0->SetPosNSigmaTOFPi(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kPion));
1029 tFemtoV0->SetNegNSigmaTOFPi(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kPion));
ce7b3d98 1030
1031 double TOFSignalPos = trackpos->GetTOFsignal();
1032 double TOFSignalNeg = trackneg->GetTOFsignal();
1033 double pidPos[5];
1034 double pidNeg[5];
1035 trackpos->GetIntegratedTimes(pidPos);
1036 trackneg->GetIntegratedTimes(pidNeg);
1037
1038 tFemtoV0->SetTOFPionTimePos(TOFSignalPos-pidPos[2]);
1039 tFemtoV0->SetTOFKaonTimePos(TOFSignalPos-pidPos[3]);
1040 tFemtoV0->SetTOFProtonTimePos(TOFSignalPos-pidPos[4]);
1041 tFemtoV0->SetTOFPionTimeNeg(TOFSignalNeg-pidNeg[2]);
1042 tFemtoV0->SetTOFKaonTimeNeg(TOFSignalNeg-pidNeg[3]);
1043 tFemtoV0->SetTOFProtonTimeNeg(TOFSignalNeg-pidNeg[4]);
973a91f8 1044 }
1045 }
1046 else
1047 {
1048 tFemtoV0->SetStatusPos(999);
1049 tFemtoV0->SetStatusNeg(999);
1050 }
1051 tFemtoV0->SetOnFlyStatusV0(tAODv0->GetOnFlyStatus());
1052}
1053
76ce4b5b 1054void AliFemtoEventReaderAOD::SetFilterBit(UInt_t ibit)
1055{
1056 fFilterBit = (1 << (ibit));
1057}
1058
64536eaf 1059
1060void AliFemtoEventReaderAOD::SetFilterMask(int ibit)
1061{
1062 fFilterMask = ibit;
1063}
1064
76ce4b5b 1065void AliFemtoEventReaderAOD::SetReadMC(unsigned char a)
1066{
1067 fReadMC = a;
1068}
1069
973a91f8 1070
1071void AliFemtoEventReaderAOD::SetReadV0(unsigned char a)
1072{
1073 fReadV0 = a;
1074}
1075
1777446b 1076void AliFemtoEventReaderAOD::SetUseMultiplicity(EstEventMult aType)
1077{
1078 fEstEventMult = aType;
1079}
1080
76ce4b5b 1081AliAODMCParticle* AliFemtoEventReaderAOD::GetParticleWithLabel(TClonesArray *mcP, Int_t aLabel)
1082{
1083 if (aLabel < 0) return 0;
1084 AliAODMCParticle *aodP;
1085 Int_t posstack = 0;
1086 if (aLabel > mcP->GetEntries())
1087 posstack = mcP->GetEntries();
1088 else
1089 posstack = aLabel;
1090
1091 aodP = (AliAODMCParticle *) mcP->At(posstack);
1092 if (aodP->GetLabel() > posstack) {
1093 do {
1094 aodP = (AliAODMCParticle *) mcP->At(posstack);
1095 if (aodP->GetLabel() == aLabel) return aodP;
1096 posstack--;
1097 }
1098 while (posstack > 0);
1099 }
1100 else {
1101 do {
1102 aodP = (AliAODMCParticle *) mcP->At(posstack);
1103 if (aodP->GetLabel() == aLabel) return aodP;
1104 posstack++;
1105 }
1106 while (posstack < mcP->GetEntries());
1107 }
1108
1109 return 0;
1110}
1111
1112void AliFemtoEventReaderAOD::CopyPIDtoFemtoTrack(AliAODTrack *tAodTrack,
1113 AliFemtoTrack *tFemtoTrack)
1114{
41cb6c1e 1115
6691ea4f 1116 // copying DCA information (taking it from global tracks gives better resolution than from TPC-only)
41cb6c1e 1117
6691ea4f 1118 double impact[2];
1119 double covimpact[3];
1120
1121 if (!tAodTrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
1122 //cout << "sth went wrong with dca propagation" << endl;
1123 tFemtoTrack->SetImpactD(-1000.0);
1124 tFemtoTrack->SetImpactZ(-1000.0);
1125
1126 }
1127 else {
1128 tFemtoTrack->SetImpactD(impact[0]);
1129 tFemtoTrack->SetImpactZ(impact[1]);
1130 }
41cb6c1e 1131
76ce4b5b 1132 double aodpid[10];
1133 tAodTrack->GetPID(aodpid);
1134 tFemtoTrack->SetPidProbElectron(aodpid[0]);
1135 tFemtoTrack->SetPidProbMuon(aodpid[1]);
1136 tFemtoTrack->SetPidProbPion(aodpid[2]);
1137 tFemtoTrack->SetPidProbKaon(aodpid[3]);
1138 tFemtoTrack->SetPidProbProton(aodpid[4]);
1139
1140 aodpid[0] = -100000.0;
1141 aodpid[1] = -100000.0;
1142 aodpid[2] = -100000.0;
1143 aodpid[3] = -100000.0;
1144 aodpid[4] = -100000.0;
6691ea4f 1145
76ce4b5b 1146 double tTOF = 0.0;
1147
4eac0b05 1148 //what is that code? for what do we need it? nsigma values are not enough?
1149 if (tAodTrack->GetStatus() & AliESDtrack::kTOFpid) { //AliESDtrack::kTOFpid=0x8000
1150 tTOF = tAodTrack->GetTOFsignal();
1151 tAodTrack->GetIntegratedTimes(aodpid);
1777446b 1152
4eac0b05 1153 tTOF -= fAODpidUtil->GetTOFResponse().GetStartTime(tAodTrack->P());
1154 }
76ce4b5b 1155
4eac0b05 1156 tFemtoTrack->SetTofExpectedTimes(tTOF-aodpid[2], tTOF-aodpid[3], tTOF-aodpid[4]);
76ce4b5b 1157
1158 ////// TPC ////////////////////////////////////////////
1159
1160 float nsigmaTPCK=-1000.;
1161 float nsigmaTPCPi=-1000.;
1162 float nsigmaTPCP=-1000.;
1163
1164 // cout<<"in reader fESDpid"<<fESDpid<<endl;
1165
1166 if (tAodTrack->IsOn(AliESDtrack::kTPCpid)){ //AliESDtrack::kTPCpid=0x0080
1167 nsigmaTPCK = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kKaon);
1168 nsigmaTPCPi = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kPion);
1169 nsigmaTPCP = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kProton);
1170 }
1171
1172 tFemtoTrack->SetNSigmaTPCPi(nsigmaTPCPi);
1173 tFemtoTrack->SetNSigmaTPCK(nsigmaTPCK);
1174 tFemtoTrack->SetNSigmaTPCP(nsigmaTPCP);
1175
1176 tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF());
1177 tFemtoTrack->SetTPCncls(tAodTrack->GetTPCNcls());
1178 tFemtoTrack->SetTPCnclsF(tAodTrack->GetTPCNcls());
1179
1180 tFemtoTrack->SetTPCsignalN(1);
1181 tFemtoTrack->SetTPCsignalS(1);
1182 tFemtoTrack->SetTPCsignal(tAodTrack->GetTPCsignal());
1183
1184 ///////TOF//////////////////////
1185
1186 float vp=-1000.;
1187 float nsigmaTOFPi=-1000.;
1188 float nsigmaTOFK=-1000.;
1189 float nsigmaTOFP=-1000.;
1190
1191 if ((tAodTrack->GetStatus() & AliESDtrack::kTOFpid) && //AliESDtrack::kTOFpid=0x8000
1192 (tAodTrack->GetStatus() & AliESDtrack::kTOFout) && //AliESDtrack::kTOFout=0x2000
63e066f7 1193 (tAodTrack->GetStatus() & AliESDtrack::kTIME) ) //AliESDtrack::kTIME=0x80000000
76ce4b5b 1194 {
1195 if(tAodTrack->IsOn(AliESDtrack::kTOFpid)) //AliESDtrack::kTOFpid=0x8000
1196 {
1197
1198 nsigmaTOFPi = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kPion);
1199 nsigmaTOFK = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kKaon);
1200 nsigmaTOFP = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kProton);
1201
1202 Double_t len=200;// esdtrack->GetIntegratedLength(); !!!!!
1203 Double_t tof=tAodTrack->GetTOFsignal();
1204 if(tof > 0.) vp=len/tof/0.03;
1205 }
1206 }
1207 tFemtoTrack->SetVTOF(vp);
1208 tFemtoTrack->SetNSigmaTOFPi(nsigmaTOFPi);
1209 tFemtoTrack->SetNSigmaTOFK(nsigmaTOFK);
1210 tFemtoTrack->SetNSigmaTOFP(nsigmaTOFP);
1211
1212
1213 //////////////////////////////////////
1214
1215}
1216
1217void AliFemtoEventReaderAOD::SetCentralityPreSelection(double min, double max)
1218{
1219 fCentRange[0] = min; fCentRange[1] = max;
1220 fUsePreCent = 1;
1777446b 1221 fEstEventMult = kCentrality;
76ce4b5b 1222}
1223
973a91f8 1224void AliFemtoEventReaderAOD::SetNoCentrality(bool anocent)
1225{
1777446b 1226 if(anocent==false) {fEstEventMult=kCentrality;}
1227 else {fEstEventMult=kReference; fUsePreCent = 0; }
973a91f8 1228}
76ce4b5b 1229
1230void AliFemtoEventReaderAOD::SetAODpidUtil(AliAODpidUtil *aAODpidUtil)
1231{
1232 fAODpidUtil = aAODpidUtil;
1233 // printf("fAODpidUtil: %x\n",fAODpidUtil);
1234}
1235
1777446b 1236void AliFemtoEventReaderAOD::SetAODheader(AliAODHeader *aAODheader)
1237{
1238 fAODheader = aAODheader;
1239}
1240
76ce4b5b 1241
ba3c23a4 1242void AliFemtoEventReaderAOD::SetMagneticFieldSign(int s)
1243{
1244 if(s>0)
1245 fMagFieldSign = 1;
1246 else if(s<0)
1247 fMagFieldSign = -1;
1248 else
1249 fMagFieldSign = 0;
1250}
1251
5e2038a9 1252void AliFemtoEventReaderAOD::SetEPVZERO(Bool_t iepvz)
1253{
1254 fisEPVZ = iepvz;
1255}
76ce4b5b 1256
ce7b3d98 1257void AliFemtoEventReaderAOD::GetGlobalPositionAtGlobalRadiiThroughTPC(AliAODTrack *track, Float_t bfield, Float_t globalPositionsAtRadii[9][3])
1258{
1259 // Gets the global position of the track at nine different radii in the TPC
1260 // track is the track you want to propagate
1261 // bfield is the magnetic field of your event
1262 // globalPositionsAtRadii is the array of global positions in the radii and xyz
1263
1264 // Initialize the array to something indicating there was no propagation
1265 for(Int_t i=0;i<9;i++){
1266 for(Int_t j=0;j<3;j++){
1267 globalPositionsAtRadii[i][j]=-9999.;
1268 }
1269 }
1270
1271 // Make a copy of the track to not change parameters of the track
1272 AliExternalTrackParam etp; etp.CopyFromVTrack(track);
1273 //printf("\nAfter CopyFromVTrack\n");
1274 //etp.Print();
1275
1276 // The global position of the the track
1277 Double_t xyz[3]={-9999.,-9999.,-9999.};
1278
1279 // Counter for which radius we want
1280 Int_t iR=0;
1281 // The radii at which we get the global positions
1282 // IROC (OROC) from 84.1 cm to 132.1 cm (134.6 cm to 246.6 cm)
1283 Float_t Rwanted[9]={85.,105.,125.,145.,165.,185.,205.,225.,245.};
1284 // The global radius we are at
1285 Float_t globalRadius=0;
1286
1287 // Propagation is done in local x of the track
1288 for (Float_t x = etp.GetX();x<247.;x+=1.){ // GetX returns local coordinates
1289 // Starts at the tracks fX and goes outwards. x = 245 is the outer radial limit
1290 // of the TPC when the track is straight, i.e. has inifinite pt and doesn't get bent.
1291 // If the track's momentum is smaller than infinite, it will develop a y-component, which
1292 // adds to the global radius
1293
1294 // Stop if the propagation was not succesful. This can happen for low pt tracks
1295 // that don't reach outer radii
1296 if(!etp.PropagateTo(x,bfield))break;
1297 etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1298 globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
1299
1300 // Roughly reached the radius we want
1301 if(globalRadius > Rwanted[iR]){
1302
1303 // Bigger loop has bad precision, we're nearly one centimeter too far, go back in small steps.
1304 while (globalRadius>Rwanted[iR]){
1305 x-=.1;
1306 // printf("propagating to x %5.2f\n",x);
1307 if(!etp.PropagateTo(x,bfield))break;
1308 etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1309 globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
1310 }
1311 //printf("At Radius:%05.2f (local x %5.2f). Setting position to x %4.1f y %4.1f z %4.1f\n",globalRadius,x,xyz[0],xyz[1],xyz[2]);
1312 globalPositionsAtRadii[iR][0]=xyz[0];
1313 globalPositionsAtRadii[iR][1]=xyz[1];
1314 globalPositionsAtRadii[iR][2]=xyz[2];
1315 // Indicate we want the next radius
1316 iR+=1;
1317 }
1318 if(iR>=8){
1319 // TPC edge reached
1320 return;
1321 }
1322 }
1323}
76ce4b5b 1324
1325