]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/FEMTOSCOPY/AliFemto/AliFemtoEventReaderAOD.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
from carlos: Adding selection variables, QA, PTmerger and higher harmonics
[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;
906b97d0 568 if((fFilterBit == (1 << (7))) || fFilterMask == 128) {//for TPC Only tracks we have to copy PID information from corresponding global tracks
87ede832 569 aodtrackpid = fEvent->GetTrack(labels[-1-fEvent->GetTrack(i)->GetID()]);
906b97d0 570 }
571 else {
87ede832 572 aodtrackpid = fEvent->GetTrack(i);
906b97d0 573 }
574
87ede832 575 CopyPIDtoFemtoTrack(aodtrackpid, trackCopy);
576
577 if (mcP) {
578 // Fill the hidden information with the simulated data
579 // Int_t pLabel = aodtrack->GetLabel();
580 AliAODMCParticle *tPart = GetParticleWithLabel(mcP, (TMath::Abs(aodtrack->GetLabel())));
581
582 AliFemtoModelGlobalHiddenInfo *tInfo = new AliFemtoModelGlobalHiddenInfo();
583 double fpx=0.0, fpy=0.0, fpz=0.0, fpt=0.0;
584 if (!tPart) {
585 fpx = fV1[0];
586 fpy = fV1[1];
587 fpz = fV1[2];
588 tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
589 tInfo->SetPDGPid(0);
590 tInfo->SetTrueMomentum(0.0, 0.0, 0.0);
591 tInfo->SetEmissionPoint(0.0, 0.0, 0.0, 0.0);
592 tInfo->SetMass(0);
593 }
594 else {
595 // Check the mother information
596
597 // Using the new way of storing the freeze-out information
598 // Final state particle is stored twice on the stack
599 // one copy (mother) is stored with original freeze-out information
600 // and is not tracked
601 // the other one (daughter) is stored with primary vertex position
602 // and is tracked
603
604 // Freeze-out coordinates
605 fpx = tPart->Xv() - fV1[0];
606 fpy = tPart->Yv() - fV1[1];
607 fpz = tPart->Zv() - fV1[2];
608 // fpt = tPart->T();
609
610 tInfo->SetGlobalEmissionPoint(fpx, fpy, fpz);
611
50bb1be0 612
87ede832 613 fpx *= 1e13;
614 fpy *= 1e13;
615 fpz *= 1e13;
616 // fpt *= 1e13;
617
618 // cout << "Looking for mother ids " << endl;
619 if (motherids[TMath::Abs(aodtrack->GetLabel())]>0) {
620 // cout << "Got mother id" << endl;
621 AliAODMCParticle *mother = GetParticleWithLabel(mcP, motherids[TMath::Abs(aodtrack->GetLabel())]);
622 // Check if this is the same particle stored twice on the stack
623 if (mother) {
624 if ((mother->GetPdgCode() == tPart->GetPdgCode() || (mother->Px() == tPart->Px()))) {
625 // It is the same particle
626 // Read in the original freeze-out information
627 // and convert it from to [fm]
628
629 // EPOS style
630 // fpx = mother->Xv()*1e13*0.197327;
631 // fpy = mother->Yv()*1e13*0.197327;
632 // fpz = mother->Zv()*1e13*0.197327;
633 // fpt = mother->T() *1e13*0.197327*0.5;
634
635
636 // Therminator style
637 fpx = mother->Xv()*1e13;
638 fpy = mother->Yv()*1e13;
639 fpz = mother->Zv()*1e13;
640 // fpt = mother->T() *1e13*3e10;
641
642 }
643 }
644 }
645
646 // if (fRotateToEventPlane) {
647 // double tPhi = TMath::ATan2(fpy, fpx);
648 // double tRad = TMath::Hypot(fpx, fpy);
649
650 // fpx = tRad*TMath::Cos(tPhi - tReactionPlane);
651 // fpy = tRad*TMath::Sin(tPhi - tReactionPlane);
652 // }
653
654 tInfo->SetPDGPid(tPart->GetPdgCode());
655
656 // if (fRotateToEventPlane) {
657 // double tPhi = TMath::ATan2(tPart->Py(), tPart->Px());
658 // double tRad = TMath::Hypot(tPart->Px(), tPart->Py());
659
660 // tInfo->SetTrueMomentum(tRad*TMath::Cos(tPhi - tReactionPlane),
661 // tRad*TMath::Sin(tPhi - tReactionPlane),
662 // tPart->Pz());
663 // }
664 // else
665 tInfo->SetTrueMomentum(tPart->Px(), tPart->Py(), tPart->Pz());
666 Double_t mass2 = (tPart->E() *tPart->E() -
667 tPart->Px()*tPart->Px() -
668 tPart->Py()*tPart->Py() -
669 tPart->Pz()*tPart->Pz());
670 if (mass2>0.0)
671 tInfo->SetMass(TMath::Sqrt(mass2));
672 else
673 tInfo->SetMass(0.0);
674
675 tInfo->SetEmissionPoint(fpx, fpy, fpz, fpt);
50bb1be0 676
ab558586 677 // // fillDCA
678 // //if (TMath::Abs(impact[0]) > 0.001) {
679 // if (tPart->IsPhysicalPrimary()){
680 // tInfo->SetPartOrigin(0);
681 // // trackCopy->SetImpactDprim(impact[0]);
290ffcb9 682 // //cout << "Read prim" << endl;
ab558586 683 // }
684 // else if (tPart->IsSecondaryFromWeakDecay()) {
685 // tInfo->SetPartOrigin(1);
686 // // trackCopy->SetImpactDweak(impact[0]);
687 // //cout << "Read wea" << endl;
688 // }
689 // else if (tPart->IsSecondaryFromMaterial()) {
690 // tInfo->SetPartOrigin(2);
691 // // trackCopy->SetImpactDmat(impact[0]);
692 // //cout << "Read mat" << endl;
693 // }
694 // //}
695 // // end fillDCA
50bb1be0 696
697
87ede832 698 }
699 trackCopy->SetHiddenInfo(tInfo);
76ce4b5b 700 }
87ede832 701
702 double pxyz[3];
703
704 //AliExternalTrackParam *param = new AliExternalTrackParam(*aodtrack->GetInnerParam());
705 trackCopy->SetInnerMomentum(aodtrack->GetTPCmomentum());
706
707 aodtrack->PxPyPz(pxyz);//reading noconstarined momentum
708 const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
709 // Check the sanity of the tracks - reject zero momentum tracks
710 if (ktP.Mag() == 0) {
711 delete trackCopy;
712 continue;
713 }
714 // }
715
316b08d6 716 tEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
87ede832 717 realnofTracks++;//real number of tracks
718 }
719
720 tEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
76ce4b5b 721 tEvent->SetNormalizedMult(tracksPrim);
722
2bcd4b96 723 if (cent) {
724 tEvent->SetCentralityV0(cent->GetCentralityPercentile("V0M"));
18757d69 725 tEvent->SetCentralityV0A(cent->GetCentralityPercentile("V0A"));
726 tEvent->SetCentralityV0C(cent->GetCentralityPercentile("V0C"));
2bcd4b96 727 tEvent->SetCentralityZNA(cent->GetCentralityPercentile("ZNA"));
18757d69 728 tEvent->SetCentralityZNC(cent->GetCentralityPercentile("ZNC"));
2bcd4b96 729 tEvent->SetCentralityCL1(cent->GetCentralityPercentile("CL1"));
18757d69 730 tEvent->SetCentralityCL0(cent->GetCentralityPercentile("CL0"));
731 tEvent->SetCentralityTKL(cent->GetCentralityPercentile("TKL"));
732 tEvent->SetCentralityFMD(cent->GetCentralityPercentile("FMD"));
733 tEvent->SetCentralityFMD(cent->GetCentralityPercentile("NPA"));
734 // tEvent->SetCentralitySPD1(cent->GetCentralityPercentile("CL1"));
607aa748 735 tEvent->SetCentralityTrk(cent->GetCentralityPercentile("TRK"));
736 tEvent->SetCentralityCND(cent->GetCentralityPercentile("CND"));
2bcd4b96 737 }
738
1777446b 739 if (fEstEventMult==kCentrality) {
2bcd4b96 740 //AliCentrality *cent = fEvent->GetCentrality();
4eac0b05 741 //cout<<"AliFemtoEventReaderAOD:"<<lrint(10*cent->GetCentralityPercentile("V0M"))<<endl;
973a91f8 742 if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("V0M")));
743 // if (cent) tEvent->SetNormalizedMult((int) cent->GetCentralityPercentile("V0M"));
c863b24a 744 }
18757d69 745 else if (fEstEventMult==kCentralityV0A) {
746 if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("V0A")));
747 }
748 else if (fEstEventMult==kCentralityV0C) {
749 if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("V0C")));
750 }
c863b24a 751 else if (fEstEventMult==kCentralityZNA) {
c863b24a 752 if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("ZNA")));
c863b24a 753 }
87ede832 754 else if (fEstEventMult==kCentralityZNC) {
18757d69 755 if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("ZNC")));
756 }
c863b24a 757 else if (fEstEventMult==kCentralityCL1) {
c863b24a 758 if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("CL1")));
76ce4b5b 759 }
18757d69 760 else if (fEstEventMult==kCentralityCL0) {
761 if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("CL0")));
762 }
607aa748 763 else if (fEstEventMult==kCentralityTRK) {
764 if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("TRK")));
765 }
18757d69 766 else if (fEstEventMult==kCentralityTKL) {
767 if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("TKL")));
768 }
607aa748 769 else if (fEstEventMult==kCentralityCND) {
770 if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("CND")));
771 }
18757d69 772 else if (fEstEventMult==kCentralityNPA) {
773 if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("NPA")));
774 }
87ede832 775 else if (fEstEventMult==kCentralityFMD) {
18757d69 776 if (cent) tEvent->SetNormalizedMult(lrint(10*cent->GetCentralityPercentile("FMD")));
777 }
1777446b 778 else if(fEstEventMult==kGlobalCount){
ba3c23a4 779 tEvent->SetNormalizedMult(tNormMult); //particles counted in the loop, trying to reproduce GetReferenceMultiplicity. If better (default) method appears it should be changed
780 }
1777446b 781 else if(fEstEventMult==kReference)
87ede832 782 {
783 tEvent->SetNormalizedMult(fAODheader->GetRefMultiplicity());
784 }
1777446b 785 else if(fEstEventMult==kTPCOnlyRef)
87ede832 786 {
787 tEvent->SetNormalizedMult(fAODheader->GetTPConlyRefMultiplicity());
788 }
1938eadf 789 else if(fEstEventMult == kVZERO)
87ede832 790 {
791 Float_t multV0 = 0;
792 for (Int_t i=0; i<64; i++)
793 multV0 += fEvent->GetVZEROData()->GetMultiplicity(i);
794 tEvent->SetNormalizedMult(multV0);
795 }
76ce4b5b 796
797 if (mcP) delete [] motherids;
798
87ede832 799 // cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
973a91f8 800
801 if(fReadV0)
87ede832 802 {
803 int count_pass = 0;
804 for (Int_t i = 0; i < fEvent->GetNumberOfV0s(); i++) {
805 AliAODv0* aodv0 = fEvent->GetV0(i);
806 if (!aodv0) continue;
807 if(aodv0->GetNDaughters()>2) continue;
808 if(aodv0->GetNProngs()>2) continue;
809 if(aodv0->GetCharge()!=0) continue;
810 if(aodv0->ChargeProng(0)==aodv0->ChargeProng(1)) continue;
811 if(aodv0->CosPointingAngle(fV1)<0.998) continue;
812 AliFemtoV0* trackCopyV0 = new AliFemtoV0();
813 count_pass++;
814 CopyAODtoFemtoV0(aodv0, trackCopyV0);
815 tEvent->V0Collection()->push_back(trackCopyV0);
816 //cout<<"Pushback v0 to v0collection"<<endl;
973a91f8 817 }
87ede832 818 }
973a91f8 819
76ce4b5b 820}
821
87ede832 822void AliFemtoEventReaderAOD::CopyAODtoFemtoTrack(AliAODTrack *tAodTrack,
823 AliFemtoTrack *tFemtoTrack
824 // AliPWG2AODTrack *tPWG2AODTrack
825 )
76ce4b5b 826{
827 // Copy the track information from the AOD into the internal AliFemtoTrack
828 // If it exists, use the additional information from the PWG2 AOD
829
830 // Primary Vertex position
87ede832 831
76ce4b5b 832 fEvent->GetPrimaryVertex()->GetPosition(fV1);
833 // fEvent->GetPrimaryVertex()->GetXYZ(fV1);
834
835 tFemtoTrack->SetCharge(tAodTrack->Charge());
87ede832 836
76ce4b5b 837 double pxyz[3];
838 tAodTrack->PxPyPz(pxyz);//reading noconstrained momentum
290ffcb9 839
76ce4b5b 840 AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
841 tFemtoTrack->SetP(v);//setting momentum
842 tFemtoTrack->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
843 const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
87ede832 844 //setting track helix
76ce4b5b 845 const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
87ede832 846 AliFmPhysicalHelixD helix(ktP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(tFemtoTrack->Charge()));
76ce4b5b 847 tFemtoTrack->SetHelix(helix);
87ede832 848
76ce4b5b 849 // Flags
850 tFemtoTrack->SetTrackId(tAodTrack->GetID());
851 tFemtoTrack->SetFlags(tAodTrack->GetFlags());
852 tFemtoTrack->SetLabel(tAodTrack->GetLabel());
87ede832 853
854 // Track quality information
76ce4b5b 855 float covmat[6];
87ede832 856 tAodTrack->GetCovMatrix(covmat);
857
50bb1be0 858 if (!fDCAglobalTrack) {
76ce4b5b 859
906b97d0 860 tFemtoTrack->SetImpactD(tAodTrack->DCA());
861 tFemtoTrack->SetImpactZ(tAodTrack->ZAtDCA());
862
76ce4b5b 863
ef838a91 864 // double impact[2];
865 // double covimpact[3];
866
867 // AliAODTrack* trk_clone = (AliAODTrack*)tAodTrack->Clone("trk_clone");
868 // // if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
869 // // if (!tAodTrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
870 // if (!trk_clone->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
871 // //cout << "sth went wrong with dca propagation" << endl;
872 // tFemtoTrack->SetImpactD(-1000.0);
873 // tFemtoTrack->SetImpactZ(-1000.0);
874
875 // }
876 // else {
877 // tFemtoTrack->SetImpactD(impact[0]);
878 // tFemtoTrack->SetImpactZ(impact[1]);
879 // }
880 // delete trk_clone;
e42e0efc 881
50bb1be0 882 }
76ce4b5b 883
884 // if (TMath::Abs(tAodTrack->Xv()) > 0.00000000001)
885 // tFemtoTrack->SetImpactD(TMath::Hypot(tAodTrack->Xv(), tAodTrack->Yv())*(tAodTrack->Xv()/TMath::Abs(tAodTrack->Xv())));
886 // else
887 // tFemtoTrack->SetImpactD(0.0);
888 // tFemtoTrack->SetImpactD(tAodTrack->DCA());
87ede832 889
76ce4b5b 890 // tFemtoTrack->SetImpactZ(tAodTrack->ZAtDCA());
891
892
893 // tFemtoTrack->SetImpactD(TMath::Hypot(tAodTrack->Xv() - fV1[0], tAodTrack->Yv() - fV1[1]));
894 // tFemtoTrack->SetImpactZ(tAodTrack->Zv() - fV1[2]);
895
896
87ede832 897 // cout
898 // << "dca" << TMath::Hypot(tAodTrack->Xv() - fV1[0], tAodTrack->Yv() - fV1[1])
899 // << "xv - fv10 = "<< tAodTrack->Xv() - fV1[0]
900 // << tAodTrack->Yv() - fV1[1]
901// << "xv = " << tAodTrack->Xv() << endl
902// << "fv1[0] = " << fV1[0] << endl
903// << "yv = " << tAodTrack->Yv() << endl
904// << "fv1[1] = " << fV1[1] << endl
905// << "zv = " << tAodTrack->Zv() << endl
906// << "fv1[2] = " << fV1[2] << endl
907// << "impact[0] = " << impact[0] << endl
908// << "impact[1] = " << impact[1] << endl
76ce4b5b 909// << endl << endl ;
910
911 tFemtoTrack->SetCdd(covmat[0]);
912 tFemtoTrack->SetCdz(covmat[1]);
913 tFemtoTrack->SetCzz(covmat[2]);
973a91f8 914 tFemtoTrack->SetITSchi2(tAodTrack->Chi2perNDF());
915 tFemtoTrack->SetITSncls(tAodTrack->GetITSNcls());
916 tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF());
917 tFemtoTrack->SetTPCncls(tAodTrack->GetTPCNcls());
918 tFemtoTrack->SetTPCnclsF(tAodTrack->GetTPCNcls());
87ede832 919 tFemtoTrack->SetTPCsignalN(1);
920 tFemtoTrack->SetTPCsignalS(1);
76ce4b5b 921 tFemtoTrack->SetTPCsignal(tAodTrack->GetTPCsignal());
922
923// if (tPWG2AODTrack) {
924// // Copy the PWG2 specific information if it exists
925// tFemtoTrack->SetTPCClusterMap(tPWG2AODTrack->GetTPCClusterMap());
926// tFemtoTrack->SetTPCSharedMap(tPWG2AODTrack->GetTPCSharedMap());
87ede832 927
76ce4b5b 928// double xtpc[3] = {0,0,0};
929// tPWG2AODTrack->GetTPCNominalEntrancePoint(xtpc);
930// tFemtoTrack->SetNominalTPCEntrancePoint(xtpc);
931// tPWG2AODTrack->GetTPCNominalExitPoint(xtpc);
932// tFemtoTrack->SetNominalTPCExitPoint(xtpc);
933// }
934// else {
87ede832 935 // If not use dummy values
76ce4b5b 936 tFemtoTrack->SetTPCClusterMap(tAodTrack->GetTPCClusterMap());
937 tFemtoTrack->SetTPCSharedMap(tAodTrack->GetTPCSharedMap());
87ede832 938
ce7b3d98 939 float globalPositionsAtRadii[9][3];
940 float bfield = 5*fMagFieldSign;
c710fd31 941
ce7b3d98 942 GetGlobalPositionAtGlobalRadiiThroughTPC(tAodTrack,bfield,globalPositionsAtRadii);
943 double tpcEntrance[3]={globalPositionsAtRadii[0][0],globalPositionsAtRadii[0][1],globalPositionsAtRadii[0][2]};
944 double **tpcPositions;
945 tpcPositions = new double*[9];
946 for(int i=0;i<9;i++)
947 tpcPositions[i] = new double[3];
948 double tpcExit[3]={globalPositionsAtRadii[8][0],globalPositionsAtRadii[8][1],globalPositionsAtRadii[8][2]};
949 for(int i=0;i<9;i++)
87ede832 950 {
951 tpcPositions[i][0] = globalPositionsAtRadii[i][0];
952 tpcPositions[i][1] = globalPositionsAtRadii[i][1];
953 tpcPositions[i][2] = globalPositionsAtRadii[i][2];
954 }
7db4e4cf 955
ce7b3d98 956 tFemtoTrack->SetNominalTPCEntrancePoint(tpcEntrance);
957 tFemtoTrack->SetNominalTPCPoints(tpcPositions);
958 tFemtoTrack->SetNominalTPCExitPoint(tpcExit);
87ede832 959 for(int i=0;i<9;i++)
960 delete [] tpcPositions[i];
961 delete [] tpcPositions;
76ce4b5b 962 // }
87ede832 963
76ce4b5b 964 // // cout << "Track has " << TMath::Hypot(tAodTrack->Xv(), tAodTrack->Yv()) << " " << tAodTrack->Zv() << " " << tAodTrack->GetTPCNcls() << endl;
87ede832 965
966
76ce4b5b 967 int indexes[3];
968 for (int ik=0; ik<3; ik++) {
969 indexes[ik] = 0;
970 }
971 tFemtoTrack->SetKinkIndexes(indexes);
4eac0b05 972
973
974 for (int ii=0; ii<6; ii++){
975 tFemtoTrack->SetITSHitOnLayer(ii,tAodTrack->HasPointOnITSLayer(ii));
976 }
977
978
76ce4b5b 979}
980
973a91f8 981void AliFemtoEventReaderAOD::CopyAODtoFemtoV0(AliAODv0 *tAODv0, AliFemtoV0 *tFemtoV0)
982{
ce7b3d98 983 tFemtoV0->SetdecayLengthV0(tAODv0->DecayLength(fV1));
973a91f8 984 tFemtoV0->SetdecayVertexV0X(tAODv0->DecayVertexV0X());
985 tFemtoV0->SetdecayVertexV0Y(tAODv0->DecayVertexV0Y());
986 tFemtoV0->SetdecayVertexV0Z(tAODv0->DecayVertexV0Z());
987 AliFemtoThreeVector decayvertex(tAODv0->DecayVertexV0X(),tAODv0->DecayVertexV0Y(),tAODv0->DecayVertexV0Z());
988 tFemtoV0->SetdecayVertexV0(decayvertex);
989 tFemtoV0->SetdcaV0Daughters(tAODv0->DcaV0Daughters());
990 tFemtoV0->SetdcaV0ToPrimVertex(tAODv0->DcaV0ToPrimVertex());
991 tFemtoV0->SetdcaPosToPrimVertex(tAODv0->DcaPosToPrimVertex());
992 tFemtoV0->SetdcaNegToPrimVertex(tAODv0->DcaNegToPrimVertex());
993 tFemtoV0->SetmomPosX(tAODv0->MomPosX());
994 tFemtoV0->SetmomPosY(tAODv0->MomPosY());
995 tFemtoV0->SetmomPosZ(tAODv0->MomPosZ());
996 AliFemtoThreeVector mompos(tAODv0->MomPosX(),tAODv0->MomPosY(),tAODv0->MomPosZ());
997 tFemtoV0->SetmomPos(mompos);
998 tFemtoV0->SetmomNegX(tAODv0->MomNegX());
999 tFemtoV0->SetmomNegY(tAODv0->MomNegY());
1000 tFemtoV0->SetmomNegZ(tAODv0->MomNegZ());
1001 AliFemtoThreeVector momneg(tAODv0->MomNegX(),tAODv0->MomNegY(),tAODv0->MomNegZ());
1002 tFemtoV0->SetmomNeg(momneg);
1003
1004 //jest cos takiego w AliFemtoV0.h czego nie ma w AliAODv0.h
87ede832 1005 //void SettpcHitsPos(const int& i);
1006 //void SettpcHitsNeg(const int& i);
973a91f8 1007
1008 //void SetTrackTopologyMapPos(unsigned int word, const unsigned long& m);
1009 //void SetTrackTopologyMapNeg(unsigned int word, const unsigned long& m);
1010
1011 tFemtoV0->SetmomV0X(tAODv0->MomV0X());
1012 tFemtoV0->SetmomV0Y(tAODv0->MomV0Y());
1013 tFemtoV0->SetmomV0Z(tAODv0->MomV0Z());
1014 AliFemtoThreeVector momv0(tAODv0->MomV0X(),tAODv0->MomV0Y(),tAODv0->MomV0Z());
1015 tFemtoV0->SetmomV0(momv0);
1016 tFemtoV0->SetalphaV0(tAODv0->AlphaV0());
1017 tFemtoV0->SetptArmV0(tAODv0->PtArmV0());
1018 tFemtoV0->SeteLambda(tAODv0->ELambda());
1019 tFemtoV0->SeteK0Short(tAODv0->EK0Short());
1020 tFemtoV0->SetePosProton(tAODv0->EPosProton());
1021 tFemtoV0->SeteNegProton(tAODv0->ENegProton());
1022 tFemtoV0->SetmassLambda(tAODv0->MassLambda());
1023 tFemtoV0->SetmassAntiLambda(tAODv0->MassAntiLambda());
1024 tFemtoV0->SetmassK0Short(tAODv0->MassK0Short());
1025 tFemtoV0->SetrapLambda(tAODv0->RapLambda());
1026 tFemtoV0->SetrapK0Short(tAODv0->RapK0Short());
87ede832 1027
1028 //void SetcTauLambda( float x);
1029 //void SetcTauK0Short( float x);
1030
ce7b3d98 1031 //tFemtoV0->SetptV0(::sqrt(tAODv0->Pt2V0())); //!
1032 tFemtoV0->SetptV0(tAODv0->Pt());
973a91f8 1033 tFemtoV0->SetptotV0(::sqrt(tAODv0->Ptot2V0()));
ce7b3d98 1034 //tFemtoV0->SetptPos(::sqrt(tAODv0->MomPosX()*tAODv0->MomPosX()+tAODv0->MomPosY()*tAODv0->MomPosY()));
1035 //tFemtoV0->SetptotPos(::sqrt(tAODv0->Ptot2Pos()));
1036 //tFemtoV0->SetptNeg(::sqrt(tAODv0->MomNegX()*tAODv0->MomNegX()+tAODv0->MomNegY()*tAODv0->MomNegY()));
1037 //tFemtoV0->SetptotNeg(::sqrt(tAODv0->Ptot2Neg()));
87ede832 1038
973a91f8 1039 tFemtoV0->SetidNeg(tAODv0->GetNegID());
1040 //cout<<"tAODv0->GetNegID(): "<<tAODv0->GetNegID()<<endl;
1041 //cout<<"tFemtoV0->IdNeg(): "<<tFemtoV0->IdNeg()<<endl;
1042 tFemtoV0->SetidPos(tAODv0->GetPosID());
1043
1044 tFemtoV0->SetEtaV0(tAODv0->Eta());
1045 tFemtoV0->SetPhiV0(tAODv0->Phi());
1046 tFemtoV0->SetCosPointingAngle(tAODv0->CosPointingAngle(fV1));
1047 //tFemtoV0->SetYV0(tAODv0->Y());
1048
1049 //void SetdedxNeg(float x);
1050 //void SeterrdedxNeg(float x);//Gael 04Fev2002
1051 //void SetlendedxNeg(float x);//Gael 04Fev2002
1052 //void SetdedxPos(float x);
1053 //void SeterrdedxPos(float x);//Gael 04Fev2002
1054 //void SetlendedxPos(float x);//Gael 04Fev2002
1055
ce7b3d98 1056 //tFemtoV0->SetEtaPos(tAODv0->PseudoRapPos());
1057 //tFemtoV0->SetEtaNeg(tAODv0->PseudoRapNeg());
973a91f8 1058
1059 AliAODTrack *trackpos = (AliAODTrack*)tAODv0->GetDaughter(0);
1060 AliAODTrack *trackneg = (AliAODTrack*)tAODv0->GetDaughter(1);
1061
1062 if(trackpos && trackneg)
87ede832 1063 {
1064 tFemtoV0->SetEtaPos(trackpos->Eta());
1065 tFemtoV0->SetEtaNeg(trackneg->Eta());
1066 tFemtoV0->SetptotPos(tAODv0->PProng(0));
1067 tFemtoV0->SetptotNeg(tAODv0->PProng(1));
1068 tFemtoV0->SetptPos(trackpos->Pt());
1069 tFemtoV0->SetptNeg(trackneg->Pt());
1070
1071 //tFemtoV0->SetEtaPos(trackpos->Eta()); //tAODv0->PseudoRapPos()
1072 //tFemtoV0->SetEtaNeg(trackneg->Eta()); //tAODv0->PseudoRapNeg()
1073 tFemtoV0->SetTPCNclsPos(trackpos->GetTPCNcls());
1074 tFemtoV0->SetTPCNclsNeg(trackneg->GetTPCNcls());
1075 tFemtoV0->SetTPCclustersPos(trackpos->GetTPCClusterMap());
1076 tFemtoV0->SetTPCclustersNeg(trackneg->GetTPCClusterMap());
1077 tFemtoV0->SetTPCsharingPos(trackpos->GetTPCSharedMap());
1078 tFemtoV0->SetTPCsharingNeg(trackneg->GetTPCSharedMap());
1079 tFemtoV0->SetNdofPos(trackpos->Chi2perNDF());
1080 tFemtoV0->SetNdofNeg(trackneg->Chi2perNDF());
1081 tFemtoV0->SetStatusPos(trackpos->GetStatus());
1082 tFemtoV0->SetStatusNeg(trackneg->GetStatus());
1083
1084 tFemtoV0->SetPosNSigmaTPCK(fAODpidUtil->NumberOfSigmasTPC(trackpos,AliPID::kKaon));
1085 tFemtoV0->SetNegNSigmaTPCK(fAODpidUtil->NumberOfSigmasTPC(trackneg,AliPID::kKaon));
1086 tFemtoV0->SetPosNSigmaTPCP(fAODpidUtil->NumberOfSigmasTPC(trackpos,AliPID::kProton));
1087 tFemtoV0->SetNegNSigmaTPCP(fAODpidUtil->NumberOfSigmasTPC(trackneg,AliPID::kProton));
1088 tFemtoV0->SetPosNSigmaTPCPi(fAODpidUtil->NumberOfSigmasTPC(trackpos,AliPID::kPion));
1089 tFemtoV0->SetNegNSigmaTPCPi(fAODpidUtil->NumberOfSigmasTPC(trackneg,AliPID::kPion));
1090
1091
1092 float bfield = 5*fMagFieldSign;
1093 float globalPositionsAtRadiiPos[9][3];
1094 GetGlobalPositionAtGlobalRadiiThroughTPC(trackpos,bfield,globalPositionsAtRadiiPos);
1095 double tpcEntrancePos[3]={globalPositionsAtRadiiPos[0][0],globalPositionsAtRadiiPos[0][1],globalPositionsAtRadiiPos[0][2]};
1096 double tpcExitPos[3]={globalPositionsAtRadiiPos[8][0],globalPositionsAtRadiiPos[8][1],globalPositionsAtRadiiPos[8][2]};
1097
1098 float globalPositionsAtRadiiNeg[9][3];
1099 GetGlobalPositionAtGlobalRadiiThroughTPC(trackneg,bfield,globalPositionsAtRadiiNeg);
1100 double tpcEntranceNeg[3]={globalPositionsAtRadiiNeg[0][0],globalPositionsAtRadiiNeg[0][1],globalPositionsAtRadiiNeg[0][2]};
1101 double tpcExitNeg[3]={globalPositionsAtRadiiNeg[8][0],globalPositionsAtRadiiNeg[8][1],globalPositionsAtRadiiNeg[8][2]};
1102
1103 AliFemtoThreeVector tmpVec;
290ffcb9 1104 tmpVec.SetX(tpcEntrancePos[0]); tmpVec.SetY(tpcEntrancePos[1]); tmpVec.SetZ(tpcEntrancePos[2]);
87ede832 1105 tFemtoV0->SetNominalTpcEntrancePointPos(tmpVec);
1106
290ffcb9 1107 tmpVec.SetX(tpcExitPos[0]); tmpVec.SetY(tpcExitPos[1]); tmpVec.SetZ(tpcExitPos[2]);
87ede832 1108 tFemtoV0->SetNominalTpcExitPointPos(tmpVec);
1109
290ffcb9 1110 tmpVec.SetX(tpcEntranceNeg[0]); tmpVec.SetY(tpcEntranceNeg[1]); tmpVec.SetZ(tpcEntranceNeg[2]);
87ede832 1111 tFemtoV0->SetNominalTpcEntrancePointNeg(tmpVec);
1112
290ffcb9 1113 tmpVec.SetX(tpcExitNeg[0]); tmpVec.SetY(tpcExitNeg[1]); tmpVec.SetZ(tpcExitNeg[2]);
87ede832 1114 tFemtoV0->SetNominalTpcExitPointNeg(tmpVec);
1115
1116 AliFemtoThreeVector vecTpcPos[9];
1117 AliFemtoThreeVector vecTpcNeg[9];
1118 for(int i=0;i<9;i++)
973a91f8 1119 {
87ede832 1120 vecTpcPos[i].SetX(globalPositionsAtRadiiPos[i][0]); vecTpcPos[i].SetY(globalPositionsAtRadiiPos[i][1]); vecTpcPos[i].SetZ(globalPositionsAtRadiiPos[i][2]);
1121 vecTpcNeg[i].SetX(globalPositionsAtRadiiNeg[i][0]); vecTpcNeg[i].SetY(globalPositionsAtRadiiNeg[i][1]); vecTpcNeg[i].SetZ(globalPositionsAtRadiiNeg[i][2]);
1122 }
1123 tFemtoV0->SetNominalTpcPointPos(vecTpcPos);
1124 tFemtoV0->SetNominalTpcPointNeg(vecTpcNeg);
ce7b3d98 1125
87ede832 1126 tFemtoV0->SetTPCMomentumPos(trackpos->GetTPCmomentum());
1127 tFemtoV0->SetTPCMomentumNeg(trackneg->GetTPCmomentum());
ce7b3d98 1128
87ede832 1129 tFemtoV0->SetdedxPos(trackpos->GetTPCsignal());
1130 tFemtoV0->SetdedxNeg(trackneg->GetTPCsignal());
ce7b3d98 1131
87ede832 1132 if((tFemtoV0->StatusPos()&AliESDtrack::kTOFpid)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTOFout)==0 )
1133 {
1134 if((tFemtoV0->StatusNeg()&AliESDtrack::kTOFpid)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTOFout)==0 )
973a91f8 1135 {
ce7b3d98 1136 tFemtoV0->SetPosNSigmaTOFK(-1000);
1137 tFemtoV0->SetNegNSigmaTOFK(-1000);
1138 tFemtoV0->SetPosNSigmaTOFP(-1000);
1139 tFemtoV0->SetNegNSigmaTOFP(-1000);
1140 tFemtoV0->SetPosNSigmaTOFPi(-1000);
1141 tFemtoV0->SetNegNSigmaTOFPi(-1000);
1142
1143 tFemtoV0->SetTOFProtonTimePos(-1000);
1144 tFemtoV0->SetTOFPionTimePos(-1000);
1145 tFemtoV0->SetTOFKaonTimePos(-1000);
1146 tFemtoV0->SetTOFProtonTimeNeg(-1000);
1147 tFemtoV0->SetTOFPionTimeNeg(-1000);
1148 tFemtoV0->SetTOFKaonTimeNeg(-1000);
973a91f8 1149 }
973a91f8 1150 }
87ede832 1151 else
973a91f8 1152 {
87ede832 1153 tFemtoV0->SetPosNSigmaTOFK(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kKaon));
1154 tFemtoV0->SetNegNSigmaTOFK(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kKaon));
1155 tFemtoV0->SetPosNSigmaTOFP(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kProton));
1156 tFemtoV0->SetNegNSigmaTOFP(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kProton));
1157 tFemtoV0->SetPosNSigmaTOFPi(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kPion));
1158 tFemtoV0->SetNegNSigmaTOFPi(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kPion));
1159
1160 double TOFSignalPos = trackpos->GetTOFsignal();
1161 double TOFSignalNeg = trackneg->GetTOFsignal();
1162 double pidPos[5];
1163 double pidNeg[5];
1164 trackpos->GetIntegratedTimes(pidPos);
1165 trackneg->GetIntegratedTimes(pidNeg);
1166
1167 tFemtoV0->SetTOFPionTimePos(TOFSignalPos-pidPos[2]);
1168 tFemtoV0->SetTOFKaonTimePos(TOFSignalPos-pidPos[3]);
1169 tFemtoV0->SetTOFProtonTimePos(TOFSignalPos-pidPos[4]);
1170 tFemtoV0->SetTOFPionTimeNeg(TOFSignalNeg-pidNeg[2]);
1171 tFemtoV0->SetTOFKaonTimeNeg(TOFSignalNeg-pidNeg[3]);
1172 tFemtoV0->SetTOFProtonTimeNeg(TOFSignalNeg-pidNeg[4]);
973a91f8 1173 }
87ede832 1174 }
1175 else
1176 {
1177 tFemtoV0->SetStatusPos(999);
1178 tFemtoV0->SetStatusNeg(999);
1179 }
973a91f8 1180 tFemtoV0->SetOnFlyStatusV0(tAODv0->GetOnFlyStatus());
1181}
1182
76ce4b5b 1183void AliFemtoEventReaderAOD::SetFilterBit(UInt_t ibit)
1184{
1185 fFilterBit = (1 << (ibit));
1186}
1187
64536eaf 1188
1189void AliFemtoEventReaderAOD::SetFilterMask(int ibit)
1190{
1191 fFilterMask = ibit;
1192}
1193
76ce4b5b 1194void AliFemtoEventReaderAOD::SetReadMC(unsigned char a)
1195{
1196 fReadMC = a;
1197}
1198
973a91f8 1199
1200void AliFemtoEventReaderAOD::SetReadV0(unsigned char a)
1201{
1202 fReadV0 = a;
1203}
1204
1777446b 1205void AliFemtoEventReaderAOD::SetUseMultiplicity(EstEventMult aType)
1206{
1207 fEstEventMult = aType;
1208}
1209
76ce4b5b 1210AliAODMCParticle* AliFemtoEventReaderAOD::GetParticleWithLabel(TClonesArray *mcP, Int_t aLabel)
1211{
1212 if (aLabel < 0) return 0;
1213 AliAODMCParticle *aodP;
1214 Int_t posstack = 0;
1215 if (aLabel > mcP->GetEntries())
1216 posstack = mcP->GetEntries();
1217 else
1218 posstack = aLabel;
1219
1220 aodP = (AliAODMCParticle *) mcP->At(posstack);
1221 if (aodP->GetLabel() > posstack) {
1222 do {
1223 aodP = (AliAODMCParticle *) mcP->At(posstack);
1224 if (aodP->GetLabel() == aLabel) return aodP;
1225 posstack--;
1226 }
1227 while (posstack > 0);
1228 }
1229 else {
1230 do {
1231 aodP = (AliAODMCParticle *) mcP->At(posstack);
1232 if (aodP->GetLabel() == aLabel) return aodP;
1233 posstack++;
1234 }
1235 while (posstack < mcP->GetEntries());
1236 }
87ede832 1237
76ce4b5b 1238 return 0;
1239}
1240
87ede832 1241void AliFemtoEventReaderAOD::CopyPIDtoFemtoTrack(AliAODTrack *tAodTrack,
1242 AliFemtoTrack *tFemtoTrack)
76ce4b5b 1243{
41cb6c1e 1244
50bb1be0 1245 if (fDCAglobalTrack) {
41cb6c1e 1246
906b97d0 1247 // code from Michael and Prabhat from AliAnalysisTaskDptDptCorrelations
1248 const AliAODVertex* vertex = (AliAODVertex*) fEvent->GetPrimaryVertex();
1249 float vertexX = -999.;
1250 float vertexY = -999.;
1251 float vertexZ = -999.;
1252
1253 if(vertex) {
1254 Double32_t fCov[6];
1255 vertex->GetCovarianceMatrix(fCov);
1256 if(vertex->GetNContributors() > 0) {
1257 if(fCov[5] != 0) {
1258 vertexX = vertex->GetX();
1259 vertexY = vertex->GetY();
1260 vertexZ = vertex->GetZ();
1261
1262 }
1263 }
1264 }
1265
1266 Double_t pos[3];
1267 tAodTrack->GetXYZ(pos);
6691ea4f 1268
906b97d0 1269 Double_t DCAX = pos[0] - vertexX;
1270 Double_t DCAY = pos[1] - vertexY;
1271 Double_t DCAZ = pos[2] - vertexZ;
e42e0efc 1272
906b97d0 1273 Double_t DCAXY = TMath::Sqrt((DCAX*DCAX) + (DCAY*DCAY));
6691ea4f 1274
906b97d0 1275 tFemtoTrack->SetImpactD(DCAXY);
1276 tFemtoTrack->SetImpactZ(DCAZ);
1277
1278
1279 // // // copying DCA information (taking it from global tracks gives better resolution than from TPC-only)
1280
1281 // double impact[2];
1282 // double covimpact[3];
1283
1284 // AliAODTrack* trk_clone = (AliAODTrack*)tAodTrack->Clone("trk_clone");
1285 // // if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
1286 // // if (!tAodTrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
1287 // if (!trk_clone->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
1288
1289 // // if (!tAodTrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
1290 // //cout << "sth went wrong with dca propagation" << endl;
1291 // tFemtoTrack->SetImpactD(-1000.0);
1292 // tFemtoTrack->SetImpactZ(-1000.0);
1293
1294 // }
1295 // else {
1296 // tFemtoTrack->SetImpactD(impact[0]);
1297 // tFemtoTrack->SetImpactZ(impact[1]);
1298 // }
1299 // delete trk_clone;
50bb1be0 1300 }
41cb6c1e 1301
76ce4b5b 1302 double aodpid[10];
1303 tAodTrack->GetPID(aodpid);
1304 tFemtoTrack->SetPidProbElectron(aodpid[0]);
1305 tFemtoTrack->SetPidProbMuon(aodpid[1]);
1306 tFemtoTrack->SetPidProbPion(aodpid[2]);
1307 tFemtoTrack->SetPidProbKaon(aodpid[3]);
1308 tFemtoTrack->SetPidProbProton(aodpid[4]);
1309
1310 aodpid[0] = -100000.0;
1311 aodpid[1] = -100000.0;
1312 aodpid[2] = -100000.0;
1313 aodpid[3] = -100000.0;
1314 aodpid[4] = -100000.0;
87ede832 1315
76ce4b5b 1316 double tTOF = 0.0;
87ede832 1317 Float_t probMis = 1.0;
76ce4b5b 1318
4eac0b05 1319 //what is that code? for what do we need it? nsigma values are not enough?
87ede832 1320 if (tAodTrack->GetStatus() & AliESDtrack::kTOFpid) { //AliESDtrack::kTOFpid=0x8000
1321 tTOF = tAodTrack->GetTOFsignal();
1322 tAodTrack->GetIntegratedTimes(aodpid);
1777446b 1323
87ede832 1324 tTOF -= fAODpidUtil->GetTOFResponse().GetStartTime(tAodTrack->P());
1325
1326 probMis = fAODpidUtil->GetTOFMismatchProbability(tAodTrack);
1327 }
1328
1329
1330
1331 tFemtoTrack->SetTofExpectedTimes(tTOF-aodpid[2], tTOF-aodpid[3], tTOF-aodpid[4]);
76ce4b5b 1332
76ce4b5b 1333 ////// TPC ////////////////////////////////////////////
1334
87ede832 1335 float nsigmaTPCK=-1000.;
1336 float nsigmaTPCPi=-1000.;
1337 float nsigmaTPCP=-1000.;
2e04885f 1338 float nsigmaTPCE=-1000.;
87ede832 1339
76ce4b5b 1340 // cout<<"in reader fESDpid"<<fESDpid<<endl;
1341
1342 if (tAodTrack->IsOn(AliESDtrack::kTPCpid)){ //AliESDtrack::kTPCpid=0x0080
1343 nsigmaTPCK = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kKaon);
1344 nsigmaTPCPi = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kPion);
1345 nsigmaTPCP = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kProton);
2e04885f 1346 nsigmaTPCE = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kElectron);
76ce4b5b 1347 }
1348
1349 tFemtoTrack->SetNSigmaTPCPi(nsigmaTPCPi);
1350 tFemtoTrack->SetNSigmaTPCK(nsigmaTPCK);
1351 tFemtoTrack->SetNSigmaTPCP(nsigmaTPCP);
2e04885f 1352 tFemtoTrack->SetNSigmaTPCE(nsigmaTPCE);
76ce4b5b 1353
87ede832 1354 tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF());
1355 tFemtoTrack->SetTPCncls(tAodTrack->GetTPCNcls());
1356 tFemtoTrack->SetTPCnclsF(tAodTrack->GetTPCNcls());
1357
1358 tFemtoTrack->SetTPCsignalN(1);
1359 tFemtoTrack->SetTPCsignalS(1);
76ce4b5b 1360 tFemtoTrack->SetTPCsignal(tAodTrack->GetTPCsignal());
76ce4b5b 1361
87ede832 1362 ///////TOF//////////////////////
76ce4b5b 1363
87ede832 1364 float vp=-1000.;
1365 float nsigmaTOFPi=-1000.;
1366 float nsigmaTOFK=-1000.;
1367 float nsigmaTOFP=-1000.;
2e04885f 1368 float nsigmaTOFE=-1000.;
87ede832 1369
1370 if ((tAodTrack->GetStatus() & AliESDtrack::kTOFpid) && //AliESDtrack::kTOFpid=0x8000
1371 (tAodTrack->GetStatus() & AliESDtrack::kTOFout) && //AliESDtrack::kTOFout=0x2000
1372 (tAodTrack->GetStatus() & AliESDtrack::kTIME) //AliESDtrack::kTIME=0x80000000
1373 && probMis < 0.01) // TOF mismatch probaility
1374 {
1375 if(tAodTrack->IsOn(AliESDtrack::kTOFpid)) //AliESDtrack::kTOFpid=0x8000
76ce4b5b 1376 {
1377
1378 nsigmaTOFPi = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kPion);
1379 nsigmaTOFK = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kKaon);
1380 nsigmaTOFP = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kProton);
2e04885f 1381 nsigmaTOFE = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kElectron);
76ce4b5b 1382
1383 Double_t len=200;// esdtrack->GetIntegratedLength(); !!!!!
1384 Double_t tof=tAodTrack->GetTOFsignal();
1385 if(tof > 0.) vp=len/tof/0.03;
1386 }
87ede832 1387 }
1388 tFemtoTrack->SetVTOF(vp);
1389 tFemtoTrack->SetNSigmaTOFPi(nsigmaTOFPi);
1390 tFemtoTrack->SetNSigmaTOFK(nsigmaTOFK);
1391 tFemtoTrack->SetNSigmaTOFP(nsigmaTOFP);
2e04885f 1392 tFemtoTrack->SetNSigmaTOFE(nsigmaTOFE);
87ede832 1393
76ce4b5b 1394
87ede832 1395 //////////////////////////////////////
76ce4b5b 1396
1397}
1398
1399void AliFemtoEventReaderAOD::SetCentralityPreSelection(double min, double max)
1400{
1401 fCentRange[0] = min; fCentRange[1] = max;
87ede832 1402 fUsePreCent = 1;
1777446b 1403 fEstEventMult = kCentrality;
76ce4b5b 1404}
1405
973a91f8 1406void AliFemtoEventReaderAOD::SetNoCentrality(bool anocent)
1407{
1777446b 1408 if(anocent==false) {fEstEventMult=kCentrality;}
1409 else {fEstEventMult=kReference; fUsePreCent = 0; }
973a91f8 1410}
76ce4b5b 1411
1412void AliFemtoEventReaderAOD::SetAODpidUtil(AliAODpidUtil *aAODpidUtil)
1413{
1414 fAODpidUtil = aAODpidUtil;
1415 // printf("fAODpidUtil: %x\n",fAODpidUtil);
1416}
1417
1777446b 1418void AliFemtoEventReaderAOD::SetAODheader(AliAODHeader *aAODheader)
1419{
1420 fAODheader = aAODheader;
1421}
1422
76ce4b5b 1423
ba3c23a4 1424void AliFemtoEventReaderAOD::SetMagneticFieldSign(int s)
1425{
1426 if(s>0)
1427 fMagFieldSign = 1;
1428 else if(s<0)
1429 fMagFieldSign = -1;
1430 else
1431 fMagFieldSign = 0;
1432}
1433
5e2038a9 1434void AliFemtoEventReaderAOD::SetEPVZERO(Bool_t iepvz)
1435{
87ede832 1436 fisEPVZ = iepvz;
5e2038a9 1437}
76ce4b5b 1438
ce7b3d98 1439void AliFemtoEventReaderAOD::GetGlobalPositionAtGlobalRadiiThroughTPC(AliAODTrack *track, Float_t bfield, Float_t globalPositionsAtRadii[9][3])
1440{
1441 // Gets the global position of the track at nine different radii in the TPC
1442 // track is the track you want to propagate
1443 // bfield is the magnetic field of your event
1444 // globalPositionsAtRadii is the array of global positions in the radii and xyz
1445
1446 // Initialize the array to something indicating there was no propagation
1447 for(Int_t i=0;i<9;i++){
1448 for(Int_t j=0;j<3;j++){
1449 globalPositionsAtRadii[i][j]=-9999.;
1450 }
1451 }
1452
1453 // Make a copy of the track to not change parameters of the track
1454 AliExternalTrackParam etp; etp.CopyFromVTrack(track);
1455 //printf("\nAfter CopyFromVTrack\n");
1456 //etp.Print();
1457
1458 // The global position of the the track
1459 Double_t xyz[3]={-9999.,-9999.,-9999.};
1460
1461 // Counter for which radius we want
1462 Int_t iR=0;
1463 // The radii at which we get the global positions
1464 // IROC (OROC) from 84.1 cm to 132.1 cm (134.6 cm to 246.6 cm)
1465 Float_t Rwanted[9]={85.,105.,125.,145.,165.,185.,205.,225.,245.};
1466 // The global radius we are at
1467 Float_t globalRadius=0;
1468
1469 // Propagation is done in local x of the track
1470 for (Float_t x = etp.GetX();x<247.;x+=1.){ // GetX returns local coordinates
1471 // Starts at the tracks fX and goes outwards. x = 245 is the outer radial limit
1472 // of the TPC when the track is straight, i.e. has inifinite pt and doesn't get bent.
1473 // If the track's momentum is smaller than infinite, it will develop a y-component, which
1474 // adds to the global radius
1475
1476 // Stop if the propagation was not succesful. This can happen for low pt tracks
1477 // that don't reach outer radii
1478 if(!etp.PropagateTo(x,bfield))break;
1479 etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1480 globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
1481
1482 // Roughly reached the radius we want
1483 if(globalRadius > Rwanted[iR]){
1484
1485 // Bigger loop has bad precision, we're nearly one centimeter too far, go back in small steps.
1486 while (globalRadius>Rwanted[iR]){
87ede832 1487 x-=.1;
1488 // printf("propagating to x %5.2f\n",x);
1489 if(!etp.PropagateTo(x,bfield))break;
1490 etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1491 globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
ce7b3d98 1492 }
1493 //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]);
1494 globalPositionsAtRadii[iR][0]=xyz[0];
1495 globalPositionsAtRadii[iR][1]=xyz[1];
1496 globalPositionsAtRadii[iR][2]=xyz[2];
1497 // Indicate we want the next radius
1498 iR+=1;
1499 }
1500 if(iR>=8){
1501 // TPC edge reached
1502 return;
1503 }
1504 }
1505}
76ce4b5b 1506
734c121a 1507void AliFemtoEventReaderAOD::SetpA2013(Bool_t pa2013)
1508{
1509 fpA2013 = pa2013;
1510}
50bb1be0 1511
1512void AliFemtoEventReaderAOD::SetDCAglobalTrack(Bool_t dcagt)
1513{
1514 fDCAglobalTrack = dcagt;
1515}
1516
316b08d6
MS		 1517
1518bool AliFemtoEventReaderAOD::RejectEventCentFlat(float MagField, float CentPercent)
1519{ // to flatten centrality distribution
1520 bool RejectEvent = kFALSE;
1521 int weightBinSign;
1522 TRandom3* fRandomNumber = new TRandom3(); //for 3D, random sign switching
1523 fRandomNumber->SetSeed(0);
1524
1525 if(MagField > 0) weightBinSign = 0;
1526 else weightBinSign = 1;
1527 float kCentWeight[2][9] = {{.878,.876,.860,.859,.859,.88,.873,.879,.894},
1528 {.828,.793,.776,.772,.775,.796,.788,.804,.839}};
1529 int weightBinCent = (int) CentPercent;
1530 if(fRandomNumber->Rndm() > kCentWeight[weightBinSign][weightBinCent]) RejectEvent = kTRUE;
1531
1532 return RejectEvent;
1533}
1534
1535void AliFemtoEventReaderAOD::SetCentralityFlattening(Bool_t dcagt)
1536{
1537 fFlatCent = dcagt;
1538}