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