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