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