]> 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
Merge branch 'feature-movesplit'
[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
53f88a61 1213 if( ((tFemtoV0->StatusPos() & AliVTrack::kTOFout) == AliVTrack::kTOFout) && ((tFemtoV0->StatusPos() & AliVTrack::kTIME) == AliVTrack::kTIME) ){
1214 // if (tFemtoV0->StatusPos() & AliESDtrack::kTOFout & AliESDtrack::kTIME) { //AliESDtrack::kTOFpid=0x8000
b6200829 1215 probMisPos = fAODpidUtil->GetTOFMismatchProbability(trackpos);
1216 }
53f88a61 1217 if( ((tFemtoV0->StatusNeg() & AliVTrack::kTOFout) == AliVTrack::kTOFout) && ((tFemtoV0->StatusNeg() & AliVTrack::kTIME) == AliVTrack::kTIME) ){
1218 // if (tFemtoV0->StatusNeg() & AliESDtrack::kTOFout & AliESDtrack::kTIME) { //AliESDtrack::kTOFpid=0x8000
b6200829 1219 probMisNeg = fAODpidUtil->GetTOFMismatchProbability(trackneg);
1220 }
1221
53f88a61 1222 // if(// (tFemtoV0->StatusPos()& AliESDtrack::kTOFpid)==0 ||
1223 // (tFemtoV0->StatusPos()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTOFout)==0 || probMisPos > 0.01)
1224
1225 if( !(((tFemtoV0->StatusPos() & AliVTrack::kTOFout) == AliVTrack::kTOFout) && ((tFemtoV0->StatusPos() & AliVTrack::kTIME) == AliVTrack::kTIME)) || probMisPos > 0.01)
1226
87ede832 1227 {
53f88a61 1228 // if(// (tFemtoV0->StatusNeg()&AliESDtrack::kTOFpid)==0 ||
1229 // (tFemtoV0->StatusNeg()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTOFout)==0 || probMisNeg > 0.01)
1230 if( !(((tFemtoV0->StatusNeg() & AliVTrack::kTOFout) == AliVTrack::kTOFout) && ((tFemtoV0->StatusNeg() & AliVTrack::kTIME) == AliVTrack::kTIME)) || probMisNeg > 0.01)
1231
973a91f8 1232 {
ce7b3d98 1233 tFemtoV0->SetPosNSigmaTOFK(-1000);
1234 tFemtoV0->SetNegNSigmaTOFK(-1000);
1235 tFemtoV0->SetPosNSigmaTOFP(-1000);
1236 tFemtoV0->SetNegNSigmaTOFP(-1000);
1237 tFemtoV0->SetPosNSigmaTOFPi(-1000);
1238 tFemtoV0->SetNegNSigmaTOFPi(-1000);
1239
1240 tFemtoV0->SetTOFProtonTimePos(-1000);
1241 tFemtoV0->SetTOFPionTimePos(-1000);
1242 tFemtoV0->SetTOFKaonTimePos(-1000);
1243 tFemtoV0->SetTOFProtonTimeNeg(-1000);
1244 tFemtoV0->SetTOFPionTimeNeg(-1000);
1245 tFemtoV0->SetTOFKaonTimeNeg(-1000);
973a91f8 1246 }
973a91f8 1247 }
87ede832 1248 else
973a91f8 1249 {
53f88a61 1250 if( ((tFemtoV0->StatusPos() & AliVTrack::kTOFout) == AliVTrack::kTOFout) && ((tFemtoV0->StatusPos() & AliVTrack::kTIME) == AliVTrack::kTIME) && probMisPos < 0.01){
1251
1252 // if(trackpos->IsOn(AliESDtrack::kTOFout & AliESDtrack::kTIME)) {
b6200829 1253 tFemtoV0->SetPosNSigmaTOFK(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kKaon));
1254 tFemtoV0->SetPosNSigmaTOFP(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kProton));
1255 tFemtoV0->SetPosNSigmaTOFPi(fAODpidUtil->NumberOfSigmasTOF(trackpos,AliPID::kPion));
1256 }
53f88a61 1257 if( ((tFemtoV0->StatusNeg() & AliVTrack::kTOFout) == AliVTrack::kTOFout) && ((tFemtoV0->StatusNeg() & AliVTrack::kTIME) == AliVTrack::kTIME) && probMisNeg < 0.01){
1258
1259 // if(trackneg->IsOn(AliESDtrack::kTOFout & AliESDtrack::kTIME)) {
b6200829 1260 tFemtoV0->SetNegNSigmaTOFK(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kKaon));
1261 tFemtoV0->SetNegNSigmaTOFP(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kProton));
1262 tFemtoV0->SetNegNSigmaTOFPi(fAODpidUtil->NumberOfSigmasTOF(trackneg,AliPID::kPion));
1263 }
87ede832 1264 double TOFSignalPos = trackpos->GetTOFsignal();
1265 double TOFSignalNeg = trackneg->GetTOFsignal();
b6200829 1266 TOFSignalPos -= fAODpidUtil->GetTOFResponse().GetStartTime(trackpos->P());
1267 TOFSignalNeg -= fAODpidUtil->GetTOFResponse().GetStartTime(trackneg->P());
87ede832 1268 double pidPos[5];
1269 double pidNeg[5];
1270 trackpos->GetIntegratedTimes(pidPos);
1271 trackneg->GetIntegratedTimes(pidNeg);
1272
1273 tFemtoV0->SetTOFPionTimePos(TOFSignalPos-pidPos[2]);
1274 tFemtoV0->SetTOFKaonTimePos(TOFSignalPos-pidPos[3]);
1275 tFemtoV0->SetTOFProtonTimePos(TOFSignalPos-pidPos[4]);
1276 tFemtoV0->SetTOFPionTimeNeg(TOFSignalNeg-pidNeg[2]);
1277 tFemtoV0->SetTOFKaonTimeNeg(TOFSignalNeg-pidNeg[3]);
1278 tFemtoV0->SetTOFProtonTimeNeg(TOFSignalNeg-pidNeg[4]);
973a91f8 1279 }
cbc1c949 1280
87ede832 1281 }
1282 else
1283 {
cbc1c949 1284
87ede832 1285 tFemtoV0->SetStatusPos(999);
1286 tFemtoV0->SetStatusNeg(999);
1287 }
cbc1c949 1288
973a91f8 1289 tFemtoV0->SetOnFlyStatusV0(tAODv0->GetOnFlyStatus());
cbc1c949 1290 return tFemtoV0;
973a91f8 1291}
1292
76ce4b5b 1293void AliFemtoEventReaderAOD::SetFilterBit(UInt_t ibit)
1294{
1295 fFilterBit = (1 << (ibit));
1296}
1297
64536eaf 1298
1299void AliFemtoEventReaderAOD::SetFilterMask(int ibit)
1300{
1301 fFilterMask = ibit;
1302}
1303
76ce4b5b 1304void AliFemtoEventReaderAOD::SetReadMC(unsigned char a)
1305{
1306 fReadMC = a;
1307}
1308
973a91f8 1309
1310void AliFemtoEventReaderAOD::SetReadV0(unsigned char a)
1311{
1312 fReadV0 = a;
1313}
1314
1777446b 1315void AliFemtoEventReaderAOD::SetUseMultiplicity(EstEventMult aType)
1316{
1317 fEstEventMult = aType;
1318}
1319
76ce4b5b 1320AliAODMCParticle* AliFemtoEventReaderAOD::GetParticleWithLabel(TClonesArray *mcP, Int_t aLabel)
1321{
1322 if (aLabel < 0) return 0;
1323 AliAODMCParticle *aodP;
1324 Int_t posstack = 0;
1325 if (aLabel > mcP->GetEntries())
1326 posstack = mcP->GetEntries();
1327 else
1328 posstack = aLabel;
1329
1330 aodP = (AliAODMCParticle *) mcP->At(posstack);
1331 if (aodP->GetLabel() > posstack) {
1332 do {
1333 aodP = (AliAODMCParticle *) mcP->At(posstack);
1334 if (aodP->GetLabel() == aLabel) return aodP;
1335 posstack--;
1336 }
1337 while (posstack > 0);
1338 }
1339 else {
1340 do {
1341 aodP = (AliAODMCParticle *) mcP->At(posstack);
1342 if (aodP->GetLabel() == aLabel) return aodP;
1343 posstack++;
1344 }
1345 while (posstack < mcP->GetEntries());
1346 }
87ede832 1347
76ce4b5b 1348 return 0;
1349}
1350
87ede832 1351void AliFemtoEventReaderAOD::CopyPIDtoFemtoTrack(AliAODTrack *tAodTrack,
1352 AliFemtoTrack *tFemtoTrack)
76ce4b5b 1353{
41cb6c1e 1354
50bb1be0 1355 if (fDCAglobalTrack) {
41cb6c1e 1356
906b97d0 1357 // code from Michael and Prabhat from AliAnalysisTaskDptDptCorrelations
1358 const AliAODVertex* vertex = (AliAODVertex*) fEvent->GetPrimaryVertex();
1359 float vertexX = -999.;
1360 float vertexY = -999.;
1361 float vertexZ = -999.;
1362
1363 if(vertex) {
1364 Double32_t fCov[6];
1365 vertex->GetCovarianceMatrix(fCov);
1366 if(vertex->GetNContributors() > 0) {
1367 if(fCov[5] != 0) {
1368 vertexX = vertex->GetX();
1369 vertexY = vertex->GetY();
1370 vertexZ = vertex->GetZ();
1371
1372 }
1373 }
1374 }
1375
1376 Double_t pos[3];
1377 tAodTrack->GetXYZ(pos);
6691ea4f 1378
906b97d0 1379 Double_t DCAX = pos[0] - vertexX;
1380 Double_t DCAY = pos[1] - vertexY;
1381 Double_t DCAZ = pos[2] - vertexZ;
e42e0efc 1382
906b97d0 1383 Double_t DCAXY = TMath::Sqrt((DCAX*DCAX) + (DCAY*DCAY));
6691ea4f 1384
906b97d0 1385 tFemtoTrack->SetImpactD(DCAXY);
1386 tFemtoTrack->SetImpactZ(DCAZ);
1387
1388
1389 // // // copying DCA information (taking it from global tracks gives better resolution than from TPC-only)
1390
1391 // double impact[2];
1392 // double covimpact[3];
1393
1394 // AliAODTrack* trk_clone = (AliAODTrack*)tAodTrack->Clone("trk_clone");
1395 // // if(!trk_clone->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
1396 // // if (!tAodTrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
1397 // if (!trk_clone->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
1398
1399 // // if (!tAodTrack->PropagateToDCA(fEvent->GetPrimaryVertex(),fEvent->GetMagneticField(),10000,impact,covimpact)) {
1400 // //cout << "sth went wrong with dca propagation" << endl;
1401 // tFemtoTrack->SetImpactD(-1000.0);
1402 // tFemtoTrack->SetImpactZ(-1000.0);
1403
1404 // }
1405 // else {
1406 // tFemtoTrack->SetImpactD(impact[0]);
1407 // tFemtoTrack->SetImpactZ(impact[1]);
1408 // }
1409 // delete trk_clone;
50bb1be0 1410 }
41cb6c1e 1411
76ce4b5b 1412 double aodpid[10];
1413 tAodTrack->GetPID(aodpid);
1414 tFemtoTrack->SetPidProbElectron(aodpid[0]);
1415 tFemtoTrack->SetPidProbMuon(aodpid[1]);
1416 tFemtoTrack->SetPidProbPion(aodpid[2]);
1417 tFemtoTrack->SetPidProbKaon(aodpid[3]);
1418 tFemtoTrack->SetPidProbProton(aodpid[4]);
1419
1420 aodpid[0] = -100000.0;
1421 aodpid[1] = -100000.0;
1422 aodpid[2] = -100000.0;
1423 aodpid[3] = -100000.0;
1424 aodpid[4] = -100000.0;
87ede832 1425
76ce4b5b 1426 double tTOF = 0.0;
87ede832 1427 Float_t probMis = 1.0;
76ce4b5b 1428
4eac0b05 1429 //what is that code? for what do we need it? nsigma values are not enough?
53f88a61 1430 // if (tAodTrack->GetStatus() & AliESDtrack::kTOFout & AliESDtrack::kTIME) { //AliESDtrack::kTOFpid=0x8000
1431
1432 ULong_t status=tAodTrack->GetStatus();
1433
1434 if( ((status & AliVTrack::kTOFout) == AliVTrack::kTOFout) && ((status & AliVTrack::kTIME) == AliVTrack::kTIME) ){
87ede832 1435 tTOF = tAodTrack->GetTOFsignal();
1436 tAodTrack->GetIntegratedTimes(aodpid);
1777446b 1437
87ede832 1438 tTOF -= fAODpidUtil->GetTOFResponse().GetStartTime(tAodTrack->P());
1439
1440 probMis = fAODpidUtil->GetTOFMismatchProbability(tAodTrack);
1441 }
1442
1443
1444
1445 tFemtoTrack->SetTofExpectedTimes(tTOF-aodpid[2], tTOF-aodpid[3], tTOF-aodpid[4]);
76ce4b5b 1446
76ce4b5b 1447 ////// TPC ////////////////////////////////////////////
1448
87ede832 1449 float nsigmaTPCK=-1000.;
1450 float nsigmaTPCPi=-1000.;
1451 float nsigmaTPCP=-1000.;
2e04885f 1452 float nsigmaTPCE=-1000.;
87ede832 1453
76ce4b5b 1454 // cout<<"in reader fESDpid"<<fESDpid<<endl;
1455
53f88a61 1456 nsigmaTPCK = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kKaon);
1457 nsigmaTPCPi = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kPion);
1458 nsigmaTPCP = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kProton);
1459 nsigmaTPCE = fAODpidUtil->NumberOfSigmasTPC(tAodTrack,AliPID::kElectron);
76ce4b5b 1460
1461 tFemtoTrack->SetNSigmaTPCPi(nsigmaTPCPi);
1462 tFemtoTrack->SetNSigmaTPCK(nsigmaTPCK);
1463 tFemtoTrack->SetNSigmaTPCP(nsigmaTPCP);
2e04885f 1464 tFemtoTrack->SetNSigmaTPCE(nsigmaTPCE);
76ce4b5b 1465
87ede832 1466 tFemtoTrack->SetTPCchi2(tAodTrack->Chi2perNDF());
1467 tFemtoTrack->SetTPCncls(tAodTrack->GetTPCNcls());
1468 tFemtoTrack->SetTPCnclsF(tAodTrack->GetTPCNcls());
1469
1470 tFemtoTrack->SetTPCsignalN(1);
1471 tFemtoTrack->SetTPCsignalS(1);
76ce4b5b 1472 tFemtoTrack->SetTPCsignal(tAodTrack->GetTPCsignal());
76ce4b5b 1473
87ede832 1474 ///////TOF//////////////////////
76ce4b5b 1475
87ede832 1476 float vp=-1000.;
1477 float nsigmaTOFPi=-1000.;
1478 float nsigmaTOFK=-1000.;
1479 float nsigmaTOFP=-1000.;
2e04885f 1480 float nsigmaTOFE=-1000.;
87ede832 1481
53f88a61 1482 if( ((status & AliVTrack::kTOFout) == AliVTrack::kTOFout) && ((status & AliVTrack::kTIME) == AliVTrack::kTIME) && probMis < 0.01){
76ce4b5b 1483
53f88a61 1484 nsigmaTOFPi = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kPion);
1485 nsigmaTOFK = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kKaon);
1486 nsigmaTOFP = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kProton);
1487 nsigmaTOFE = fAODpidUtil->NumberOfSigmasTOF(tAodTrack,AliPID::kElectron);
76ce4b5b 1488
53f88a61 1489 Double_t len=200;// esdtrack->GetIntegratedLength(); !!!!!
1490 Double_t tof=tAodTrack->GetTOFsignal();
1491 if(tof > 0.) vp=len/tof/0.03;
87ede832 1492 }
53f88a61 1493
87ede832 1494 tFemtoTrack->SetVTOF(vp);
1495 tFemtoTrack->SetNSigmaTOFPi(nsigmaTOFPi);
1496 tFemtoTrack->SetNSigmaTOFK(nsigmaTOFK);
1497 tFemtoTrack->SetNSigmaTOFP(nsigmaTOFP);
2e04885f 1498 tFemtoTrack->SetNSigmaTOFE(nsigmaTOFE);
87ede832 1499
76ce4b5b 1500
87ede832 1501 //////////////////////////////////////
76ce4b5b 1502
1503}
1504
1505void AliFemtoEventReaderAOD::SetCentralityPreSelection(double min, double max)
1506{
1507 fCentRange[0] = min; fCentRange[1] = max;
87ede832 1508 fUsePreCent = 1;
1777446b 1509 fEstEventMult = kCentrality;
76ce4b5b 1510}
1511
973a91f8 1512void AliFemtoEventReaderAOD::SetNoCentrality(bool anocent)
1513{
1777446b 1514 if(anocent==false) {fEstEventMult=kCentrality;}
1515 else {fEstEventMult=kReference; fUsePreCent = 0; }
973a91f8 1516}
76ce4b5b 1517
1518void AliFemtoEventReaderAOD::SetAODpidUtil(AliAODpidUtil *aAODpidUtil)
1519{
1520 fAODpidUtil = aAODpidUtil;
1521 // printf("fAODpidUtil: %x\n",fAODpidUtil);
1522}
1523
1777446b 1524void AliFemtoEventReaderAOD::SetAODheader(AliAODHeader *aAODheader)
1525{
1526 fAODheader = aAODheader;
1527}
1528
76ce4b5b 1529
ba3c23a4 1530void AliFemtoEventReaderAOD::SetMagneticFieldSign(int s)
1531{
1532 if(s>0)
1533 fMagFieldSign = 1;
1534 else if(s<0)
1535 fMagFieldSign = -1;
1536 else
1537 fMagFieldSign = 0;
1538}
1539
5e2038a9 1540void AliFemtoEventReaderAOD::SetEPVZERO(Bool_t iepvz)
1541{
87ede832 1542 fisEPVZ = iepvz;
5e2038a9 1543}
76ce4b5b 1544
ce7b3d98 1545void AliFemtoEventReaderAOD::GetGlobalPositionAtGlobalRadiiThroughTPC(AliAODTrack *track, Float_t bfield, Float_t globalPositionsAtRadii[9][3])
1546{
1547 // Gets the global position of the track at nine different radii in the TPC
1548 // track is the track you want to propagate
1549 // bfield is the magnetic field of your event
1550 // globalPositionsAtRadii is the array of global positions in the radii and xyz
1551
1552 // Initialize the array to something indicating there was no propagation
1553 for(Int_t i=0;i<9;i++){
1554 for(Int_t j=0;j<3;j++){
1555 globalPositionsAtRadii[i][j]=-9999.;
1556 }
1557 }
1558
1559 // Make a copy of the track to not change parameters of the track
1560 AliExternalTrackParam etp; etp.CopyFromVTrack(track);
1561 //printf("\nAfter CopyFromVTrack\n");
1562 //etp.Print();
1563
1564 // The global position of the the track
1565 Double_t xyz[3]={-9999.,-9999.,-9999.};
1566
1567 // Counter for which radius we want
1568 Int_t iR=0;
1569 // The radii at which we get the global positions
1570 // IROC (OROC) from 84.1 cm to 132.1 cm (134.6 cm to 246.6 cm)
1571 Float_t Rwanted[9]={85.,105.,125.,145.,165.,185.,205.,225.,245.};
1572 // The global radius we are at
1573 Float_t globalRadius=0;
1574
1575 // Propagation is done in local x of the track
1576 for (Float_t x = etp.GetX();x<247.;x+=1.){ // GetX returns local coordinates
1577 // Starts at the tracks fX and goes outwards. x = 245 is the outer radial limit
1578 // of the TPC when the track is straight, i.e. has inifinite pt and doesn't get bent.
1579 // If the track's momentum is smaller than infinite, it will develop a y-component, which
1580 // adds to the global radius
1581
1582 // Stop if the propagation was not succesful. This can happen for low pt tracks
1583 // that don't reach outer radii
1584 if(!etp.PropagateTo(x,bfield))break;
1585 etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1586 globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
1587
1588 // Roughly reached the radius we want
1589 if(globalRadius > Rwanted[iR]){
1590
1591 // Bigger loop has bad precision, we're nearly one centimeter too far, go back in small steps.
1592 while (globalRadius>Rwanted[iR]){
87ede832 1593 x-=.1;
1594 // printf("propagating to x %5.2f\n",x);
1595 if(!etp.PropagateTo(x,bfield))break;
1596 etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1597 globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
ce7b3d98 1598 }
1599 //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]);
1600 globalPositionsAtRadii[iR][0]=xyz[0];
1601 globalPositionsAtRadii[iR][1]=xyz[1];
1602 globalPositionsAtRadii[iR][2]=xyz[2];
1603 // Indicate we want the next radius
1604 iR+=1;
1605 }
1606 if(iR>=8){
1607 // TPC edge reached
1608 return;
1609 }
1610 }
1611}
76ce4b5b 1612
734c121a 1613void AliFemtoEventReaderAOD::SetpA2013(Bool_t pa2013)
1614{
1615 fpA2013 = pa2013;
1616}
50bb1be0 1617
7dec7414 1618void AliFemtoEventReaderAOD::SetUseMVPlpSelection(Bool_t mvplp)
1619{
1620 fMVPlp = mvplp;
1621}
1622
1623void AliFemtoEventReaderAOD::SetIsPileUpEvent(Bool_t ispileup)
1624{
1625 fisPileUp = ispileup;
1626}
1627
1628
1629
50bb1be0 1630void AliFemtoEventReaderAOD::SetDCAglobalTrack(Bool_t dcagt)
1631{
1632 fDCAglobalTrack = dcagt;
1633}
1634
316b08d6
MS		 1635
1636bool AliFemtoEventReaderAOD::RejectEventCentFlat(float MagField, float CentPercent)
1637{ // to flatten centrality distribution
1638 bool RejectEvent = kFALSE;
1639 int weightBinSign;
1640 TRandom3* fRandomNumber = new TRandom3(); //for 3D, random sign switching
1641 fRandomNumber->SetSeed(0);
1642
1643 if(MagField > 0) weightBinSign = 0;
1644 else weightBinSign = 1;
1645 float kCentWeight[2][9] = {{.878,.876,.860,.859,.859,.88,.873,.879,.894},
1646 {.828,.793,.776,.772,.775,.796,.788,.804,.839}};
1647 int weightBinCent = (int) CentPercent;
1648 if(fRandomNumber->Rndm() > kCentWeight[weightBinSign][weightBinCent]) RejectEvent = kTRUE;
24135013 1649 delete fRandomNumber;
316b08d6
MS		 1650 return RejectEvent;
1651}
1652
1653void AliFemtoEventReaderAOD::SetCentralityFlattening(Bool_t dcagt)
1654{
1655 fFlatCent = dcagt;
1656}