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