]> 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
Fix in order to allow to use the flatenning for stripped AOD (Julian)
[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
2bcd4b96 650 if (cent) {
651 tEvent->SetCentralityV0(cent->GetCentralityPercentile("V0M"));
652 tEvent->SetCentralityZNA(cent->GetCentralityPercentile("ZNA"));
653 tEvent->SetCentralityCL1(cent->GetCentralityPercentile("CL1"));
654 // tEvent->SetCentralityFMD(cent->GetCentralityPercentile("FMD"));
607aa748 655 tEvent->SetCentralityTrk(cent->GetCentralityPercentile("TRK"));
656 tEvent->SetCentralityCND(cent->GetCentralityPercentile("CND"));
2bcd4b96 657 }
658
1777446b 659 if (fEstEventMult==kCentrality) {
2bcd4b96 660 //AliCentrality *cent = fEvent->GetCentrality();
4eac0b05 661 //cout<<"AliFemtoEventReaderAOD:"<<lrint(10*cent->GetCentralityPercentile("V0M"))<<endl;
973a91f8 662 if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("V0M")));
663 // if (cent) tEvent->SetNormalizedMult((int) cent->GetCentralityPercentile("V0M"));
c863b24a 664 }
665 else if (fEstEventMult==kCentralityZNA) {
c863b24a 666 if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("ZNA")));
c863b24a 667 }
668 else if (fEstEventMult==kCentralityCL1) {
c863b24a 669 if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("CL1")));
76ce4b5b 670 }
607aa748 671 else if (fEstEventMult==kCentralityTRK) {
672 if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("TRK")));
673 }
674 else if (fEstEventMult==kCentralityCND) {
675 if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("CND")));
676 }
1777446b 677 else if(fEstEventMult==kGlobalCount){
ba3c23a4 678 tEvent->SetNormalizedMult(tNormMult); //particles counted in the loop, trying to reproduce GetReferenceMultiplicity. If better (default) method appears it should be changed
679 }
1777446b 680 else if(fEstEventMult==kReference)
681 {
682 tEvent->SetNormalizedMult(fAODheader->GetRefMultiplicity());
683 }
684 else if(fEstEventMult==kTPCOnlyRef)
685 {
686 tEvent->SetNormalizedMult(fAODheader->GetTPConlyRefMultiplicity());
687 }
1938eadf 688 else if(fEstEventMult == kVZERO)
689 {
690 Float_t multV0 = 0;
691 for (Int_t i=0; i<64; i++)
692 multV0 += fEvent->GetVZEROData()->GetMultiplicity(i);
693 tEvent->SetNormalizedMult(multV0);
694 }
76ce4b5b 695
696 if (mcP) delete [] motherids;
697
5e2038a9 698 // cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
973a91f8 699
700 if(fReadV0)
701 {
ce7b3d98 702 int count_pass = 0;
973a91f8 703 for (Int_t i = 0; i < fEvent->GetNumberOfV0s(); i++) {
704 AliAODv0* aodv0 = fEvent->GetV0(i);
705 if (!aodv0) continue;
706 if(aodv0->GetNDaughters()>2) continue;
707 if(aodv0->GetNProngs()>2) continue;
708 if(aodv0->GetCharge()!=0) continue;
709 if(aodv0->ChargeProng(0)==aodv0->ChargeProng(1)) continue;
ce7b3d98 710 if(aodv0->CosPointingAngle(fV1)<0.998) continue;
973a91f8 711 AliFemtoV0* trackCopyV0 = new AliFemtoV0();
ce7b3d98 712 count_pass++;
973a91f8 713 CopyAODtoFemtoV0(aodv0, trackCopyV0);
714 tEvent->V0Collection()->push_back(trackCopyV0);
715 //cout<<"Pushback v0 to v0collection"<<endl;
716 }
717 }
718
76ce4b5b 719}
720
721void AliFemtoEventReaderAOD::CopyAODtoFemtoTrack(AliAODTrack *tAodTrack,
722 AliFemtoTrack *tFemtoTrack
723 // AliPWG2AODTrack *tPWG2AODTrack
724 )
725{
726 // Copy the track information from the AOD into the internal AliFemtoTrack
727 // If it exists, use the additional information from the PWG2 AOD
728
729 // Primary Vertex position
973a91f8 730
76ce4b5b 731 fEvent->GetPrimaryVertex()->GetPosition(fV1);
732 // fEvent->GetPrimaryVertex()->GetXYZ(fV1);
733
734 tFemtoTrack->SetCharge(tAodTrack->Charge());
735
736 double pxyz[3];
737 tAodTrack->PxPyPz(pxyz);//reading noconstrained momentum
738 AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
739 tFemtoTrack->SetP(v);//setting momentum
740 tFemtoTrack->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
741 const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
742 //setting track helix
743 const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
744 AliFmPhysicalHelixD helix(ktP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(tFemtoTrack->Charge()));
745 tFemtoTrack->SetHelix(helix);
746
747 // Flags
748 tFemtoTrack->SetTrackId(tAodTrack->GetID());
749 tFemtoTrack->SetFlags(tAodTrack->GetFlags());
750 tFemtoTrack->SetLabel(tAodTrack->GetLabel());
751
752 // Track quality information
753 float covmat[6];
754 tAodTrack->GetCovMatrix(covmat);
755
41cb6c1e 756 // ! DCA information is done in CopyPIDtoFemtoTrack()
76ce4b5b 757
41cb6c1e 758 // double impact[2];
759 // double covimpact[3];
76ce4b5b 760
41cb6c1e 761 // if (!tAodTrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
762 // //cout << "sth went wrong with dca propagation" << endl;
763 // tFemtoTrack->SetImpactD(-1000.0);
764 // tFemtoTrack->SetImpactZ(-1000.0);
765
766 // }
767 // else {
768 // tFemtoTrack->SetImpactD(impact[0]);
769 // tFemtoTrack->SetImpactZ(impact[1]+fV1[2]);
770 // }
76ce4b5b 771
772 // if (TMath::Abs(tAodTrack->Xv()) > 0.00000000001)
773 // tFemtoTrack->SetImpactD(TMath::Hypot(tAodTrack->Xv(), tAodTrack->Yv())*(tAodTrack->Xv()/TMath::Abs(tAodTrack->Xv())));
774 // else
775 // tFemtoTrack->SetImpactD(0.0);
776 // tFemtoTrack->SetImpactD(tAodTrack->DCA());
777
778 // tFemtoTrack->SetImpactZ(tAodTrack->ZAtDCA());
779
780
781 // tFemtoTrack->SetImpactD(TMath::Hypot(tAodTrack->Xv() - fV1[0], tAodTrack->Yv() - fV1[1]));
782 // tFemtoTrack->SetImpactZ(tAodTrack->Zv() - fV1[2]);
783
784
785 // cout
786 // << "dca" << TMath::Hypot(tAodTrack->Xv() - fV1[0], tAodTrack->Yv() - fV1[1])
787 // << "xv - fv10 = "<< tAodTrack->Xv() - fV1[0]
788 // << tAodTrack->Yv() - fV1[1]
789// << "xv = " << tAodTrack->Xv() << endl
790// << "fv1[0] = " << fV1[0] << endl
791// << "yv = " << tAodTrack->Yv() << endl
792// << "fv1[1] = " << fV1[1] << endl
793// << "zv = " << tAodTrack->Zv() << endl
794// << "fv1[2] = " << fV1[2] << endl
795// << "impact[0] = " << impact[0] << endl
796// << "impact[1] = " << impact[1] << endl
797// << endl << endl ;
798
799 tFemtoTrack->SetCdd(covmat[0]);
800 tFemtoTrack->SetCdz(covmat[1]);
801 tFemtoTrack->SetCzz(covmat[2]);
973a91f8 802 tFemtoTrack->SetITSchi2(tAodTrack->Chi2perNDF());
803 tFemtoTrack->SetITSncls(tAodTrack->GetITSNcls());
804 tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF());
805 tFemtoTrack->SetTPCncls(tAodTrack->GetTPCNcls());
806 tFemtoTrack->SetTPCnclsF(tAodTrack->GetTPCNcls());
76ce4b5b 807 tFemtoTrack->SetTPCsignalN(1);
808 tFemtoTrack->SetTPCsignalS(1);
809 tFemtoTrack->SetTPCsignal(tAodTrack->GetTPCsignal());
810
811// if (tPWG2AODTrack) {
812// // Copy the PWG2 specific information if it exists
813// tFemtoTrack->SetTPCClusterMap(tPWG2AODTrack->GetTPCClusterMap());
814// tFemtoTrack->SetTPCSharedMap(tPWG2AODTrack->GetTPCSharedMap());
815
816// double xtpc[3] = {0,0,0};
817// tPWG2AODTrack->GetTPCNominalEntrancePoint(xtpc);
818// tFemtoTrack->SetNominalTPCEntrancePoint(xtpc);
819// tPWG2AODTrack->GetTPCNominalExitPoint(xtpc);
820// tFemtoTrack->SetNominalTPCExitPoint(xtpc);
821// }
822// else {
823 // If not use dummy values
824 tFemtoTrack->SetTPCClusterMap(tAodTrack->GetTPCClusterMap());
825 tFemtoTrack->SetTPCSharedMap(tAodTrack->GetTPCSharedMap());
826
ce7b3d98 827
828 float globalPositionsAtRadii[9][3];
829 float bfield = 5*fMagFieldSign;
830 GetGlobalPositionAtGlobalRadiiThroughTPC(tAodTrack,bfield,globalPositionsAtRadii);
831 double tpcEntrance[3]={globalPositionsAtRadii[0][0],globalPositionsAtRadii[0][1],globalPositionsAtRadii[0][2]};
832 double **tpcPositions;
833 tpcPositions = new double*[9];
834 for(int i=0;i<9;i++)
835 tpcPositions[i] = new double[3];
836 double tpcExit[3]={globalPositionsAtRadii[8][0],globalPositionsAtRadii[8][1],globalPositionsAtRadii[8][2]};
837 for(int i=0;i<9;i++)
838 {
839 tpcPositions[i][0] = globalPositionsAtRadii[i][0];
840 tpcPositions[i][1] = globalPositionsAtRadii[i][1];
841 tpcPositions[i][2] = globalPositionsAtRadii[i][2];
842 }
843 tFemtoTrack->SetNominalTPCEntrancePoint(tpcEntrance);
844 tFemtoTrack->SetNominalTPCPoints(tpcPositions);
845 tFemtoTrack->SetNominalTPCExitPoint(tpcExit);
ba3c23a4 846
76ce4b5b 847 // }
848
849 // // cout << "Track has " << TMath::Hypot(tAodTrack->Xv(), tAodTrack->Yv()) << " " << tAodTrack->Zv() << " " << tAodTrack->GetTPCNcls() << endl;
850
851
852 int indexes[3];
853 for (int ik=0; ik<3; ik++) {
854 indexes[ik] = 0;
855 }
856 tFemtoTrack->SetKinkIndexes(indexes);
4eac0b05 857
858
859 for (int ii=0; ii<6; ii++){
860 tFemtoTrack->SetITSHitOnLayer(ii,tAodTrack->HasPointOnITSLayer(ii));
861 }
862
863
76ce4b5b 864}
865
973a91f8 866void AliFemtoEventReaderAOD::CopyAODtoFemtoV0(AliAODv0 *tAODv0, AliFemtoV0 *tFemtoV0)
867{
ce7b3d98 868 tFemtoV0->SetdecayLengthV0(tAODv0->DecayLength(fV1));
973a91f8 869 tFemtoV0->SetdecayVertexV0X(tAODv0->DecayVertexV0X());
870 tFemtoV0->SetdecayVertexV0Y(tAODv0->DecayVertexV0Y());
871 tFemtoV0->SetdecayVertexV0Z(tAODv0->DecayVertexV0Z());
872 AliFemtoThreeVector decayvertex(tAODv0->DecayVertexV0X(),tAODv0->DecayVertexV0Y(),tAODv0->DecayVertexV0Z());
873 tFemtoV0->SetdecayVertexV0(decayvertex);
874 tFemtoV0->SetdcaV0Daughters(tAODv0->DcaV0Daughters());
875 tFemtoV0->SetdcaV0ToPrimVertex(tAODv0->DcaV0ToPrimVertex());
876 tFemtoV0->SetdcaPosToPrimVertex(tAODv0->DcaPosToPrimVertex());
877 tFemtoV0->SetdcaNegToPrimVertex(tAODv0->DcaNegToPrimVertex());
878 tFemtoV0->SetmomPosX(tAODv0->MomPosX());
879 tFemtoV0->SetmomPosY(tAODv0->MomPosY());
880 tFemtoV0->SetmomPosZ(tAODv0->MomPosZ());
881 AliFemtoThreeVector mompos(tAODv0->MomPosX(),tAODv0->MomPosY(),tAODv0->MomPosZ());
882 tFemtoV0->SetmomPos(mompos);
883 tFemtoV0->SetmomNegX(tAODv0->MomNegX());
884 tFemtoV0->SetmomNegY(tAODv0->MomNegY());
885 tFemtoV0->SetmomNegZ(tAODv0->MomNegZ());
886 AliFemtoThreeVector momneg(tAODv0->MomNegX(),tAODv0->MomNegY(),tAODv0->MomNegZ());
887 tFemtoV0->SetmomNeg(momneg);
888
889 //jest cos takiego w AliFemtoV0.h czego nie ma w AliAODv0.h
890 //void SettpcHitsPos(const int& i);
891 //void SettpcHitsNeg(const int& i);
892
893 //void SetTrackTopologyMapPos(unsigned int word, const unsigned long& m);
894 //void SetTrackTopologyMapNeg(unsigned int word, const unsigned long& m);
895
896 tFemtoV0->SetmomV0X(tAODv0->MomV0X());
897 tFemtoV0->SetmomV0Y(tAODv0->MomV0Y());
898 tFemtoV0->SetmomV0Z(tAODv0->MomV0Z());
899 AliFemtoThreeVector momv0(tAODv0->MomV0X(),tAODv0->MomV0Y(),tAODv0->MomV0Z());
900 tFemtoV0->SetmomV0(momv0);
901 tFemtoV0->SetalphaV0(tAODv0->AlphaV0());
902 tFemtoV0->SetptArmV0(tAODv0->PtArmV0());
903 tFemtoV0->SeteLambda(tAODv0->ELambda());
904 tFemtoV0->SeteK0Short(tAODv0->EK0Short());
905 tFemtoV0->SetePosProton(tAODv0->EPosProton());
906 tFemtoV0->SeteNegProton(tAODv0->ENegProton());
907 tFemtoV0->SetmassLambda(tAODv0->MassLambda());
908 tFemtoV0->SetmassAntiLambda(tAODv0->MassAntiLambda());
909 tFemtoV0->SetmassK0Short(tAODv0->MassK0Short());
910 tFemtoV0->SetrapLambda(tAODv0->RapLambda());
911 tFemtoV0->SetrapK0Short(tAODv0->RapK0Short());
912
913 //void SetcTauLambda( float x);
914 //void SetcTauK0Short( float x);
915
ce7b3d98 916 //tFemtoV0->SetptV0(::sqrt(tAODv0->Pt2V0())); //!
917 tFemtoV0->SetptV0(tAODv0->Pt());
973a91f8 918 tFemtoV0->SetptotV0(::sqrt(tAODv0->Ptot2V0()));
ce7b3d98 919 //tFemtoV0->SetptPos(::sqrt(tAODv0->MomPosX()*tAODv0->MomPosX()+tAODv0->MomPosY()*tAODv0->MomPosY()));
920 //tFemtoV0->SetptotPos(::sqrt(tAODv0->Ptot2Pos()));
921 //tFemtoV0->SetptNeg(::sqrt(tAODv0->MomNegX()*tAODv0->MomNegX()+tAODv0->MomNegY()*tAODv0->MomNegY()));
922 //tFemtoV0->SetptotNeg(::sqrt(tAODv0->Ptot2Neg()));
973a91f8 923
924 tFemtoV0->SetidNeg(tAODv0->GetNegID());
925 //cout<<"tAODv0->GetNegID(): "<<tAODv0->GetNegID()<<endl;
926 //cout<<"tFemtoV0->IdNeg(): "<<tFemtoV0->IdNeg()<<endl;
927 tFemtoV0->SetidPos(tAODv0->GetPosID());
928
929 tFemtoV0->SetEtaV0(tAODv0->Eta());
930 tFemtoV0->SetPhiV0(tAODv0->Phi());
931 tFemtoV0->SetCosPointingAngle(tAODv0->CosPointingAngle(fV1));
932 //tFemtoV0->SetYV0(tAODv0->Y());
933
934 //void SetdedxNeg(float x);
935 //void SeterrdedxNeg(float x);//Gael 04Fev2002
936 //void SetlendedxNeg(float x);//Gael 04Fev2002
937 //void SetdedxPos(float x);
938 //void SeterrdedxPos(float x);//Gael 04Fev2002
939 //void SetlendedxPos(float x);//Gael 04Fev2002
940
ce7b3d98 941 //tFemtoV0->SetEtaPos(tAODv0->PseudoRapPos());
942 //tFemtoV0->SetEtaNeg(tAODv0->PseudoRapNeg());
973a91f8 943
944 AliAODTrack *trackpos = (AliAODTrack*)tAODv0->GetDaughter(0);
945 AliAODTrack *trackneg = (AliAODTrack*)tAODv0->GetDaughter(1);
946
947 if(trackpos && trackneg)
948 {
ce7b3d98 949 tFemtoV0->SetEtaPos(trackpos->Eta());
950 tFemtoV0->SetEtaNeg(trackneg->Eta());
951 tFemtoV0->SetptotPos(tAODv0->PProng(0));
952 tFemtoV0->SetptotNeg(tAODv0->PProng(1));
953 tFemtoV0->SetptPos(trackpos->Pt());
954 tFemtoV0->SetptNeg(trackneg->Pt());
955
973a91f8 956 //tFemtoV0->SetEtaPos(trackpos->Eta()); //tAODv0->PseudoRapPos()
957 //tFemtoV0->SetEtaNeg(trackneg->Eta()); //tAODv0->PseudoRapNeg()
958 tFemtoV0->SetTPCNclsPos(trackpos->GetTPCNcls());
959 tFemtoV0->SetTPCNclsNeg(trackneg->GetTPCNcls());
960 tFemtoV0->SetTPCclustersPos(trackpos->GetTPCClusterMap());
961 tFemtoV0->SetTPCclustersNeg(trackneg->GetTPCClusterMap());
962 tFemtoV0->SetTPCsharingPos(trackpos->GetTPCSharedMap());
963 tFemtoV0->SetTPCsharingNeg(trackneg->GetTPCSharedMap());
964 tFemtoV0->SetNdofPos(trackpos->Chi2perNDF());
965 tFemtoV0->SetNdofNeg(trackneg->Chi2perNDF());
966 tFemtoV0->SetStatusPos(trackpos->GetStatus());
967 tFemtoV0->SetStatusNeg(trackneg->GetStatus());
968
969 tFemtoV0->SetPosNSigmaTPCK(fAODpidUtil->NumberOfSigmasTPC(trackpos,AliPID::kKaon));
970 tFemtoV0->SetNegNSigmaTPCK(fAODpidUtil->NumberOfSigmasTPC(trackneg,AliPID::kKaon));
971 tFemtoV0->SetPosNSigmaTPCP(fAODpidUtil->NumberOfSigmasTPC(trackpos,AliPID::kProton));
972 tFemtoV0->SetNegNSigmaTPCP(fAODpidUtil->NumberOfSigmasTPC(trackneg,AliPID::kProton));
973 tFemtoV0->SetPosNSigmaTPCPi(fAODpidUtil->NumberOfSigmasTPC(trackpos,AliPID::kPion));
974 tFemtoV0->SetNegNSigmaTPCPi(fAODpidUtil->NumberOfSigmasTPC(trackneg,AliPID::kPion));
975
ce7b3d98 976
977 float bfield = 5*fMagFieldSign;
978 float globalPositionsAtRadiiPos[9][3];
979 GetGlobalPositionAtGlobalRadiiThroughTPC(trackpos,bfield,globalPositionsAtRadiiPos);
980 double tpcEntrancePos[3]={globalPositionsAtRadiiPos[0][0],globalPositionsAtRadiiPos[0][1],globalPositionsAtRadiiPos[0][2]};
981 double tpcExitPos[3]={globalPositionsAtRadiiPos[8][0],globalPositionsAtRadiiPos[8][1],globalPositionsAtRadiiPos[8][2]};
982
983 float globalPositionsAtRadiiNeg[9][3];
984 GetGlobalPositionAtGlobalRadiiThroughTPC(trackneg,bfield,globalPositionsAtRadiiNeg);
985 double tpcEntranceNeg[3]={globalPositionsAtRadiiNeg[0][0],globalPositionsAtRadiiNeg[0][1],globalPositionsAtRadiiNeg[0][2]};
986 double tpcExitNeg[3]={globalPositionsAtRadiiNeg[8][0],globalPositionsAtRadiiNeg[8][1],globalPositionsAtRadiiNeg[8][2]};
987
988 AliFemtoThreeVector tmpVec;
989 tmpVec.SetX(tpcEntrancePos[0]); tmpVec.SetX(tpcEntrancePos[1]); tmpVec.SetX(tpcEntrancePos[2]);
ba3c23a4 990 tFemtoV0->SetNominalTpcEntrancePointPos(tmpVec);
ce7b3d98 991
992 tmpVec.SetX(tpcExitPos[0]); tmpVec.SetX(tpcExitPos[1]); tmpVec.SetX(tpcExitPos[2]);
ba3c23a4 993 tFemtoV0->SetNominalTpcExitPointPos(tmpVec);
ce7b3d98 994
995 tmpVec.SetX(tpcEntranceNeg[0]); tmpVec.SetX(tpcEntranceNeg[1]); tmpVec.SetX(tpcEntranceNeg[2]);
996 tFemtoV0->SetNominalTpcEntrancePointNeg(tmpVec);
997
998 tmpVec.SetX(tpcExitNeg[0]); tmpVec.SetX(tpcExitNeg[1]); tmpVec.SetX(tpcExitNeg[2]);
ba3c23a4 999 tFemtoV0->SetNominalTpcExitPointNeg(tmpVec);
973a91f8 1000
ce7b3d98 1001 AliFemtoThreeVector vecTpcPos[9];
1002 AliFemtoThreeVector vecTpcNeg[9];
1003 for(int i=0;i<9;i++)
1004 {
1005 vecTpcPos[i].SetX(globalPositionsAtRadiiPos[i][0]); vecTpcPos[i].SetY(globalPositionsAtRadiiPos[i][1]); vecTpcPos[i].SetZ(globalPositionsAtRadiiPos[i][2]);
1006 vecTpcNeg[i].SetX(globalPositionsAtRadiiNeg[i][0]); vecTpcNeg[i].SetY(globalPositionsAtRadiiNeg[i][1]); vecTpcNeg[i].SetZ(globalPositionsAtRadiiNeg[i][2]);
1007 }
1008 tFemtoV0->SetNominalTpcPointPos(vecTpcPos);
1009 tFemtoV0->SetNominalTpcPointNeg(vecTpcNeg);
1010
1011 tFemtoV0->SetTPCMomentumPos(trackpos->GetTPCmomentum());
1012 tFemtoV0->SetTPCMomentumNeg(trackneg->GetTPCmomentum());
1013
1014 tFemtoV0->SetdedxPos(trackpos->GetTPCsignal());
1015 tFemtoV0->SetdedxNeg(trackneg->GetTPCsignal());
1016
63e066f7 1017 if((tFemtoV0->StatusPos()&AliESDtrack::kTOFpid)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTOFout)==0 )
973a91f8 1018 {
63e066f7 1019 if((tFemtoV0->StatusNeg()&AliESDtrack::kTOFpid)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTOFout)==0 )
973a91f8 1020 {
ce7b3d98 1021 tFemtoV0->SetPosNSigmaTOFK(-1000);
1022 tFemtoV0->SetNegNSigmaTOFK(-1000);
1023 tFemtoV0->SetPosNSigmaTOFP(-1000);
1024 tFemtoV0->SetNegNSigmaTOFP(-1000);
1025 tFemtoV0->SetPosNSigmaTOFPi(-1000);
1026 tFemtoV0->SetNegNSigmaTOFPi(-1000);
1027
1028 tFemtoV0->SetTOFProtonTimePos(-1000);
1029 tFemtoV0->SetTOFPionTimePos(-1000);
1030 tFemtoV0->SetTOFKaonTimePos(-1000);
1031 tFemtoV0->SetTOFProtonTimeNeg(-1000);
1032 tFemtoV0->SetTOFPionTimeNeg(-1000);
1033 tFemtoV0->SetTOFKaonTimeNeg(-1000);
973a91f8 1034 }
1035 }
1036 else
1037 {
1038 tFemtoV0->SetPosNSigmaTOFK(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kKaon));
1039 tFemtoV0->SetNegNSigmaTOFK(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kKaon));
1040 tFemtoV0->SetPosNSigmaTOFP(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kProton));
1041 tFemtoV0->SetNegNSigmaTOFP(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kProton));
1042 tFemtoV0->SetPosNSigmaTOFPi(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kPion));
1043 tFemtoV0->SetNegNSigmaTOFPi(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kPion));
ce7b3d98 1044
1045 double TOFSignalPos = trackpos->GetTOFsignal();
1046 double TOFSignalNeg = trackneg->GetTOFsignal();
1047 double pidPos[5];
1048 double pidNeg[5];
1049 trackpos->GetIntegratedTimes(pidPos);
1050 trackneg->GetIntegratedTimes(pidNeg);
1051
1052 tFemtoV0->SetTOFPionTimePos(TOFSignalPos-pidPos[2]);
1053 tFemtoV0->SetTOFKaonTimePos(TOFSignalPos-pidPos[3]);
1054 tFemtoV0->SetTOFProtonTimePos(TOFSignalPos-pidPos[4]);
1055 tFemtoV0->SetTOFPionTimeNeg(TOFSignalNeg-pidNeg[2]);
1056 tFemtoV0->SetTOFKaonTimeNeg(TOFSignalNeg-pidNeg[3]);
1057 tFemtoV0->SetTOFProtonTimeNeg(TOFSignalNeg-pidNeg[4]);
973a91f8 1058 }
1059 }
1060 else
1061 {
1062 tFemtoV0->SetStatusPos(999);
1063 tFemtoV0->SetStatusNeg(999);
1064 }
1065 tFemtoV0->SetOnFlyStatusV0(tAODv0->GetOnFlyStatus());
1066}
1067
76ce4b5b 1068void AliFemtoEventReaderAOD::SetFilterBit(UInt_t ibit)
1069{
1070 fFilterBit = (1 << (ibit));
1071}
1072
64536eaf 1073
1074void AliFemtoEventReaderAOD::SetFilterMask(int ibit)
1075{
1076 fFilterMask = ibit;
1077}
1078
76ce4b5b 1079void AliFemtoEventReaderAOD::SetReadMC(unsigned char a)
1080{
1081 fReadMC = a;
1082}
1083
973a91f8 1084
1085void AliFemtoEventReaderAOD::SetReadV0(unsigned char a)
1086{
1087 fReadV0 = a;
1088}
1089
1777446b 1090void AliFemtoEventReaderAOD::SetUseMultiplicity(EstEventMult aType)
1091{
1092 fEstEventMult = aType;
1093}
1094
76ce4b5b 1095AliAODMCParticle* AliFemtoEventReaderAOD::GetParticleWithLabel(TClonesArray *mcP, Int_t aLabel)
1096{
1097 if (aLabel < 0) return 0;
1098 AliAODMCParticle *aodP;
1099 Int_t posstack = 0;
1100 if (aLabel > mcP->GetEntries())
1101 posstack = mcP->GetEntries();
1102 else
1103 posstack = aLabel;
1104
1105 aodP = (AliAODMCParticle *) mcP->At(posstack);
1106 if (aodP->GetLabel() > posstack) {
1107 do {
1108 aodP = (AliAODMCParticle *) mcP->At(posstack);
1109 if (aodP->GetLabel() == aLabel) return aodP;
1110 posstack--;
1111 }
1112 while (posstack > 0);
1113 }
1114 else {
1115 do {
1116 aodP = (AliAODMCParticle *) mcP->At(posstack);
1117 if (aodP->GetLabel() == aLabel) return aodP;
1118 posstack++;
1119 }
1120 while (posstack < mcP->GetEntries());
1121 }
1122
1123 return 0;
1124}
1125
1126void AliFemtoEventReaderAOD::CopyPIDtoFemtoTrack(AliAODTrack *tAodTrack,
1127 AliFemtoTrack *tFemtoTrack)
1128{
41cb6c1e 1129
6691ea4f 1130 // copying DCA information (taking it from global tracks gives better resolution than from TPC-only)
41cb6c1e 1131
6691ea4f 1132 double impact[2];
1133 double covimpact[3];
1134
1135 if (!tAodTrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
1136 //cout << "sth went wrong with dca propagation" << endl;
1137 tFemtoTrack->SetImpactD(-1000.0);
1138 tFemtoTrack->SetImpactZ(-1000.0);
1139
1140 }
1141 else {
1142 tFemtoTrack->SetImpactD(impact[0]);
1143 tFemtoTrack->SetImpactZ(impact[1]);
1144 }
41cb6c1e 1145
76ce4b5b 1146 double aodpid[10];
1147 tAodTrack->GetPID(aodpid);
1148 tFemtoTrack->SetPidProbElectron(aodpid[0]);
1149 tFemtoTrack->SetPidProbMuon(aodpid[1]);
1150 tFemtoTrack->SetPidProbPion(aodpid[2]);
1151 tFemtoTrack->SetPidProbKaon(aodpid[3]);
1152 tFemtoTrack->SetPidProbProton(aodpid[4]);
1153
1154 aodpid[0] = -100000.0;
1155 aodpid[1] = -100000.0;
1156 aodpid[2] = -100000.0;
1157 aodpid[3] = -100000.0;
1158 aodpid[4] = -100000.0;
6691ea4f 1159
76ce4b5b 1160 double tTOF = 0.0;
1161
4eac0b05 1162 //what is that code? for what do we need it? nsigma values are not enough?
1163 if (tAodTrack->GetStatus() & AliESDtrack::kTOFpid) { //AliESDtrack::kTOFpid=0x8000
1164 tTOF = tAodTrack->GetTOFsignal();
1165 tAodTrack->GetIntegratedTimes(aodpid);
1777446b 1166
4eac0b05 1167 tTOF -= fAODpidUtil->GetTOFResponse().GetStartTime(tAodTrack->P());
1168 }
76ce4b5b 1169
4eac0b05 1170 tFemtoTrack->SetTofExpectedTimes(tTOF-aodpid[2], tTOF-aodpid[3], tTOF-aodpid[4]);
76ce4b5b 1171
1172 ////// TPC ////////////////////////////////////////////
1173
1174 float nsigmaTPCK=-1000.;
1175 float nsigmaTPCPi=-1000.;
1176 float nsigmaTPCP=-1000.;
1177
1178 // cout<<"in reader fESDpid"<<fESDpid<<endl;
1179
1180 if (tAodTrack->IsOn(AliESDtrack::kTPCpid)){ //AliESDtrack::kTPCpid=0x0080
1181 nsigmaTPCK = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kKaon);
1182 nsigmaTPCPi = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kPion);
1183 nsigmaTPCP = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kProton);
1184 }
1185
1186 tFemtoTrack->SetNSigmaTPCPi(nsigmaTPCPi);
1187 tFemtoTrack->SetNSigmaTPCK(nsigmaTPCK);
1188 tFemtoTrack->SetNSigmaTPCP(nsigmaTPCP);
1189
1190 tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF());
1191 tFemtoTrack->SetTPCncls(tAodTrack->GetTPCNcls());
1192 tFemtoTrack->SetTPCnclsF(tAodTrack->GetTPCNcls());
1193
1194 tFemtoTrack->SetTPCsignalN(1);
1195 tFemtoTrack->SetTPCsignalS(1);
1196 tFemtoTrack->SetTPCsignal(tAodTrack->GetTPCsignal());
1197
1198 ///////TOF//////////////////////
1199
1200 float vp=-1000.;
1201 float nsigmaTOFPi=-1000.;
1202 float nsigmaTOFK=-1000.;
1203 float nsigmaTOFP=-1000.;
1204
1205 if ((tAodTrack->GetStatus() & AliESDtrack::kTOFpid) && //AliESDtrack::kTOFpid=0x8000
1206 (tAodTrack->GetStatus() & AliESDtrack::kTOFout) && //AliESDtrack::kTOFout=0x2000
63e066f7 1207 (tAodTrack->GetStatus() & AliESDtrack::kTIME) ) //AliESDtrack::kTIME=0x80000000
76ce4b5b 1208 {
1209 if(tAodTrack->IsOn(AliESDtrack::kTOFpid)) //AliESDtrack::kTOFpid=0x8000
1210 {
1211
1212 nsigmaTOFPi = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kPion);
1213 nsigmaTOFK = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kKaon);
1214 nsigmaTOFP = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kProton);
1215
1216 Double_t len=200;// esdtrack->GetIntegratedLength(); !!!!!
1217 Double_t tof=tAodTrack->GetTOFsignal();
1218 if(tof > 0.) vp=len/tof/0.03;
1219 }
1220 }
1221 tFemtoTrack->SetVTOF(vp);
1222 tFemtoTrack->SetNSigmaTOFPi(nsigmaTOFPi);
1223 tFemtoTrack->SetNSigmaTOFK(nsigmaTOFK);
1224 tFemtoTrack->SetNSigmaTOFP(nsigmaTOFP);
1225
1226
1227 //////////////////////////////////////
1228
1229}
1230
1231void AliFemtoEventReaderAOD::SetCentralityPreSelection(double min, double max)
1232{
1233 fCentRange[0] = min; fCentRange[1] = max;
1234 fUsePreCent = 1;
1777446b 1235 fEstEventMult = kCentrality;
76ce4b5b 1236}
1237
973a91f8 1238void AliFemtoEventReaderAOD::SetNoCentrality(bool anocent)
1239{
1777446b 1240 if(anocent==false) {fEstEventMult=kCentrality;}
1241 else {fEstEventMult=kReference; fUsePreCent = 0; }
973a91f8 1242}
76ce4b5b 1243
1244void AliFemtoEventReaderAOD::SetAODpidUtil(AliAODpidUtil *aAODpidUtil)
1245{
1246 fAODpidUtil = aAODpidUtil;
1247 // printf("fAODpidUtil: %x\n",fAODpidUtil);
1248}
1249
1777446b 1250void AliFemtoEventReaderAOD::SetAODheader(AliAODHeader *aAODheader)
1251{
1252 fAODheader = aAODheader;
1253}
1254
76ce4b5b 1255
ba3c23a4 1256void AliFemtoEventReaderAOD::SetMagneticFieldSign(int s)
1257{
1258 if(s>0)
1259 fMagFieldSign = 1;
1260 else if(s<0)
1261 fMagFieldSign = -1;
1262 else
1263 fMagFieldSign = 0;
1264}
1265
5e2038a9 1266void AliFemtoEventReaderAOD::SetEPVZERO(Bool_t iepvz)
1267{
1268 fisEPVZ = iepvz;
1269}
76ce4b5b 1270
ce7b3d98 1271void AliFemtoEventReaderAOD::GetGlobalPositionAtGlobalRadiiThroughTPC(AliAODTrack *track, Float_t bfield, Float_t globalPositionsAtRadii[9][3])
1272{
1273 // Gets the global position of the track at nine different radii in the TPC
1274 // track is the track you want to propagate
1275 // bfield is the magnetic field of your event
1276 // globalPositionsAtRadii is the array of global positions in the radii and xyz
1277
1278 // Initialize the array to something indicating there was no propagation
1279 for(Int_t i=0;i<9;i++){
1280 for(Int_t j=0;j<3;j++){
1281 globalPositionsAtRadii[i][j]=-9999.;
1282 }
1283 }
1284
1285 // Make a copy of the track to not change parameters of the track
1286 AliExternalTrackParam etp; etp.CopyFromVTrack(track);
1287 //printf("\nAfter CopyFromVTrack\n");
1288 //etp.Print();
1289
1290 // The global position of the the track
1291 Double_t xyz[3]={-9999.,-9999.,-9999.};
1292
1293 // Counter for which radius we want
1294 Int_t iR=0;
1295 // The radii at which we get the global positions
1296 // IROC (OROC) from 84.1 cm to 132.1 cm (134.6 cm to 246.6 cm)
1297 Float_t Rwanted[9]={85.,105.,125.,145.,165.,185.,205.,225.,245.};
1298 // The global radius we are at
1299 Float_t globalRadius=0;
1300
1301 // Propagation is done in local x of the track
1302 for (Float_t x = etp.GetX();x<247.;x+=1.){ // GetX returns local coordinates
1303 // Starts at the tracks fX and goes outwards. x = 245 is the outer radial limit
1304 // of the TPC when the track is straight, i.e. has inifinite pt and doesn't get bent.
1305 // If the track's momentum is smaller than infinite, it will develop a y-component, which
1306 // adds to the global radius
1307
1308 // Stop if the propagation was not succesful. This can happen for low pt tracks
1309 // that don't reach outer radii
1310 if(!etp.PropagateTo(x,bfield))break;
1311 etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1312 globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
1313
1314 // Roughly reached the radius we want
1315 if(globalRadius > Rwanted[iR]){
1316
1317 // Bigger loop has bad precision, we're nearly one centimeter too far, go back in small steps.
1318 while (globalRadius>Rwanted[iR]){
1319 x-=.1;
1320 // printf("propagating to x %5.2f\n",x);
1321 if(!etp.PropagateTo(x,bfield))break;
1322 etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1323 globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
1324 }
1325 //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]);
1326 globalPositionsAtRadii[iR][0]=xyz[0];
1327 globalPositionsAtRadii[iR][1]=xyz[1];
1328 globalPositionsAtRadii[iR][2]=xyz[2];
1329 // Indicate we want the next radius
1330 iR+=1;
1331 }
1332 if(iR>=8){
1333 // TPC edge reached
1334 return;
1335 }
1336 }
1337}
76ce4b5b 1338
1339