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