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