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