]> 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
Adding TOF calib task for calibration of problematic channels
[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;
3a817989 831 AliFemtoV0* trackCopyV0 = CopyAODtoFemtoV0(aodv0);
87ede832 832 tEvent->V0Collection()->push_back(trackCopyV0);
3a817989 833 count_pass++;
87ede832 834 //cout<<"Pushback v0 to v0collection"<<endl;
973a91f8 835 }
87ede832 836 }
973a91f8 837
ae746ffd 838 return tEvent;
76ce4b5b 839}
840
ae746ffd 841AliFemtoTrack* AliFemtoEventReaderAOD::CopyAODtoFemtoTrack(AliAODTrack *tAodTrack
87ede832 842 // AliPWG2AODTrack *tPWG2AODTrack
843 )
76ce4b5b 844{
845 // Copy the track information from the AOD into the internal AliFemtoTrack
846 // If it exists, use the additional information from the PWG2 AOD
ae746ffd 847 AliFemtoTrack *tFemtoTrack = new AliFemtoTrack();
76ce4b5b 848
849 // Primary Vertex position
87ede832 850
76ce4b5b 851 fEvent->GetPrimaryVertex()->GetPosition(fV1);
852 // fEvent->GetPrimaryVertex()->GetXYZ(fV1);
eef80284 853 tFemtoTrack->SetPrimaryVertex(fV1);
76ce4b5b 854
855 tFemtoTrack->SetCharge(tAodTrack->Charge());
87ede832 856
76ce4b5b 857 double pxyz[3];
858 tAodTrack->PxPyPz(pxyz);//reading noconstrained momentum
290ffcb9 859
76ce4b5b 860 AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
861 tFemtoTrack->SetP(v);//setting momentum
862 tFemtoTrack->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
863 const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
87ede832 864 //setting track helix
76ce4b5b 865 const AliFmThreeVectorD ktP(pxyz[0],pxyz[1],pxyz[2]);
87ede832 866 AliFmPhysicalHelixD helix(ktP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(tFemtoTrack->Charge()));
76ce4b5b 867 tFemtoTrack->SetHelix(helix);
87ede832 868
76ce4b5b 869 // Flags
870 tFemtoTrack->SetTrackId(tAodTrack->GetID());
871 tFemtoTrack->SetFlags(tAodTrack->GetFlags());
872 tFemtoTrack->SetLabel(tAodTrack->GetLabel());
87ede832 873
874 // Track quality information
76ce4b5b 875 float covmat[6];
87ede832 876 tAodTrack->GetCovMatrix(covmat);
877
50bb1be0 878 if (!fDCAglobalTrack) {
76ce4b5b 879
906b97d0 880 tFemtoTrack->SetImpactD(tAodTrack->DCA());
881 tFemtoTrack->SetImpactZ(tAodTrack->ZAtDCA());
882
76ce4b5b 883
ef838a91 884 // double impact[2];
885 // double covimpact[3];
886
887 // AliAODTrack* trk_clone = (AliAODTrack*)tAodTrack->Clone("trk_clone");
888 // // if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
889 // // if (!tAodTrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
890 // if (!trk_clone->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
891 // //cout << "sth went wrong with dca propagation" << endl;
892 // tFemtoTrack->SetImpactD(-1000.0);
893 // tFemtoTrack->SetImpactZ(-1000.0);
894
895 // }
896 // else {
897 // tFemtoTrack->SetImpactD(impact[0]);
898 // tFemtoTrack->SetImpactZ(impact[1]);
899 // }
900 // delete trk_clone;
e42e0efc 901
50bb1be0 902 }
76ce4b5b 903
904 // if (TMath::Abs(tAodTrack->Xv()) > 0.00000000001)
905 // tFemtoTrack->SetImpactD(TMath::Hypot(tAodTrack->Xv(), tAodTrack->Yv())*(tAodTrack->Xv()/TMath::Abs(tAodTrack->Xv())));
906 // else
907 // tFemtoTrack->SetImpactD(0.0);
908 // tFemtoTrack->SetImpactD(tAodTrack->DCA());
87ede832 909
76ce4b5b 910 // tFemtoTrack->SetImpactZ(tAodTrack->ZAtDCA());
911
912
913 // tFemtoTrack->SetImpactD(TMath::Hypot(tAodTrack->Xv() - fV1[0], tAodTrack->Yv() - fV1[1]));
914 // tFemtoTrack->SetImpactZ(tAodTrack->Zv() - fV1[2]);
915
916
87ede832 917 // cout
918 // << "dca" << TMath::Hypot(tAodTrack->Xv() - fV1[0], tAodTrack->Yv() - fV1[1])
919 // << "xv - fv10 = "<< tAodTrack->Xv() - fV1[0]
920 // << tAodTrack->Yv() - fV1[1]
921// << "xv = " << tAodTrack->Xv() << endl
922// << "fv1[0] = " << fV1[0] << endl
923// << "yv = " << tAodTrack->Yv() << endl
924// << "fv1[1] = " << fV1[1] << endl
925// << "zv = " << tAodTrack->Zv() << endl
926// << "fv1[2] = " << fV1[2] << endl
927// << "impact[0] = " << impact[0] << endl
928// << "impact[1] = " << impact[1] << endl
76ce4b5b 929// << endl << endl ;
930
931 tFemtoTrack->SetCdd(covmat[0]);
932 tFemtoTrack->SetCdz(covmat[1]);
933 tFemtoTrack->SetCzz(covmat[2]);
973a91f8 934 tFemtoTrack->SetITSchi2(tAodTrack->Chi2perNDF());
935 tFemtoTrack->SetITSncls(tAodTrack->GetITSNcls());
936 tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF());
937 tFemtoTrack->SetTPCncls(tAodTrack->GetTPCNcls());
938 tFemtoTrack->SetTPCnclsF(tAodTrack->GetTPCNcls());
87ede832 939 tFemtoTrack->SetTPCsignalN(1);
940 tFemtoTrack->SetTPCsignalS(1);
76ce4b5b 941 tFemtoTrack->SetTPCsignal(tAodTrack->GetTPCsignal());
942
943// if (tPWG2AODTrack) {
944// // Copy the PWG2 specific information if it exists
945// tFemtoTrack->SetTPCClusterMap(tPWG2AODTrack->GetTPCClusterMap());
946// tFemtoTrack->SetTPCSharedMap(tPWG2AODTrack->GetTPCSharedMap());
87ede832 947
76ce4b5b 948// double xtpc[3] = {0,0,0};
949// tPWG2AODTrack->GetTPCNominalEntrancePoint(xtpc);
950// tFemtoTrack->SetNominalTPCEntrancePoint(xtpc);
951// tPWG2AODTrack->GetTPCNominalExitPoint(xtpc);
952// tFemtoTrack->SetNominalTPCExitPoint(xtpc);
953// }
954// else {
87ede832 955 // If not use dummy values
76ce4b5b 956 tFemtoTrack->SetTPCClusterMap(tAodTrack->GetTPCClusterMap());
957 tFemtoTrack->SetTPCSharedMap(tAodTrack->GetTPCSharedMap());
87ede832 958
ce7b3d98 959 float globalPositionsAtRadii[9][3];
960 float bfield = 5*fMagFieldSign;
c710fd31 961
ce7b3d98 962 GetGlobalPositionAtGlobalRadiiThroughTPC(tAodTrack,bfield,globalPositionsAtRadii);
963 double tpcEntrance[3]={globalPositionsAtRadii[0][0],globalPositionsAtRadii[0][1],globalPositionsAtRadii[0][2]};
964 double **tpcPositions;
965 tpcPositions = new double*[9];
966 for(int i=0;i<9;i++)
967 tpcPositions[i] = new double[3];
968 double tpcExit[3]={globalPositionsAtRadii[8][0],globalPositionsAtRadii[8][1],globalPositionsAtRadii[8][2]};
969 for(int i=0;i<9;i++)
87ede832 970 {
971 tpcPositions[i][0] = globalPositionsAtRadii[i][0];
972 tpcPositions[i][1] = globalPositionsAtRadii[i][1];
973 tpcPositions[i][2] = globalPositionsAtRadii[i][2];
974 }
7db4e4cf 975
ce7b3d98 976 tFemtoTrack->SetNominalTPCEntrancePoint(tpcEntrance);
977 tFemtoTrack->SetNominalTPCPoints(tpcPositions);
978 tFemtoTrack->SetNominalTPCExitPoint(tpcExit);
87ede832 979 for(int i=0;i<9;i++)
980 delete [] tpcPositions[i];
981 delete [] tpcPositions;
76ce4b5b 982 // }
87ede832 983
76ce4b5b 984 // // cout << "Track has " << TMath::Hypot(tAodTrack->Xv(), tAodTrack->Yv()) << " " << tAodTrack->Zv() << " " << tAodTrack->GetTPCNcls() << endl;
87ede832 985
986
76ce4b5b 987 int indexes[3];
988 for (int ik=0; ik<3; ik++) {
989 indexes[ik] = 0;
990 }
991 tFemtoTrack->SetKinkIndexes(indexes);
4eac0b05 992
993
994 for (int ii=0; ii<6; ii++){
995 tFemtoTrack->SetITSHitOnLayer(ii,tAodTrack->HasPointOnITSLayer(ii));
996 }
997
ae746ffd 998 return tFemtoTrack;
76ce4b5b 999}
1000
cbc1c949 1001AliFemtoV0* AliFemtoEventReaderAOD::CopyAODtoFemtoV0(AliAODv0 *tAODv0 )
973a91f8 1002{
cbc1c949 1003 AliFemtoV0 *tFemtoV0 = new AliFemtoV0();
1004
ce7b3d98 1005 tFemtoV0->SetdecayLengthV0(tAODv0->DecayLength(fV1));
973a91f8 1006 tFemtoV0->SetdecayVertexV0X(tAODv0->DecayVertexV0X());
1007 tFemtoV0->SetdecayVertexV0Y(tAODv0->DecayVertexV0Y());
1008 tFemtoV0->SetdecayVertexV0Z(tAODv0->DecayVertexV0Z());
1009 AliFemtoThreeVector decayvertex(tAODv0->DecayVertexV0X(),tAODv0->DecayVertexV0Y(),tAODv0->DecayVertexV0Z());
1010 tFemtoV0->SetdecayVertexV0(decayvertex);
1011 tFemtoV0->SetdcaV0Daughters(tAODv0->DcaV0Daughters());
1012 tFemtoV0->SetdcaV0ToPrimVertex(tAODv0->DcaV0ToPrimVertex());
1013 tFemtoV0->SetdcaPosToPrimVertex(tAODv0->DcaPosToPrimVertex());
1014 tFemtoV0->SetdcaNegToPrimVertex(tAODv0->DcaNegToPrimVertex());
1015 tFemtoV0->SetmomPosX(tAODv0->MomPosX());
1016 tFemtoV0->SetmomPosY(tAODv0->MomPosY());
1017 tFemtoV0->SetmomPosZ(tAODv0->MomPosZ());
1018 AliFemtoThreeVector mompos(tAODv0->MomPosX(),tAODv0->MomPosY(),tAODv0->MomPosZ());
1019 tFemtoV0->SetmomPos(mompos);
1020 tFemtoV0->SetmomNegX(tAODv0->MomNegX());
1021 tFemtoV0->SetmomNegY(tAODv0->MomNegY());
1022 tFemtoV0->SetmomNegZ(tAODv0->MomNegZ());
1023 AliFemtoThreeVector momneg(tAODv0->MomNegX(),tAODv0->MomNegY(),tAODv0->MomNegZ());
1024 tFemtoV0->SetmomNeg(momneg);
1025
1026 //jest cos takiego w AliFemtoV0.h czego nie ma w AliAODv0.h
87ede832 1027 //void SettpcHitsPos(const int& i);
1028 //void SettpcHitsNeg(const int& i);
973a91f8 1029
1030 //void SetTrackTopologyMapPos(unsigned int word, const unsigned long& m);
1031 //void SetTrackTopologyMapNeg(unsigned int word, const unsigned long& m);
1032
cbc1c949 1033
973a91f8 1034 tFemtoV0->SetmomV0X(tAODv0->MomV0X());
1035 tFemtoV0->SetmomV0Y(tAODv0->MomV0Y());
1036 tFemtoV0->SetmomV0Z(tAODv0->MomV0Z());
1037 AliFemtoThreeVector momv0(tAODv0->MomV0X(),tAODv0->MomV0Y(),tAODv0->MomV0Z());
1038 tFemtoV0->SetmomV0(momv0);
1039 tFemtoV0->SetalphaV0(tAODv0->AlphaV0());
1040 tFemtoV0->SetptArmV0(tAODv0->PtArmV0());
1041 tFemtoV0->SeteLambda(tAODv0->ELambda());
1042 tFemtoV0->SeteK0Short(tAODv0->EK0Short());
1043 tFemtoV0->SetePosProton(tAODv0->EPosProton());
1044 tFemtoV0->SeteNegProton(tAODv0->ENegProton());
1045 tFemtoV0->SetmassLambda(tAODv0->MassLambda());
1046 tFemtoV0->SetmassAntiLambda(tAODv0->MassAntiLambda());
1047 tFemtoV0->SetmassK0Short(tAODv0->MassK0Short());
1048 tFemtoV0->SetrapLambda(tAODv0->RapLambda());
1049 tFemtoV0->SetrapK0Short(tAODv0->RapK0Short());
87ede832 1050
1051 //void SetcTauLambda( float x);
1052 //void SetcTauK0Short( float x);
1053
ce7b3d98 1054 //tFemtoV0->SetptV0(::sqrt(tAODv0->Pt2V0())); //!
1055 tFemtoV0->SetptV0(tAODv0->Pt());
973a91f8 1056 tFemtoV0->SetptotV0(::sqrt(tAODv0->Ptot2V0()));
ce7b3d98 1057 //tFemtoV0->SetptPos(::sqrt(tAODv0->MomPosX()*tAODv0->MomPosX()+tAODv0->MomPosY()*tAODv0->MomPosY()));
1058 //tFemtoV0->SetptotPos(::sqrt(tAODv0->Ptot2Pos()));
1059 //tFemtoV0->SetptNeg(::sqrt(tAODv0->MomNegX()*tAODv0->MomNegX()+tAODv0->MomNegY()*tAODv0->MomNegY()));
1060 //tFemtoV0->SetptotNeg(::sqrt(tAODv0->Ptot2Neg()));
87ede832 1061
973a91f8 1062 tFemtoV0->SetidNeg(tAODv0->GetNegID());
1063 //cout<<"tAODv0->GetNegID(): "<<tAODv0->GetNegID()<<endl;
1064 //cout<<"tFemtoV0->IdNeg(): "<<tFemtoV0->IdNeg()<<endl;
1065 tFemtoV0->SetidPos(tAODv0->GetPosID());
1066
1067 tFemtoV0->SetEtaV0(tAODv0->Eta());
1068 tFemtoV0->SetPhiV0(tAODv0->Phi());
1069 tFemtoV0->SetCosPointingAngle(tAODv0->CosPointingAngle(fV1));
1070 //tFemtoV0->SetYV0(tAODv0->Y());
1071
cbc1c949 1072
973a91f8 1073 //void SetdedxNeg(float x);
1074 //void SeterrdedxNeg(float x);//Gael 04Fev2002
1075 //void SetlendedxNeg(float x);//Gael 04Fev2002
1076 //void SetdedxPos(float x);
1077 //void SeterrdedxPos(float x);//Gael 04Fev2002
1078 //void SetlendedxPos(float x);//Gael 04Fev2002
1079
ce7b3d98 1080 //tFemtoV0->SetEtaPos(tAODv0->PseudoRapPos());
1081 //tFemtoV0->SetEtaNeg(tAODv0->PseudoRapNeg());
973a91f8 1082
1083 AliAODTrack *trackpos = (AliAODTrack*)tAODv0->GetDaughter(0);
1084 AliAODTrack *trackneg = (AliAODTrack*)tAODv0->GetDaughter(1);
1085
cbc1c949 1086
973a91f8 1087 if(trackpos && trackneg)
87ede832 1088 {
1089 tFemtoV0->SetEtaPos(trackpos->Eta());
1090 tFemtoV0->SetEtaNeg(trackneg->Eta());
1091 tFemtoV0->SetptotPos(tAODv0->PProng(0));
1092 tFemtoV0->SetptotNeg(tAODv0->PProng(1));
1093 tFemtoV0->SetptPos(trackpos->Pt());
1094 tFemtoV0->SetptNeg(trackneg->Pt());
1095
1096 //tFemtoV0->SetEtaPos(trackpos->Eta()); //tAODv0->PseudoRapPos()
1097 //tFemtoV0->SetEtaNeg(trackneg->Eta()); //tAODv0->PseudoRapNeg()
1098 tFemtoV0->SetTPCNclsPos(trackpos->GetTPCNcls());
1099 tFemtoV0->SetTPCNclsNeg(trackneg->GetTPCNcls());
1100 tFemtoV0->SetTPCclustersPos(trackpos->GetTPCClusterMap());
1101 tFemtoV0->SetTPCclustersNeg(trackneg->GetTPCClusterMap());
1102 tFemtoV0->SetTPCsharingPos(trackpos->GetTPCSharedMap());
1103 tFemtoV0->SetTPCsharingNeg(trackneg->GetTPCSharedMap());
1104 tFemtoV0->SetNdofPos(trackpos->Chi2perNDF());
1105 tFemtoV0->SetNdofNeg(trackneg->Chi2perNDF());
1106 tFemtoV0->SetStatusPos(trackpos->GetStatus());
1107 tFemtoV0->SetStatusNeg(trackneg->GetStatus());
1108
1109 tFemtoV0->SetPosNSigmaTPCK(fAODpidUtil->NumberOfSigmasTPC(trackpos,AliPID::kKaon));
1110 tFemtoV0->SetNegNSigmaTPCK(fAODpidUtil->NumberOfSigmasTPC(trackneg,AliPID::kKaon));
1111 tFemtoV0->SetPosNSigmaTPCP(fAODpidUtil->NumberOfSigmasTPC(trackpos,AliPID::kProton));
1112 tFemtoV0->SetNegNSigmaTPCP(fAODpidUtil->NumberOfSigmasTPC(trackneg,AliPID::kProton));
1113 tFemtoV0->SetPosNSigmaTPCPi(fAODpidUtil->NumberOfSigmasTPC(trackpos,AliPID::kPion));
1114 tFemtoV0->SetNegNSigmaTPCPi(fAODpidUtil->NumberOfSigmasTPC(trackneg,AliPID::kPion));
1115
1116
1117 float bfield = 5*fMagFieldSign;
1118 float globalPositionsAtRadiiPos[9][3];
1119 GetGlobalPositionAtGlobalRadiiThroughTPC(trackpos,bfield,globalPositionsAtRadiiPos);
1120 double tpcEntrancePos[3]={globalPositionsAtRadiiPos[0][0],globalPositionsAtRadiiPos[0][1],globalPositionsAtRadiiPos[0][2]};
1121 double tpcExitPos[3]={globalPositionsAtRadiiPos[8][0],globalPositionsAtRadiiPos[8][1],globalPositionsAtRadiiPos[8][2]};
1122
1123 float globalPositionsAtRadiiNeg[9][3];
1124 GetGlobalPositionAtGlobalRadiiThroughTPC(trackneg,bfield,globalPositionsAtRadiiNeg);
1125 double tpcEntranceNeg[3]={globalPositionsAtRadiiNeg[0][0],globalPositionsAtRadiiNeg[0][1],globalPositionsAtRadiiNeg[0][2]};
1126 double tpcExitNeg[3]={globalPositionsAtRadiiNeg[8][0],globalPositionsAtRadiiNeg[8][1],globalPositionsAtRadiiNeg[8][2]};
1127
1128 AliFemtoThreeVector tmpVec;
290ffcb9 1129 tmpVec.SetX(tpcEntrancePos[0]); tmpVec.SetY(tpcEntrancePos[1]); tmpVec.SetZ(tpcEntrancePos[2]);
87ede832 1130 tFemtoV0->SetNominalTpcEntrancePointPos(tmpVec);
1131
290ffcb9 1132 tmpVec.SetX(tpcExitPos[0]); tmpVec.SetY(tpcExitPos[1]); tmpVec.SetZ(tpcExitPos[2]);
87ede832 1133 tFemtoV0->SetNominalTpcExitPointPos(tmpVec);
1134
290ffcb9 1135 tmpVec.SetX(tpcEntranceNeg[0]); tmpVec.SetY(tpcEntranceNeg[1]); tmpVec.SetZ(tpcEntranceNeg[2]);
87ede832 1136 tFemtoV0->SetNominalTpcEntrancePointNeg(tmpVec);
1137
290ffcb9 1138 tmpVec.SetX(tpcExitNeg[0]); tmpVec.SetY(tpcExitNeg[1]); tmpVec.SetZ(tpcExitNeg[2]);
87ede832 1139 tFemtoV0->SetNominalTpcExitPointNeg(tmpVec);
1140
cbc1c949 1141
87ede832 1142 AliFemtoThreeVector vecTpcPos[9];
1143 AliFemtoThreeVector vecTpcNeg[9];
1144 for(int i=0;i<9;i++)
973a91f8 1145 {
87ede832 1146 vecTpcPos[i].SetX(globalPositionsAtRadiiPos[i][0]); vecTpcPos[i].SetY(globalPositionsAtRadiiPos[i][1]); vecTpcPos[i].SetZ(globalPositionsAtRadiiPos[i][2]);
1147 vecTpcNeg[i].SetX(globalPositionsAtRadiiNeg[i][0]); vecTpcNeg[i].SetY(globalPositionsAtRadiiNeg[i][1]); vecTpcNeg[i].SetZ(globalPositionsAtRadiiNeg[i][2]);
1148 }
1149 tFemtoV0->SetNominalTpcPointPos(vecTpcPos);
1150 tFemtoV0->SetNominalTpcPointNeg(vecTpcNeg);
ce7b3d98 1151
87ede832 1152 tFemtoV0->SetTPCMomentumPos(trackpos->GetTPCmomentum());
1153 tFemtoV0->SetTPCMomentumNeg(trackneg->GetTPCmomentum());
ce7b3d98 1154
87ede832 1155 tFemtoV0->SetdedxPos(trackpos->GetTPCsignal());
1156 tFemtoV0->SetdedxNeg(trackneg->GetTPCsignal());
ce7b3d98 1157
b6200829 1158
1159 Float_t probMisPos = 1.0;
1160 Float_t probMisNeg = 1.0;
1161
1162 if (tFemtoV0->StatusPos() & AliESDtrack::kTOFpid) { //AliESDtrack::kTOFpid=0x8000
1163 probMisPos = fAODpidUtil->GetTOFMismatchProbability(trackpos);
1164 }
1165 if (tFemtoV0->StatusNeg() & AliESDtrack::kTOFpid) { //AliESDtrack::kTOFpid=0x8000
1166 probMisNeg = fAODpidUtil->GetTOFMismatchProbability(trackneg);
1167 }
1168
1169 if((tFemtoV0->StatusPos()&AliESDtrack::kTOFpid)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTOFout)==0 || probMisPos > 0.01)
87ede832 1170 {
b6200829 1171 if((tFemtoV0->StatusNeg()&AliESDtrack::kTOFpid)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTOFout)==0 || probMisNeg > 0.01)
973a91f8 1172 {
ce7b3d98 1173 tFemtoV0->SetPosNSigmaTOFK(-1000);
1174 tFemtoV0->SetNegNSigmaTOFK(-1000);
1175 tFemtoV0->SetPosNSigmaTOFP(-1000);
1176 tFemtoV0->SetNegNSigmaTOFP(-1000);
1177 tFemtoV0->SetPosNSigmaTOFPi(-1000);
1178 tFemtoV0->SetNegNSigmaTOFPi(-1000);
1179
1180 tFemtoV0->SetTOFProtonTimePos(-1000);
1181 tFemtoV0->SetTOFPionTimePos(-1000);
1182 tFemtoV0->SetTOFKaonTimePos(-1000);
1183 tFemtoV0->SetTOFProtonTimeNeg(-1000);
1184 tFemtoV0->SetTOFPionTimeNeg(-1000);
1185 tFemtoV0->SetTOFKaonTimeNeg(-1000);
973a91f8 1186 }
973a91f8 1187 }
87ede832 1188 else
973a91f8 1189 {
b6200829 1190 if(trackpos->IsOn(AliESDtrack::kTOFpid)) {
1191 tFemtoV0->SetPosNSigmaTOFK(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kKaon));
1192 tFemtoV0->SetPosNSigmaTOFP(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kProton));
1193 tFemtoV0->SetPosNSigmaTOFPi(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kPion));
1194 }
1195 if(trackneg->IsOn(AliESDtrack::kTOFpid)) {
1196 tFemtoV0->SetNegNSigmaTOFK(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kKaon));
1197 tFemtoV0->SetNegNSigmaTOFP(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kProton));
1198 tFemtoV0->SetNegNSigmaTOFPi(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kPion));
1199 }
87ede832 1200 double TOFSignalPos = trackpos->GetTOFsignal();
1201 double TOFSignalNeg = trackneg->GetTOFsignal();
b6200829 1202 TOFSignalPos -= fAODpidUtil->GetTOFResponse().GetStartTime(trackpos->P());
1203 TOFSignalNeg -= fAODpidUtil->GetTOFResponse().GetStartTime(trackneg->P());
87ede832 1204 double pidPos[5];
1205 double pidNeg[5];
1206 trackpos->GetIntegratedTimes(pidPos);
1207 trackneg->GetIntegratedTimes(pidNeg);
1208
1209 tFemtoV0->SetTOFPionTimePos(TOFSignalPos-pidPos[2]);
1210 tFemtoV0->SetTOFKaonTimePos(TOFSignalPos-pidPos[3]);
1211 tFemtoV0->SetTOFProtonTimePos(TOFSignalPos-pidPos[4]);
1212 tFemtoV0->SetTOFPionTimeNeg(TOFSignalNeg-pidNeg[2]);
1213 tFemtoV0->SetTOFKaonTimeNeg(TOFSignalNeg-pidNeg[3]);
1214 tFemtoV0->SetTOFProtonTimeNeg(TOFSignalNeg-pidNeg[4]);
973a91f8 1215 }
cbc1c949 1216
87ede832 1217 }
1218 else
1219 {
cbc1c949 1220
87ede832 1221 tFemtoV0->SetStatusPos(999);
1222 tFemtoV0->SetStatusNeg(999);
1223 }
cbc1c949 1224
973a91f8 1225 tFemtoV0->SetOnFlyStatusV0(tAODv0->GetOnFlyStatus());
cbc1c949 1226 return tFemtoV0;
973a91f8 1227}
1228
76ce4b5b 1229void AliFemtoEventReaderAOD::SetFilterBit(UInt_t ibit)
1230{
1231 fFilterBit = (1 << (ibit));
1232}
1233
64536eaf 1234
1235void AliFemtoEventReaderAOD::SetFilterMask(int ibit)
1236{
1237 fFilterMask = ibit;
1238}
1239
76ce4b5b 1240void AliFemtoEventReaderAOD::SetReadMC(unsigned char a)
1241{
1242 fReadMC = a;
1243}
1244
973a91f8 1245
1246void AliFemtoEventReaderAOD::SetReadV0(unsigned char a)
1247{
1248 fReadV0 = a;
1249}
1250
1777446b 1251void AliFemtoEventReaderAOD::SetUseMultiplicity(EstEventMult aType)
1252{
1253 fEstEventMult = aType;
1254}
1255
76ce4b5b 1256AliAODMCParticle* AliFemtoEventReaderAOD::GetParticleWithLabel(TClonesArray *mcP, Int_t aLabel)
1257{
1258 if (aLabel < 0) return 0;
1259 AliAODMCParticle *aodP;
1260 Int_t posstack = 0;
1261 if (aLabel > mcP->GetEntries())
1262 posstack = mcP->GetEntries();
1263 else
1264 posstack = aLabel;
1265
1266 aodP = (AliAODMCParticle *) mcP->At(posstack);
1267 if (aodP->GetLabel() > posstack) {
1268 do {
1269 aodP = (AliAODMCParticle *) mcP->At(posstack);
1270 if (aodP->GetLabel() == aLabel) return aodP;
1271 posstack--;
1272 }
1273 while (posstack > 0);
1274 }
1275 else {
1276 do {
1277 aodP = (AliAODMCParticle *) mcP->At(posstack);
1278 if (aodP->GetLabel() == aLabel) return aodP;
1279 posstack++;
1280 }
1281 while (posstack < mcP->GetEntries());
1282 }
87ede832 1283
76ce4b5b 1284 return 0;
1285}
1286
87ede832 1287void AliFemtoEventReaderAOD::CopyPIDtoFemtoTrack(AliAODTrack *tAodTrack,
1288 AliFemtoTrack *tFemtoTrack)
76ce4b5b 1289{
41cb6c1e 1290
50bb1be0 1291 if (fDCAglobalTrack) {
41cb6c1e 1292
906b97d0 1293 // code from Michael and Prabhat from AliAnalysisTaskDptDptCorrelations
1294 const AliAODVertex* vertex = (AliAODVertex*) fEvent->GetPrimaryVertex();
1295 float vertexX = -999.;
1296 float vertexY = -999.;
1297 float vertexZ = -999.;
1298
1299 if(vertex) {
1300 Double32_t fCov[6];
1301 vertex->GetCovarianceMatrix(fCov);
1302 if(vertex->GetNContributors() > 0) {
1303 if(fCov[5] != 0) {
1304 vertexX = vertex->GetX();
1305 vertexY = vertex->GetY();
1306 vertexZ = vertex->GetZ();
1307
1308 }
1309 }
1310 }
1311
1312 Double_t pos[3];
1313 tAodTrack->GetXYZ(pos);
6691ea4f 1314
906b97d0 1315 Double_t DCAX = pos[0] - vertexX;
1316 Double_t DCAY = pos[1] - vertexY;
1317 Double_t DCAZ = pos[2] - vertexZ;
e42e0efc 1318
906b97d0 1319 Double_t DCAXY = TMath::Sqrt((DCAX*DCAX) + (DCAY*DCAY));
6691ea4f 1320
906b97d0 1321 tFemtoTrack->SetImpactD(DCAXY);
1322 tFemtoTrack->SetImpactZ(DCAZ);
1323
1324
1325 // // // copying DCA information (taking it from global tracks gives better resolution than from TPC-only)
1326
1327 // double impact[2];
1328 // double covimpact[3];
1329
1330 // AliAODTrack* trk_clone = (AliAODTrack*)tAodTrack->Clone("trk_clone");
1331 // // if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
1332 // // if (!tAodTrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
1333 // if (!trk_clone->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
1334
1335 // // if (!tAodTrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
1336 // //cout << "sth went wrong with dca propagation" << endl;
1337 // tFemtoTrack->SetImpactD(-1000.0);
1338 // tFemtoTrack->SetImpactZ(-1000.0);
1339
1340 // }
1341 // else {
1342 // tFemtoTrack->SetImpactD(impact[0]);
1343 // tFemtoTrack->SetImpactZ(impact[1]);
1344 // }
1345 // delete trk_clone;
50bb1be0 1346 }
41cb6c1e 1347
76ce4b5b 1348 double aodpid[10];
1349 tAodTrack->GetPID(aodpid);
1350 tFemtoTrack->SetPidProbElectron(aodpid[0]);
1351 tFemtoTrack->SetPidProbMuon(aodpid[1]);
1352 tFemtoTrack->SetPidProbPion(aodpid[2]);
1353 tFemtoTrack->SetPidProbKaon(aodpid[3]);
1354 tFemtoTrack->SetPidProbProton(aodpid[4]);
1355
1356 aodpid[0] = -100000.0;
1357 aodpid[1] = -100000.0;
1358 aodpid[2] = -100000.0;
1359 aodpid[3] = -100000.0;
1360 aodpid[4] = -100000.0;
87ede832 1361
76ce4b5b 1362 double tTOF = 0.0;
87ede832 1363 Float_t probMis = 1.0;
76ce4b5b 1364
4eac0b05 1365 //what is that code? for what do we need it? nsigma values are not enough?
87ede832 1366 if (tAodTrack->GetStatus() & AliESDtrack::kTOFpid) { //AliESDtrack::kTOFpid=0x8000
1367 tTOF = tAodTrack->GetTOFsignal();
1368 tAodTrack->GetIntegratedTimes(aodpid);
1777446b 1369
87ede832 1370 tTOF -= fAODpidUtil->GetTOFResponse().GetStartTime(tAodTrack->P());
1371
1372 probMis = fAODpidUtil->GetTOFMismatchProbability(tAodTrack);
1373 }
1374
1375
1376
1377 tFemtoTrack->SetTofExpectedTimes(tTOF-aodpid[2], tTOF-aodpid[3], tTOF-aodpid[4]);
76ce4b5b 1378
76ce4b5b 1379 ////// TPC ////////////////////////////////////////////
1380
87ede832 1381 float nsigmaTPCK=-1000.;
1382 float nsigmaTPCPi=-1000.;
1383 float nsigmaTPCP=-1000.;
2e04885f 1384 float nsigmaTPCE=-1000.;
87ede832 1385
76ce4b5b 1386 // cout<<"in reader fESDpid"<<fESDpid<<endl;
1387
1388 if (tAodTrack->IsOn(AliESDtrack::kTPCpid)){ //AliESDtrack::kTPCpid=0x0080
1389 nsigmaTPCK = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kKaon);
1390 nsigmaTPCPi = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kPion);
1391 nsigmaTPCP = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kProton);
2e04885f 1392 nsigmaTPCE = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kElectron);
76ce4b5b 1393 }
1394
1395 tFemtoTrack->SetNSigmaTPCPi(nsigmaTPCPi);
1396 tFemtoTrack->SetNSigmaTPCK(nsigmaTPCK);
1397 tFemtoTrack->SetNSigmaTPCP(nsigmaTPCP);
2e04885f 1398 tFemtoTrack->SetNSigmaTPCE(nsigmaTPCE);
76ce4b5b 1399
87ede832 1400 tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF());
1401 tFemtoTrack->SetTPCncls(tAodTrack->GetTPCNcls());
1402 tFemtoTrack->SetTPCnclsF(tAodTrack->GetTPCNcls());
1403
1404 tFemtoTrack->SetTPCsignalN(1);
1405 tFemtoTrack->SetTPCsignalS(1);
76ce4b5b 1406 tFemtoTrack->SetTPCsignal(tAodTrack->GetTPCsignal());
76ce4b5b 1407
87ede832 1408 ///////TOF//////////////////////
76ce4b5b 1409
87ede832 1410 float vp=-1000.;
1411 float nsigmaTOFPi=-1000.;
1412 float nsigmaTOFK=-1000.;
1413 float nsigmaTOFP=-1000.;
2e04885f 1414 float nsigmaTOFE=-1000.;
87ede832 1415
1416 if ((tAodTrack->GetStatus() & AliESDtrack::kTOFpid) && //AliESDtrack::kTOFpid=0x8000
1417 (tAodTrack->GetStatus() & AliESDtrack::kTOFout) && //AliESDtrack::kTOFout=0x2000
1418 (tAodTrack->GetStatus() & AliESDtrack::kTIME) //AliESDtrack::kTIME=0x80000000
1419 && probMis < 0.01) // TOF mismatch probaility
1420 {
1421 if(tAodTrack->IsOn(AliESDtrack::kTOFpid)) //AliESDtrack::kTOFpid=0x8000
76ce4b5b 1422 {
1423
1424 nsigmaTOFPi = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kPion);
1425 nsigmaTOFK = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kKaon);
1426 nsigmaTOFP = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kProton);
2e04885f 1427 nsigmaTOFE = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kElectron);
76ce4b5b 1428
1429 Double_t len=200;// esdtrack->GetIntegratedLength(); !!!!!
1430 Double_t tof=tAodTrack->GetTOFsignal();
1431 if(tof > 0.) vp=len/tof/0.03;
1432 }
87ede832 1433 }
1434 tFemtoTrack->SetVTOF(vp);
1435 tFemtoTrack->SetNSigmaTOFPi(nsigmaTOFPi);
1436 tFemtoTrack->SetNSigmaTOFK(nsigmaTOFK);
1437 tFemtoTrack->SetNSigmaTOFP(nsigmaTOFP);
2e04885f 1438 tFemtoTrack->SetNSigmaTOFE(nsigmaTOFE);
87ede832 1439
76ce4b5b 1440
87ede832 1441 //////////////////////////////////////
76ce4b5b 1442
1443}
1444
1445void AliFemtoEventReaderAOD::SetCentralityPreSelection(double min, double max)
1446{
1447 fCentRange[0] = min; fCentRange[1] = max;
87ede832 1448 fUsePreCent = 1;
1777446b 1449 fEstEventMult = kCentrality;
76ce4b5b 1450}
1451
973a91f8 1452void AliFemtoEventReaderAOD::SetNoCentrality(bool anocent)
1453{
1777446b 1454 if(anocent==false) {fEstEventMult=kCentrality;}
1455 else {fEstEventMult=kReference; fUsePreCent = 0; }
973a91f8 1456}
76ce4b5b 1457
1458void AliFemtoEventReaderAOD::SetAODpidUtil(AliAODpidUtil *aAODpidUtil)
1459{
1460 fAODpidUtil = aAODpidUtil;
1461 // printf("fAODpidUtil: %x\n",fAODpidUtil);
1462}
1463
1777446b 1464void AliFemtoEventReaderAOD::SetAODheader(AliAODHeader *aAODheader)
1465{
1466 fAODheader = aAODheader;
1467}
1468
76ce4b5b 1469
ba3c23a4 1470void AliFemtoEventReaderAOD::SetMagneticFieldSign(int s)
1471{
1472 if(s>0)
1473 fMagFieldSign = 1;
1474 else if(s<0)
1475 fMagFieldSign = -1;
1476 else
1477 fMagFieldSign = 0;
1478}
1479
5e2038a9 1480void AliFemtoEventReaderAOD::SetEPVZERO(Bool_t iepvz)
1481{
87ede832 1482 fisEPVZ = iepvz;
5e2038a9 1483}
76ce4b5b 1484
ce7b3d98 1485void AliFemtoEventReaderAOD::GetGlobalPositionAtGlobalRadiiThroughTPC(AliAODTrack *track, Float_t bfield, Float_t globalPositionsAtRadii[9][3])
1486{
1487 // Gets the global position of the track at nine different radii in the TPC
1488 // track is the track you want to propagate
1489 // bfield is the magnetic field of your event
1490 // globalPositionsAtRadii is the array of global positions in the radii and xyz
1491
1492 // Initialize the array to something indicating there was no propagation
1493 for(Int_t i=0;i<9;i++){
1494 for(Int_t j=0;j<3;j++){
1495 globalPositionsAtRadii[i][j]=-9999.;
1496 }
1497 }
1498
1499 // Make a copy of the track to not change parameters of the track
1500 AliExternalTrackParam etp; etp.CopyFromVTrack(track);
1501 //printf("\nAfter CopyFromVTrack\n");
1502 //etp.Print();
1503
1504 // The global position of the the track
1505 Double_t xyz[3]={-9999.,-9999.,-9999.};
1506
1507 // Counter for which radius we want
1508 Int_t iR=0;
1509 // The radii at which we get the global positions
1510 // IROC (OROC) from 84.1 cm to 132.1 cm (134.6 cm to 246.6 cm)
1511 Float_t Rwanted[9]={85.,105.,125.,145.,165.,185.,205.,225.,245.};
1512 // The global radius we are at
1513 Float_t globalRadius=0;
1514
1515 // Propagation is done in local x of the track
1516 for (Float_t x = etp.GetX();x<247.;x+=1.){ // GetX returns local coordinates
1517 // Starts at the tracks fX and goes outwards. x = 245 is the outer radial limit
1518 // of the TPC when the track is straight, i.e. has inifinite pt and doesn't get bent.
1519 // If the track's momentum is smaller than infinite, it will develop a y-component, which
1520 // adds to the global radius
1521
1522 // Stop if the propagation was not succesful. This can happen for low pt tracks
1523 // that don't reach outer radii
1524 if(!etp.PropagateTo(x,bfield))break;
1525 etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1526 globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
1527
1528 // Roughly reached the radius we want
1529 if(globalRadius > Rwanted[iR]){
1530
1531 // Bigger loop has bad precision, we're nearly one centimeter too far, go back in small steps.
1532 while (globalRadius>Rwanted[iR]){
87ede832 1533 x-=.1;
1534 // printf("propagating to x %5.2f\n",x);
1535 if(!etp.PropagateTo(x,bfield))break;
1536 etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1537 globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
ce7b3d98 1538 }
1539 //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]);
1540 globalPositionsAtRadii[iR][0]=xyz[0];
1541 globalPositionsAtRadii[iR][1]=xyz[1];
1542 globalPositionsAtRadii[iR][2]=xyz[2];
1543 // Indicate we want the next radius
1544 iR+=1;
1545 }
1546 if(iR>=8){
1547 // TPC edge reached
1548 return;
1549 }
1550 }
1551}
76ce4b5b 1552
734c121a 1553void AliFemtoEventReaderAOD::SetpA2013(Bool_t pa2013)
1554{
1555 fpA2013 = pa2013;
1556}
50bb1be0 1557
7dec7414 1558void AliFemtoEventReaderAOD::SetUseMVPlpSelection(Bool_t mvplp)
1559{
1560 fMVPlp = mvplp;
1561}
1562
1563void AliFemtoEventReaderAOD::SetIsPileUpEvent(Bool_t ispileup)
1564{
1565 fisPileUp = ispileup;
1566}
1567
1568
1569
50bb1be0 1570void AliFemtoEventReaderAOD::SetDCAglobalTrack(Bool_t dcagt)
1571{
1572 fDCAglobalTrack = dcagt;
1573}
1574
316b08d6
MS		 1575
1576bool AliFemtoEventReaderAOD::RejectEventCentFlat(float MagField, float CentPercent)
1577{ // to flatten centrality distribution
1578 bool RejectEvent = kFALSE;
1579 int weightBinSign;
1580 TRandom3* fRandomNumber = new TRandom3(); //for 3D, random sign switching
1581 fRandomNumber->SetSeed(0);
1582
1583 if(MagField > 0) weightBinSign = 0;
1584 else weightBinSign = 1;
1585 float kCentWeight[2][9] = {{.878,.876,.860,.859,.859,.88,.873,.879,.894},
1586 {.828,.793,.776,.772,.775,.796,.788,.804,.839}};
1587 int weightBinCent = (int) CentPercent;
1588 if(fRandomNumber->Rndm() > kCentWeight[weightBinSign][weightBinCent]) RejectEvent = kTRUE;
24135013 1589 delete fRandomNumber;
316b08d6
MS		 1590 return RejectEvent;
1591}
1592
1593void AliFemtoEventReaderAOD::SetCentralityFlattening(Bool_t dcagt)
1594{
1595 fFlatCent = dcagt;
1596}