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