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