]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/FEMTOSCOPY/AliFemto/AliFemtoEventReaderESDChain.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
replacing kTOFpid by kTOFout && kTIME as requested by TOF PID group
[u/mrichter/AliRoot.git] / PWGCF / FEMTOSCOPY / AliFemto / AliFemtoEventReaderESDChain.cxx
CommitLineData
76ce4b5b 1////////////////////////////////////////////////////////////////////////////////
2/// ///
3/// AliFemtoEventReaderESDChain - the reader class for the Alice ESD ///
4/// tailored for the Task framework ///
5/// Reads in AliESDfriend to create shared hit/quality information ///
6/// Authors: Adam Kisiel kisiel@mps.ohio-state.edu ///
7/// ///
8////////////////////////////////////////////////////////////////////////////////
9#include "AliFemtoEventReaderESDChain.h"
10
11#include "TFile.h"
12#include "TTree.h"
13#include "AliESDEvent.h"
14#include "AliESDtrack.h"
15#include "AliESDVertex.h"
16#include "AliMultiplicity.h"
17#include "AliCentrality.h"
18#include "AliEventplane.h"
19#include "AliESDVZERO.h"
20#include "AliFmPhysicalHelixD.h"
21#include "AliFmThreeVectorF.h"
22#include "SystemOfUnits.h"
23#include "AliFemtoEvent.h"
24#include "AliFemtoModelHiddenInfo.h"
25#include "AliPID.h"
7dec7414 26#include "AliAnalysisUtils.h"
76ce4b5b 27
28ClassImp(AliFemtoEventReaderESDChain)
29
30#if !(ST_NO_NAMESPACES)
31 using namespace units;
32#endif
33
34using namespace std;
973a91f8 35
76ce4b5b 36//____________________________
37AliFemtoEventReaderESDChain::AliFemtoEventReaderESDChain():
38 fFileName(" "),
39 fConstrained(true),
40 fReadInner(false),
41 fUseTPCOnly(false),
42 fNumberofEvent(0),
43 fCurEvent(0),
44 fCurFile(0),
45 fEvent(0x0),
76ce4b5b 46 fTrackType(kGlobal),
4eac0b05 47 fEstEventMult(kReferenceITSTPC),
76ce4b5b 48 fEventTrig(AliVEvent::kMB), //trigger
49 fESDpid(0),
973a91f8 50 fIsPidOwner(0),
ce7b3d98 51 fReadV0(0),
734c121a 52 fMagFieldSign(0),
7dec7414 53 fpA2013(kFALSE),
54 fisPileUp(kFALSE),
7d3e2025 55 fMVPlp(kFALSE),
56 fMinVtxContr(0),
57 fMinPlpContribMV(0),
58 fMinPlpContribSPD(0)
76ce4b5b 59{
60 //constructor with 0 parameters , look at default settings
61 // fClusterPerPadrow = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
62 // for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
63 // fClusterPerPadrow[tPad] = new list<Int_t>();
64 // }
65 // fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
66 // for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
67 // fSharedList[tPad] = new list<Int_t>();
68 // }
69}
70
71//__________________
72AliFemtoEventReaderESDChain::AliFemtoEventReaderESDChain(const AliFemtoEventReaderESDChain& aReader):
73 AliFemtoEventReader(aReader),
74 fFileName(" "),
75 fConstrained(true),
76 fReadInner(false),
77 fUseTPCOnly(false),
78 fNumberofEvent(0),
79 fCurEvent(0),
80 fCurFile(0),
81 fEvent(0x0),
76ce4b5b 82 fTrackType(kGlobal),
4eac0b05 83 fEstEventMult(kReferenceITSTPC),
76ce4b5b 84 fEventTrig(AliVEvent::kMB), //trigger
85 fESDpid(0),
973a91f8 86 fIsPidOwner(0),
ce7b3d98 87 fReadV0(0),
734c121a 88 fMagFieldSign(0),
7dec7414 89 fpA2013(kFALSE),
90 fisPileUp(kFALSE),
7d3e2025 91 fMVPlp(kFALSE),
92 fMinVtxContr(0),
93 fMinPlpContribMV(0),
94 fMinPlpContribSPD(0)
76ce4b5b 95{
96 // Copy constructor
97 fConstrained = aReader.fConstrained;
98 fReadInner = aReader.fReadInner;
99 fUseTPCOnly = aReader.fUseTPCOnly;
100 fNumberofEvent = aReader.fNumberofEvent;
101 fCurEvent = aReader.fCurEvent;
102 fCurFile = aReader.fCurFile;
103 // fEvent = new AliESD(*aReader.fEvent);
104 fEvent = new AliESDEvent();
76ce4b5b 105 fTrackType = aReader.fTrackType;
106 fEstEventMult = aReader.fEstEventMult;
107 fEventTrig = aReader.fEventTrig; //trigger
ce7b3d98 108 fReadV0 = aReader.fReadV0;
109 fMagFieldSign = aReader.fMagFieldSign;
734c121a 110 fpA2013 = aReader.fpA2013;
7dec7414 111 fisPileUp = aReader.fisPileUp;
112 fMVPlp = aReader.fMVPlp;
7d3e2025 113 fMinVtxContr = aReader.fMinVtxContr;
114 fMinPlpContribMV = aReader.fMinPlpContribMV;
115 fMinPlpContribSPD = aReader.fMinPlpContribSPD;
116
76ce4b5b 117 // fEventFriend = aReader.fEventFriend;
118 // fClusterPerPadrow = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
119 // for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
120 // fClusterPerPadrow[tPad] = new list<Int_t>();
121 // list<Int_t>::iterator iter;
122 // for (iter=aReader.fClusterPerPadrow[tPad]->begin(); iter!=aReader.fClusterPerPadrow[tPad]->end(); iter++) {
123 // fClusterPerPadrow[tPad]->push_back(*iter);
124 // }
125 // }
126 // fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
127 // for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
128 // fSharedList[tPad] = new list<Int_t>();
129 // list<Int_t>::iterator iter;
130 // for (iter=aReader.fSharedList[tPad]->begin(); iter!=aReader.fSharedList[tPad]->end(); iter++) {
131 // fSharedList[tPad]->push_back(*iter);
132 // }
133 // }
134}
135//__________________
136AliFemtoEventReaderESDChain::~AliFemtoEventReaderESDChain()
137{
138 //Destructor
139 delete fEvent;
140
141 // for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
142 // fClusterPerPadrow[tPad]->clear();
143 // delete fClusterPerPadrow[tPad];
144 // }
145 // delete [] fClusterPerPadrow;
146 // for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
147 // fSharedList[tPad]->clear();
148 // delete fSharedList[tPad];
149 // }
150 // delete [] fSharedList;
4eac0b05 151
76ce4b5b 152}
153
154//__________________
155AliFemtoEventReaderESDChain& AliFemtoEventReaderESDChain::operator=(const AliFemtoEventReaderESDChain& aReader)
156{
157 // Assignment operator
158 if (this == &aReader)
159 return *this;
160
161 fConstrained = aReader.fConstrained;
162 fReadInner = aReader.fReadInner;
163 fUseTPCOnly = aReader.fUseTPCOnly;
164 fNumberofEvent = aReader.fNumberofEvent;
165 fCurEvent = aReader.fCurEvent;
166 fCurFile = aReader.fCurFile;
167 if (fEvent) delete fEvent;
168 fEvent = new AliESDEvent();
169 fTrackType = aReader.fTrackType;
170 fEstEventMult = aReader.fEstEventMult;
171
973a91f8 172 fReadV0 = aReader.fReadV0;
ce7b3d98 173 fMagFieldSign = aReader.fMagFieldSign;
734c121a 174 fpA2013 = aReader.fpA2013;
7dec7414 175 fisPileUp = aReader.fisPileUp;
176 fMVPlp = aReader.fMVPlp;
76ce4b5b 177 // fEventFriend = aReader.fEventFriend;
178
179 // if (fClusterPerPadrow) {
180 // for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
181 // fClusterPerPadrow[tPad]->clear();
182 // delete fClusterPerPadrow[tPad];
183 // }
184 // delete [] fClusterPerPadrow;
185 // }
186
187 // if (fSharedList) {
188 // for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
189 // fSharedList[tPad]->clear();
190 // delete fSharedList[tPad];
191 // }
192 // delete [] fSharedList;
193 // }
194
195 // fClusterPerPadrow = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
196 // for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
197 // fClusterPerPadrow[tPad] = new list<Int_t>();
198 // list<Int_t>::iterator iter;
199 // for (iter=aReader.fClusterPerPadrow[tPad]->begin(); iter!=aReader.fClusterPerPadrow[tPad]->end(); iter++) {
200 // fClusterPerPadrow[tPad]->push_back(*iter);
201 // }
202 // }
203 // fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
204 // for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
205 // fSharedList[tPad] = new list<Int_t>();
206 // list<Int_t>::iterator iter;
207 // for (iter=aReader.fSharedList[tPad]->begin(); iter!=aReader.fSharedList[tPad]->end(); iter++) {
208 // fSharedList[tPad]->push_back(*iter);
209 // }
210 // }
211
212 return *this;
213}
214//__________________
215// Simple report
216AliFemtoString AliFemtoEventReaderESDChain::Report()
217{
218 AliFemtoString temp = "\n This is the AliFemtoEventReaderESDChain\n";
219 return temp;
220}
221
222//__________________
223void AliFemtoEventReaderESDChain::SetConstrained(const bool constrained)
224{
225 // Select whether to read constrained or not constrained momentum
226 fConstrained=constrained;
227}
228//__________________
229bool AliFemtoEventReaderESDChain::GetConstrained() const
230{
231 // Check whether we read constrained or not constrained momentum
232 return fConstrained;
233}
234//__________________
235void AliFemtoEventReaderESDChain::SetReadTPCInner(const bool readinner)
236{
237 fReadInner=readinner;
238}
239
240bool AliFemtoEventReaderESDChain::GetReadTPCInner() const
241{
242 return fReadInner;
243}
244
245//__________________
246void AliFemtoEventReaderESDChain::SetUseTPCOnly(const bool usetpconly)
247{
248 fUseTPCOnly=usetpconly;
249}
250
251bool AliFemtoEventReaderESDChain::GetUseTPCOnly() const
252{
253 return fUseTPCOnly;
254}
255
76ce4b5b 256void AliFemtoEventReaderESDChain::SetUseMultiplicity(EstEventMult aType)
257{
258 fEstEventMult = aType;
259}
260
261AliFemtoEvent* AliFemtoEventReaderESDChain::ReturnHbtEvent()
262{
263 // Get the event, read all the relevant information
264 // and fill the AliFemtoEvent class
265 // Returns a valid AliFemtoEvent
266 AliFemtoEvent *hbtEvent = 0;
734c121a 267
76ce4b5b 268
734c121a 269 hbtEvent = new AliFemtoEvent;
270
271 CopyESDtoFemtoEvent(hbtEvent);
272
273 fCurEvent++;
274
275
276 return hbtEvent;
277}
278
279void AliFemtoEventReaderESDChain::CopyESDtoFemtoEvent(AliFemtoEvent *hbtEvent)
280{
281
282
283 //string tFriendFileName;
76ce4b5b 284 // Get the friend information
285 if (Debug()>1) cout<<"starting to read event "<<fCurEvent<<endl;
286 // fEvent->SetESDfriend(fEventFriend);
287 if(fEvent->GetAliESDOld())fEvent->CopyFromOldESD();
288
734c121a 289
76ce4b5b 290
4eac0b05 291 //setting basic things
76ce4b5b 292 // hbtEvent->SetEventNumber(fEvent->GetEventNumber());
293 hbtEvent->SetRunNumber(fEvent->GetRunNumber());
294 //hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
295 hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
296 hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
297 hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
298 hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
299 hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
300 hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
301 hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
302 hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
303 // hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
304
305 if ((fEvent->IsTriggerClassFired("CINT1WU-B-NOPF-ALL")) ||
306 (fEvent->IsTriggerClassFired("CINT1B-ABCE-NOPF-ALL")) ||
307 (fEvent->IsTriggerClassFired("CINT1-B-NOPF-ALLNOTRD")) ||
308 (fEvent->IsTriggerClassFired("CINT1-B-NOPF-FASTNOTRD")))
309 hbtEvent->SetTriggerCluster(1);
310 else if ((fEvent->IsTriggerClassFired("CSH1WU-B-NOPF-ALL")) ||
311 (fEvent->IsTriggerClassFired("CSH1-B-NOPF-ALLNOTRD")))
312 hbtEvent->SetTriggerCluster(2);
313 else
314 hbtEvent->SetTriggerCluster(0);
315
316 //Vertex
317 double fV1[3];
318 double fVCov[6];
734c121a 319
7dec7414 320 //AliAnalysisUtils
321 if(fisPileUp||fpA2013)
734c121a 322 {
7dec7414 323 AliAnalysisUtils *anaUtil=new AliAnalysisUtils();
7d3e2025 324 if(fMinVtxContr)
325 anaUtil->SetMinVtxContr(fMinVtxContr);
7dec7414 326 if(fpA2013)
327 if(anaUtil->IsVertexSelected2013pA(fEvent)==kFALSE) return; //Vertex rejection for pA analysis.
328 if(fMVPlp) anaUtil->SetUseMVPlpSelection(kTRUE);
329 else anaUtil->SetUseMVPlpSelection(kFALSE);
7d3e2025 330 if(fMinPlpContribMV) anaUtil->SetMinPlpContribMV(fMinPlpContribMV);
331 if(fMinPlpContribSPD) anaUtil->SetMinPlpContribSPD(fMinPlpContribSPD);
7dec7414 332 if(fisPileUp)
333 if(anaUtil->IsPileUpEvent(fEvent)) return; //Pile-up rejection.
334 delete anaUtil;
734c121a 335 }
7dec7414 336
337
76ce4b5b 338 if (fUseTPCOnly) {
7dec7414 339 const AliESDVertex* esdvertex = (AliESDVertex*) fEvent->GetPrimaryVertexTPC();
340 if(!esdvertex || esdvertex->GetNContributors() < 1) return; //Bad vertex, skip event.
341
342 esdvertex->GetXYZ(fV1);
343 esdvertex->GetCovMatrix(fVCov);
344 if (!esdvertex->GetStatus())
76ce4b5b 345 fVCov[4] = -1001.0;
346 }
347 else {
7dec7414 348 const AliESDVertex* esdvertex = (AliESDVertex*) fEvent->GetPrimaryVertex();
349 if(!esdvertex || esdvertex->GetNContributors() < 1) return; //Bad vertex, skip event.
350
351 esdvertex->GetXYZ(fV1);
352 esdvertex->GetCovMatrix(fVCov);
353 if (!esdvertex->GetStatus())
76ce4b5b 354 fVCov[4] = -1001.0;
355 }
734c121a 356
76ce4b5b 357 AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
358 hbtEvent->SetPrimVertPos(vertex);
359 hbtEvent->SetPrimVertCov(fVCov);
360
361 Int_t spdetaonecount = 0;
362
363 for (int iter=0; iter<fEvent->GetMultiplicity()->GetNumberOfTracklets(); iter++)
364 if (fabs(fEvent->GetMultiplicity()->GetEta(iter)) < 1.0)
365 spdetaonecount++;
366
367 // hbtEvent->SetSPDMult(fEvent->GetMultiplicity()->GetNumberOfTracklets());
368 hbtEvent->SetSPDMult(spdetaonecount);
369
370 //starting to reading tracks
371 int nofTracks=0; //number of reconstructed tracks in event
372 nofTracks=fEvent->GetNumberOfTracks();
373 int realnofTracks=0;//number of track which we use ina analysis
374
375 // // Clear the shared cluster list
376 // for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
377 // fClusterPerPadrow[tPad]->clear();
378 // }
379 // for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
380 // fSharedList[tPad]->clear();
381 // }
382
383
384 // for (int i=0;i<nofTracks;i++) {
385 // const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
386
387 // list<Int_t>::iterator tClustIter;
388
389 // Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster];
390 // Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices);
391 // for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
392 // if (tTrackIndices[tNcl] >= 0) {
393 // tClustIter = find(fClusterPerPadrow[tNcl]->begin(), fClusterPerPadrow[tNcl]->end(), tTrackIndices[tNcl]);
394 // if (tClustIter == fClusterPerPadrow[tNcl]->end()) {
395 // fClusterPerPadrow[tNcl]->push_back(tTrackIndices[tNcl]);
396 // }
397 // else {
398 // fSharedList[tNcl]->push_back(tTrackIndices[tNcl]);
399 // }
400 // }
401 // }
402
403 // }
404
405 int tNormMult = 0;
406 int tNormMultPos = 0;
407 int tNormMultNeg = 0;
408
409 Float_t tTotalPt = 0.0;
410
411 Float_t b[2];
412 Float_t bCov[3];
413
1938eadf 414 Int_t tTracklet=0, tITSTPC=0;
76ce4b5b 415
4eac0b05 416 //W-AliESDEvent::EstimateMultiplicity: This obsolete method will be eliminated soon. Use AliESDtrackCuts::GetReferenceMultiplicity
417 //fEvent->EstimateMultiplicity(tTracklet, tITSTPC, tITSPure, 1.2);
76ce4b5b 418
419 hbtEvent->SetMultiplicityEstimateITSTPC(tITSTPC);
420 hbtEvent->SetMultiplicityEstimateTracklets(tTracklet);
421 // hbtEvent->SetMultiplicityEstimateITSPure(tITSPure);
422 hbtEvent->SetMultiplicityEstimateITSPure(fEvent->GetMultiplicity()->GetNumberOfITSClusters(1));
423
424 for (int i=0;i<nofTracks;i++)
425 {
426 bool tGoodMomentum=true; //flaga to chcek if we can read momentum of this track
427
428 const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
429 // const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
430
431 if ((esdtrack->GetStatus() & AliESDtrack::kTPCrefit) &&
432 (esdtrack->GetStatus() & AliESDtrack::kITSrefit)) {
433 if (esdtrack->GetTPCNcls() > 70)
434 if (esdtrack->GetTPCchi2()/esdtrack->GetTPCNcls() < 4.0) {
1777446b 435 if (esdtrack->Pt() > 0.15 && esdtrack->Pt() < 20)
436 if (TMath::Abs(esdtrack->Eta()) < 0.8) {
437 esdtrack->GetImpactParameters(b,bCov);
438 if ((b[0]<0.2) && (b[1] < 2.0)) {
439 tNormMult++;
440 tTotalPt += esdtrack->Pt();
441 }
76ce4b5b 442 }
76ce4b5b 443 }
444 }
1777446b 445
76ce4b5b 446 hbtEvent->SetZDCEMEnergy(tTotalPt);
447 // if (esdtrack->GetStatus() & AliESDtrack::kTPCrefit)
448 // if (esdtrack->GetTPCNcls() > 80)
449 // if (esdtrack->GetTPCchi2()/esdtrack->GetTPCNcls() < 6.0)
450 // if (esdtrack->GetConstrainedParam())
451 // if (fabs(esdtrack->GetConstrainedParam()->Eta()) < 0.5)
452 // if (esdtrack->GetConstrainedParam()->Pt() < 1.0) {
453 // if (esdtrack->GetSign() > 0)
454 // tNormMultPos++;
455 // else if (esdtrack->GetSign() < 0)
456 // tNormMultNeg--;
457 // }
458
459 // If reading ITS-only tracks, reject all with TPC
460 if (fTrackType == kITSOnly) {
461 if (esdtrack->GetStatus() & AliESDtrack::kTPCrefit) continue;
462 if (!(esdtrack->GetStatus() & AliESDtrack::kITSrefit)) continue;
463 if (esdtrack->GetStatus() & AliESDtrack::kTPCin) continue;
464 UChar_t iclm = esdtrack->GetITSClusterMap();
465 Int_t incls = 0;
466 for (int iter=0; iter<6; iter++) if (iclm&(1<<iter)) incls++;
467 if (incls<=3) {
468 if (Debug()>1) cout << "Rejecting track with " << incls << " clusters" << endl;
469 continue;
470 }
471 }
472
473 AliFemtoTrack* trackCopy = new AliFemtoTrack();
474 trackCopy->SetCharge((short)esdtrack->GetSign());
475
476 //in aliroot we have AliPID
477 //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
478 //we use only 5 first
479 double esdpid[5];
480 // esdtrack->GetESDpid(esdpid);
481 esdtrack->GetTPCpid(esdpid);
482 trackCopy->SetPidProbElectron(esdpid[0]);
483 trackCopy->SetPidProbMuon(esdpid[1]);
484 trackCopy->SetPidProbPion(esdpid[2]);
485 trackCopy->SetPidProbKaon(esdpid[3]);
486 trackCopy->SetPidProbProton(esdpid[4]);
2952eb1d 487
76ce4b5b 488 esdpid[0] = -100000.0;
489 esdpid[1] = -100000.0;
490 esdpid[2] = -100000.0;
491 esdpid[3] = -100000.0;
492 esdpid[4] = -100000.0;
2952eb1d 493
76ce4b5b 494 double tTOF = 0.0;
495
2952eb1d 496 if (esdtrack->GetStatus()&AliESDtrack::kTOFout&AliESDtrack::kTIME) {
76ce4b5b 497 tTOF = esdtrack->GetTOFsignal();
498 esdtrack->GetIntegratedTimes(esdpid);
499 }
500
501 trackCopy->SetTofExpectedTimes(tTOF-esdpid[2], tTOF-esdpid[3], tTOF-esdpid[4]);
502
503 ////// TPC ////////////////////////////////////////////
504
2952eb1d 505 float nsigmaTPCK=-1000.;
506 float nsigmaTPCPi=-1000.;
507 float nsigmaTPCP=-1000.;
508
509
76ce4b5b 510 if ((fESDpid) && (esdtrack->IsOn(AliESDtrack::kTPCpid))){
511 nsigmaTPCK = fESDpid->NumberOfSigmasTPC(esdtrack,AliPID::kKaon);
512 nsigmaTPCPi = fESDpid->NumberOfSigmasTPC(esdtrack,AliPID::kPion);
513 nsigmaTPCP = fESDpid->NumberOfSigmasTPC(esdtrack,AliPID::kProton);
2952eb1d 514
76ce4b5b 515 }
516 trackCopy->SetNSigmaTPCPi(nsigmaTPCPi);
517 trackCopy->SetNSigmaTPCK(nsigmaTPCK);
518 trackCopy->SetNSigmaTPCP(nsigmaTPCP);
2952eb1d 519
76ce4b5b 520 ///// TOF ///////////////////////////////////////////////
2952eb1d 521
76ce4b5b 522 float vp=-1000.;
523 float nsigmaTOFPi=-1000.;
524 float nsigmaTOFK=-1000.;
525 float nsigmaTOFP=-1000.;
2952eb1d 526
527 if (// (esdtrack->GetStatus()&AliESDtrack::kTOFpid) &&
76ce4b5b 528 (esdtrack->GetStatus()&AliESDtrack::kTOFout) &&
1938eadf 529 (esdtrack->GetStatus()&AliESDtrack::kTIME))
76ce4b5b 530 {
531
532 //if ((esdtrack->GetStatus()&AliESDtrack::kTOFpid) &&
533 //(esdtrack->GetStatus()&AliESDtrack::kTOFout) &&
534 //(esdtrack->GetStatus()&AliESDtrack::kTIME)){
535 // collect info from ESDpid class
2952eb1d 536
537 if ((fESDpid) && (esdtrack->IsOn(AliESDtrack::kTOFout & AliESDtrack::kTIME))) {
538
539
76ce4b5b 540 double tZero = fESDpid->GetTOFResponse().GetStartTime(esdtrack->P());
2952eb1d 541
76ce4b5b 542 nsigmaTOFPi = fESDpid->NumberOfSigmasTOF(esdtrack,AliPID::kPion,tZero);
543 nsigmaTOFK = fESDpid->NumberOfSigmasTOF(esdtrack,AliPID::kKaon,tZero);
544 nsigmaTOFP = fESDpid->NumberOfSigmasTOF(esdtrack,AliPID::kProton,tZero);
2952eb1d 545
76ce4b5b 546 Double_t len=esdtrack->GetIntegratedLength();
547 Double_t tof=esdtrack->GetTOFsignal();
548 if(tof > 0.) vp=len/tof/0.03;
549 }
550 }
2952eb1d 551
76ce4b5b 552 trackCopy->SetVTOF(vp);
553 trackCopy->SetNSigmaTOFPi(nsigmaTOFPi);
554 trackCopy->SetNSigmaTOFK(nsigmaTOFK);
555 trackCopy->SetNSigmaTOFP(nsigmaTOFP);
2952eb1d 556
76ce4b5b 557 double pxyz[3];
558 double rxyz[3];
559 double impact[2];
560 double covimpact[3];
2952eb1d 561
76ce4b5b 562 if (fUseTPCOnly) {
563 if (!esdtrack->GetTPCInnerParam()) {
564 delete trackCopy;
565 continue;
566 }
567
568
569 AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
570 param->GetXYZ(rxyz);
571 param->PropagateToDCA(fEvent->GetPrimaryVertexTPC(), (fEvent->GetMagneticField()), 10000, impact, covimpact);
572 param->GetPxPyPz(pxyz);//reading noconstarined momentum
573
574 if (fReadInner == true) {
575 AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
576 tInfo->SetPDGPid(211);
577 tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]);
578 tInfo->SetMass(0.13957);
579 // tInfo->SetEmissionPoint(rxyz[0], rxyz[1], rxyz[2], 0.0);
580 // tInfo->SetEmissionPoint(fV1[0], fV1[1], fV1[2], 0.0);
581 tInfo->SetEmissionPoint(rxyz[0]-fV1[0], rxyz[1]-fV1[1], rxyz[2]-fV1[2], 0.0);
582 trackCopy->SetHiddenInfo(tInfo);
583 }
584
585 AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
586 if (v.Mag() < 0.0001) {
587 // cout << "Found 0 momentum ???? " <<endl;
588 delete trackCopy;
589 continue;
590 }
591 trackCopy->SetP(v);//setting momentum
592 trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
593
594 const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
595 const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
596 //setting helix I do not if it is ok
597 AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
598 trackCopy->SetHelix(helix);
599
600 //some stuff which could be useful
601 trackCopy->SetImpactD(impact[0]);
602 trackCopy->SetImpactZ(impact[1]);
603 trackCopy->SetCdd(covimpact[0]);
604 trackCopy->SetCdz(covimpact[1]);
605 trackCopy->SetCzz(covimpact[2]);
606 trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact, covimpact));
607
608 delete param;
609 }
610 else {
611 if (fReadInner == true) {
612
613 if (esdtrack->GetTPCInnerParam()) {
7eb9f5d0 614 AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetInnerParam());
615 //trackCopy->SetInnerMomentum(param->P());
616 trackCopy->SetInnerMomentum(esdtrack->GetTPCmomentum());
76ce4b5b 617 param->GetXYZ(rxyz);
618 // param->PropagateToDCA(fEvent->GetPrimaryVertex(), (fEvent->GetMagneticField()), 10000);
619 param->GetPxPyPz(pxyz);//reading noconstarined momentum
620 delete param;
621
622 AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
623 tInfo->SetPDGPid(211);
624 tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]);
625 tInfo->SetMass(0.13957);
626 // tInfo->SetEmissionPoint(rxyz[0], rxyz[1], rxyz[2], 0.0);
627 //tInfo->SetEmissionPoint(fV1[0], fV1[1], fV1[2], 0.0);
628 tInfo->SetEmissionPoint(rxyz[0]-fV1[0], rxyz[1]-fV1[1], rxyz[2]-fV1[2], 0.0);
629 trackCopy->SetHiddenInfo(tInfo);
630 }
631 }
632
633 if (fTrackType == kGlobal) {
634 if (fConstrained==true)
635 tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
636 else
637 tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
638 }
639 else if (fTrackType == kTPCOnly) {
640 if (esdtrack->GetTPCInnerParam())
641 esdtrack->GetTPCInnerParam()->GetPxPyPz(pxyz);
642 else {
643 delete trackCopy;
644 continue;
645 }
646 }
647 else if (fTrackType == kITSOnly) {
648 if (fConstrained==true)
649 tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
650 else
651 tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
652 }
653
654
655 AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
656 if (v.Mag() < 0.0001) {
657 // cout << "Found 0 momentum ???? " <<endl;
658 delete trackCopy;
659 continue;
660 }
661 trackCopy->SetP(v);//setting momentum
662 trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
663 const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
664 const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
665 //setting helix I do not if it is ok
666 AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
667 trackCopy->SetHelix(helix);
668
669 //some stuff which could be useful
670 float imp[2];
671 float cim[3];
672 // if (fTrackType == kTPCOnly) {
673 // esdtrack->GetTPCInnerParam()->GetImpactParameters(imp,cim);
674 // }
675 // else {
676 esdtrack->GetImpactParameters(imp,cim);
677 // }
678
679 impact[0] = imp[0];
680 impact[1] = imp[1];
681 covimpact[0] = cim[0];
682 covimpact[1] = cim[1];
683 covimpact[2] = cim[2];
684
685 trackCopy->SetImpactD(impact[0]);
686 trackCopy->SetImpactZ(impact[1]);
687 trackCopy->SetCdd(covimpact[0]);
688 trackCopy->SetCdz(covimpact[1]);
689 trackCopy->SetCzz(covimpact[2]);
690 trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact,covimpact));
691 }
692
693 trackCopy->SetTrackId(esdtrack->GetID());
694 trackCopy->SetFlags(esdtrack->GetStatus());
695 trackCopy->SetLabel(esdtrack->GetLabel());
696
697 trackCopy->SetITSchi2(esdtrack->GetITSchi2());
698 if (esdtrack->GetITSFakeFlag())
699 trackCopy->SetITSncls(-esdtrack->GetNcls(0));
700 else
701 trackCopy->SetITSncls(esdtrack->GetNcls(0));
702 trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());
703 trackCopy->SetTPCncls(esdtrack->GetTPCNcls());
704 trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());
705 trackCopy->SetTPCsignal(esdtrack->GetTPCsignal());
706 trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class
707 trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma());
708
709 trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap());
710 trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap());
711
712 double xtpc[3];
713 esdtrack->GetInnerXYZ(xtpc);
714 xtpc[2] -= fV1[2];
715 trackCopy->SetNominalTPCEntrancePoint(xtpc);
716
717 esdtrack->GetOuterXYZ(xtpc);
718 xtpc[2] -= fV1[2];
719 trackCopy->SetNominalTPCExitPoint(xtpc);
720
721 int indexes[3];
722 for (int ik=0; ik<3; ik++) {
723 indexes[ik] = esdtrack->GetKinkIndex(ik);
724 }
725 trackCopy->SetKinkIndexes(indexes);
4eac0b05 726
727 for (int ii=0; ii<6; ii++){
728 trackCopy->SetITSHitOnLayer(ii,esdtrack->HasPointOnITSLayer(ii));
729 }
730
76ce4b5b 731 //decision if we want this track
732 //if we using diffrent labels we want that this label was use for first time
733 //if we use hidden info we want to have match between sim data and ESD
734 if (tGoodMomentum==true)
735 {
736 hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
737 realnofTracks++;//real number of tracks
738 // delete trackCopy;
739 }
740 else
741 {
742 delete trackCopy;
743 }
744
745 }
746
747 hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
748
749 AliCentrality *cent = fEvent->GetCentrality();
750 if (cent) {
751 hbtEvent->SetCentralityV0(cent->GetCentralityPercentile("V0M"));
18757d69 752 hbtEvent->SetCentralityV0A(cent->GetCentralityPercentile("V0A"));
753 hbtEvent->SetCentralityV0C(cent->GetCentralityPercentile("V0C"));
c863b24a 754 hbtEvent->SetCentralityZNA(cent->GetCentralityPercentile("ZNA"));
18757d69 755 hbtEvent->SetCentralityZNC(cent->GetCentralityPercentile("ZNC"));
c863b24a 756 hbtEvent->SetCentralityCL1(cent->GetCentralityPercentile("CL1"));
18757d69 757 hbtEvent->SetCentralityCL0(cent->GetCentralityPercentile("CL0"));
758 hbtEvent->SetCentralityTKL(cent->GetCentralityPercentile("TKL"));
759 hbtEvent->SetCentralityFMD(cent->GetCentralityPercentile("FMD"));
760 hbtEvent->SetCentralityFMD(cent->GetCentralityPercentile("NPA"));
2bcd4b96 761 // hbtEvent->SetCentralitySPD1(cent->GetCentralityPercentile("CL1"));
607aa748 762 hbtEvent->SetCentralityTrk(cent->GetCentralityPercentile("TRK"));
763 hbtEvent->SetCentralityCND(cent->GetCentralityPercentile("CND"));
76ce4b5b 764
765 if (Debug()>1) printf(" FemtoReader Got Event with %f %f %f %f\n", cent->GetCentralityPercentile("V0M"), 0.0, cent->GetCentralityPercentile("CL1"), 0.0);
766 }
767
768 if (fEstEventMult == kGlobalCount)
769 hbtEvent->SetNormalizedMult(tNormMult);
ce7b3d98 770 else if (fEstEventMult == kReferenceITSTPC)
7eb9f5d0 771 hbtEvent->SetNormalizedMult(AliESDtrackCuts::GetReferenceMultiplicity(fEvent,AliESDtrackCuts::kTrackletsITSTPC,1.0));
ce7b3d98 772 else if(fEstEventMult == kReferenceITSSA)
7eb9f5d0 773 hbtEvent->SetNormalizedMult(AliESDtrackCuts::GetReferenceMultiplicity(fEvent,AliESDtrackCuts::kTrackletsITSSA,1.0));
ce7b3d98 774 else if(fEstEventMult == kReferenceTracklets)
7eb9f5d0 775 hbtEvent->SetNormalizedMult(AliESDtrackCuts::GetReferenceMultiplicity(fEvent,AliESDtrackCuts::kTracklets,1.0));
76ce4b5b 776 else if (fEstEventMult == kSPDLayer1)
777 hbtEvent->SetNormalizedMult(fEvent->GetMultiplicity()->GetNumberOfITSClusters(1));
4eac0b05 778 else if (fEstEventMult == kVZERO)
779 {
780 Float_t multV0 = 0;
781 for (Int_t i=0; i<64; i++)
782 multV0 += fEvent->GetVZEROData()->GetMultiplicity(i);
783 hbtEvent->SetNormalizedMult(multV0);
784 }
1938eadf 785 else if (fEstEventMult == kCentrality) {
76ce4b5b 786 // centrality between 0 (central) and 1 (very peripheral)
787
788 if (cent) {
789 if (cent->GetCentralityPercentile("V0M") < 0.00001)
790 hbtEvent->SetNormalizedMult(-1);
791 else
792 hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("V0M")));
793 if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
794 10.0*cent->GetCentralityPercentile("V0M"), lrint(10.0*cent->GetCentralityPercentile("V0M")));
795 }
796 }
18757d69 797 else if (fEstEventMult == kCentralityV0A) {
798 // centrality between 0 (central) and 1 (very peripheral)
799
800 if (cent) {
801 if (cent->GetCentralityPercentile("V0A") < 0.00001)
802 hbtEvent->SetNormalizedMult(-1);
803 else
804 hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("V0A")));
805 if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
806 10.0*cent->GetCentralityPercentile("V0A"), lrint(10.0*cent->GetCentralityPercentile("V0A")));
807 }
808 }
809 else if (fEstEventMult == kCentralityV0C) {
810 // centrality between 0 (central) and 1 (very peripheral)
811
812 if (cent) {
813 if (cent->GetCentralityPercentile("V0C") < 0.00001)
814 hbtEvent->SetNormalizedMult(-1);
815 else
816 hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("V0C")));
817 if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
818 10.0*cent->GetCentralityPercentile("V0C"), lrint(10.0*cent->GetCentralityPercentile("V0C")));
819 }
820 }
c863b24a 821 else if (fEstEventMult == kCentralityZNA) {
822 // centrality between 0 (central) and 1 (very peripheral)
823
824 if (cent) {
825 if (cent->GetCentralityPercentile("ZNA") < 0.00001)
826 hbtEvent->SetNormalizedMult(-1);
827 else
828 hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("ZNA")));
829 if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
830 10.0*cent->GetCentralityPercentile("ZNA"), lrint(10.0*cent->GetCentralityPercentile("ZNA")));
831 }
832 }
18757d69 833 else if (fEstEventMult == kCentralityZNC) {
834 // centrality between 0 (central) and 1 (very peripheral)
835
836 if (cent) {
837 if (cent->GetCentralityPercentile("ZNC") < 0.00001)
838 hbtEvent->SetNormalizedMult(-1);
839 else
840 hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("ZNC")));
841 if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
842 10.0*cent->GetCentralityPercentile("ZNC"), lrint(10.0*cent->GetCentralityPercentile("ZNC")));
843 }
844 }
c863b24a 845 else if (fEstEventMult == kCentralityCL1) {
846 // centrality between 0 (central) and 1 (very peripheral)
847
848 if (cent) {
849 if (cent->GetCentralityPercentile("CL1") < 0.00001)
850 hbtEvent->SetNormalizedMult(-1);
851 else
852 hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("CL1")));
853 if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
854 10.0*cent->GetCentralityPercentile("CL1"), lrint(10.0*cent->GetCentralityPercentile("CL1")));
855 }
856 }
18757d69 857 else if (fEstEventMult == kCentralityCL0) {
858 // centrality between 0 (central) and 1 (very peripheral)
859
860 if (cent) {
861 if (cent->GetCentralityPercentile("CL0") < 0.00001)
862 hbtEvent->SetNormalizedMult(-1);
863 else
864 hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("CL0")));
865 if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
866 10.0*cent->GetCentralityPercentile("CL0"), lrint(10.0*cent->GetCentralityPercentile("CL0")));
867 }
868 }
607aa748 869 else if (fEstEventMult == kCentralityTRK) {
870 // centrality between 0 (central) and 1 (very peripheral)
871
872 if (cent) {
873 if (cent->GetCentralityPercentile("TRK") < 0.00001)
874 hbtEvent->SetNormalizedMult(-1);
875 else
876 hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("TRK")));
877 if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
878 10.0*cent->GetCentralityPercentile("TRK"), lrint(10.0*cent->GetCentralityPercentile("TRK")));
879 }
880 }
18757d69 881 else if (fEstEventMult == kCentralityTKL) {
882 // centrality between 0 (central) and 1 (very peripheral)
883
884 if (cent) {
885 if (cent->GetCentralityPercentile("TKL") < 0.00001)
886 hbtEvent->SetNormalizedMult(-1);
887 else
888 hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("TKL")));
889 if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
890 10.0*cent->GetCentralityPercentile("TKL"), lrint(10.0*cent->GetCentralityPercentile("TKL")));
891 }
892 }
893 else if (fEstEventMult == kCentralityNPA) {
894 // centrality between 0 (central) and 1 (very peripheral)
895
896 if (cent) {
897 if (cent->GetCentralityPercentile("NPA") < 0.00001)
898 hbtEvent->SetNormalizedMult(-1);
899 else
900 hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("NPA")));
901 if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
902 10.0*cent->GetCentralityPercentile("NPA"), lrint(10.0*cent->GetCentralityPercentile("NPA")));
903 }
904 }
607aa748 905 else if (fEstEventMult == kCentralityCND) {
906 // centrality between 0 (central) and 1 (very peripheral)
907
908 if (cent) {
909 if (cent->GetCentralityPercentile("CND") < 0.00001)
910 hbtEvent->SetNormalizedMult(-1);
911 else
912 hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("CND")));
913 if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
914 10.0*cent->GetCentralityPercentile("CND"), lrint(10.0*cent->GetCentralityPercentile("CND")));
915 }
916 }
18757d69 917 else if (fEstEventMult == kCentralityFMD) {
918 // centrality between 0 (central) and 1 (very peripheral)
919
920 if (cent) {
921 if (cent->GetCentralityPercentile("FMD") < 0.00001)
922 hbtEvent->SetNormalizedMult(-1);
923 else
924 hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("FMD")));
925 if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(),
926 10.0*cent->GetCentralityPercentile("FMD"), lrint(10.0*cent->GetCentralityPercentile("FMD")));
927 }
928 }
607aa748 929
930
76ce4b5b 931
932 if (tNormMultPos > tNormMultNeg)
933 hbtEvent->SetZDCParticipants(tNormMultPos);
934 else
935 hbtEvent->SetZDCParticipants(tNormMultNeg);
936
937 AliEventplane* ep = fEvent->GetEventplane();
938 if (ep) {
939 hbtEvent->SetEP(ep);
940 hbtEvent->SetReactionPlaneAngle(ep->GetEventplane("Q"));
941 }
942
973a91f8 943 //V0
944 if(fReadV0)
945 {
946 for (Int_t i = 0; i < fEvent->GetNumberOfV0s(); i++) {
947 AliESDv0* esdv0 = fEvent->GetV0(i);
948 if (!esdv0) continue;
949 //if(esdv0->GetNDaughters()>2) continue;
950 //if(esdv0->GetNProngs()>2) continue;
951 if(esdv0->Charge()!=0) continue;
952 AliESDtrack *trackPos = fEvent->GetTrack(esdv0->GetPindex());
953 if(!trackPos) continue;
954 AliESDtrack *trackNeg = fEvent->GetTrack(esdv0->GetNindex());
955 if(!trackNeg) continue;
956 if(trackPos->Charge()==trackNeg->Charge()) continue;
957 AliFemtoV0* trackCopyV0 = new AliFemtoV0();
958 CopyESDtoFemtoV0(esdv0, trackCopyV0, fEvent);
959 hbtEvent->V0Collection()->push_back(trackCopyV0);
960 //cout<<"Pushback v0 to v0collection"<<endl;
961 }
962 }
963
964
76ce4b5b 965 if (Debug()>1) cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
76ce4b5b 966}
734c121a 967
76ce4b5b 968//___________________
969void AliFemtoEventReaderESDChain::SetESDSource(AliESDEvent *aESD)
970{
971 // The chain loads the ESD for us
972 // You must provide the address where it can be found
973 fEvent = aESD;
974}
975//___________________
976// void AliFemtoEventReaderESDChain::SetESDfriendSource(AliESDfriend *aFriend)
977// {
978// // We need the ESD tree to obtain
979// // information about the friend file location
980// fEventFriend = aFriend;
981// }
982
983//____________________________________________________________________
984Float_t AliFemtoEventReaderESDChain::GetSigmaToVertex(double *impact, double *covar)
985{
986 // Calculates the number of sigma to the vertex.
987
988 Float_t b[2];
989 Float_t bRes[2];
990 Float_t bCov[3];
991
992 b[0] = impact[0];
993 b[1] = impact[1];
994 bCov[0] = covar[0];
995 bCov[1] = covar[1];
996 bCov[2] = covar[2];
997
998 bRes[0] = TMath::Sqrt(bCov[0]);
999 bRes[1] = TMath::Sqrt(bCov[2]);
1000
1001 // -----------------------------------
1002 // How to get to a n-sigma cut?
1003 //
1004 // The accumulated statistics from 0 to d is
1005 //
1006 // -> Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
1007 // -> 1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
1008 //
1009 // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
1010 // Can this be expressed in a different way?
1011
1012 if (bRes[0] == 0 || bRes[1] ==0)
1013 return -1;
1014
1015 Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
1016
1017 // stupid rounding problem screws up everything:
1018 // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
1019 if (TMath::Exp(-d * d / 2) < 1e-10)
1020 return 1000;
1021
1022 d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
1023 return d;
1024}
1025
973a91f8 1026void AliFemtoEventReaderESDChain::CopyESDtoFemtoV0(AliESDv0 *tESDv0, AliFemtoV0 *tFemtoV0, AliESDEvent *tESDevent)
1027{
ce7b3d98 1028 double fPrimaryVtxPosition[3];
1029 tESDevent->GetPrimaryVertex()->GetXYZ(fPrimaryVtxPosition);
1030 tFemtoV0->SetdcaV0ToPrimVertex(tESDv0->GetD(fPrimaryVtxPosition[0],fPrimaryVtxPosition[1],fPrimaryVtxPosition[2]));
973a91f8 1031
1032 //tFemtoV0->SetdecayLengthV0(tESDv0->DecayLengthV0()); //wrocic do tego
1033 //tFemtoV0->SetdecayVertexV0X(tESDv0->DecayVertexV0X());
1034 //tFemtoV0->SetdecayVertexV0Y(tESDv0->DecayVertexV0Y());
1035 //tFemtoV0->SetdecayVertexV0Z(tESDv0->DecayVertexV0Z()); //nie ma w AliESDv0
1036 //AliFemtoThreeVector decayvertex(tESDv0->DecayVertexV0X(),tESDv0->DecayVertexV0Y(),tESDv0->DecayVertexV0Z());
1037 //tFemtoV0->SetdecayVertexV0(decayvertex);
1038 tFemtoV0->SetdcaV0Daughters(tESDv0->GetDcaV0Daughters());
973a91f8 1039 tFemtoV0->SetmomV0X(tESDv0->Px());
1040 tFemtoV0->SetmomV0Y(tESDv0->Py());
1041 tFemtoV0->SetmomV0Z(tESDv0->Pz());
1042 AliFemtoThreeVector momv0(tESDv0->Px(),tESDv0->Py(),tESDv0->Pz());
1043 tFemtoV0->SetmomV0(momv0);
1044 tFemtoV0->SetalphaV0(tESDv0->AlphaV0());
1045 tFemtoV0->SetptArmV0(tESDv0->PtArmV0());
1046 //tFemtoV0->SeteLambda(tESDv0->ELambda());
1047 //tFemtoV0->SeteK0Short(tESDv0->EK0Short());
1048 //tFemtoV0->SetePosProton(tESDv0->EPosProton());
1049 //tFemtoV0->SeteNegProton(tESDv0->ENegProton());
1050 tFemtoV0->SetmassLambda(tESDv0->GetEffMass(4,2));
1051 tFemtoV0->SetmassAntiLambda(tESDv0->GetEffMass(2,4));
1052 tFemtoV0->SetmassK0Short(tESDv0->GetEffMass(2,2));
1053 //tFemtoV0->SetrapLambda(tESDv0->RapLambda());
1054 //tFemtoV0->SetrapK0Short(tESDv0->RapK0Short());
973a91f8 1055 tFemtoV0->SetptV0(tESDv0->Pt());
1056 tFemtoV0->SetptotV0(tESDv0->P());
973a91f8 1057 tFemtoV0->SetEtaV0(tESDv0->Eta());
1058 tFemtoV0->SetPhiV0(tESDv0->Phi());
1059 tFemtoV0->SetCosPointingAngle(tESDv0->GetV0CosineOfPointingAngle(fPrimaryVtxPosition[0],fPrimaryVtxPosition[1], fPrimaryVtxPosition[2]));
ce7b3d98 1060 tFemtoV0->SetYV0(tESDv0->Y());
973a91f8 1061
ce7b3d98 1062 AliESDtrack *trackpos = tESDevent->GetTrack(tESDv0->GetPindex()); //AliAODTrack *trackpos = (AliAODTrack*)tESDv0->GetDaughter(0);
1063 AliESDtrack *trackneg = tESDevent->GetTrack(tESDv0->GetNindex()); //AliAODTrack *trackneg = (AliAODTrack*)tESDv0->GetDaughter(1);
973a91f8 1064
ce7b3d98 1065 if(trackpos && trackneg)
973a91f8 1066 {
ce7b3d98 1067 tFemtoV0->SetdcaPosToPrimVertex(TMath::Abs(trackpos->GetD(fPrimaryVtxPosition[0],fPrimaryVtxPosition[1],tESDevent->GetMagneticField())));
1068 tFemtoV0->SetdcaNegToPrimVertex(TMath::Abs(trackneg->GetD(fPrimaryVtxPosition[0],fPrimaryVtxPosition[1],tESDevent->GetMagneticField())));
1e100162 1069 double MomPos[3];
ce7b3d98 1070 trackpos->PxPyPz(MomPos);
1e100162 1071 tFemtoV0->SetmomPosX(MomPos[0]);
1072 tFemtoV0->SetmomPosY(MomPos[1]);
1073 tFemtoV0->SetmomPosZ(MomPos[2]);
1074 AliFemtoThreeVector mompos(MomPos[0],MomPos[1],MomPos[2]);
1075 tFemtoV0->SetmomPos(mompos);
1076
1077 double MomNeg[3];
ce7b3d98 1078 trackneg->PxPyPz(MomNeg);
1e100162 1079 tFemtoV0->SetmomNegX(MomNeg[0]);
1080 tFemtoV0->SetmomNegY(MomNeg[1]);
1081 tFemtoV0->SetmomNegZ(MomNeg[2]);
1082 AliFemtoThreeVector momneg(MomNeg[0],MomNeg[1],MomNeg[2]);
1083 tFemtoV0->SetmomNeg(momneg);
1084
ce7b3d98 1085 tFemtoV0->SetptPos(trackpos->Pt());
1086 tFemtoV0->SetptotPos(trackpos->P());
1087 tFemtoV0->SetptNeg(trackneg->Pt());
1088 tFemtoV0->SetptotNeg(trackneg->P());
1e100162 1089
ce7b3d98 1090 tFemtoV0->SetidNeg(trackneg->GetID());
1e100162 1091 //cout<<"tESDv0->GetNegID(): "<<tESDv0->GetNegID()<<endl;
1092 //cout<<"tFemtoV0->IdNeg(): "<<tFemtoV0->IdNeg()<<endl;
ce7b3d98 1093 tFemtoV0->SetidPos(trackpos->GetID());
1e100162 1094
ce7b3d98 1095 tFemtoV0->SetEtaPos(trackpos->Eta());
1096 tFemtoV0->SetEtaNeg(trackneg->Eta());
1097
1098 tFemtoV0->SetEtaPos(trackpos->Eta()); //tESDv0->PseudoRapPos()
1099 tFemtoV0->SetEtaNeg(trackneg->Eta()); //tESDv0->PseudoRapNeg()
1100 tFemtoV0->SetTPCNclsPos(trackpos->GetTPCNcls());
1101 tFemtoV0->SetTPCNclsNeg(trackneg->GetTPCNcls());
1102 tFemtoV0->SetTPCclustersPos(trackpos->GetTPCClusterMap());
1103 tFemtoV0->SetTPCclustersNeg(trackneg->GetTPCClusterMap());
1104 tFemtoV0->SetTPCsharingPos(trackpos->GetTPCSharedMap());
1105 tFemtoV0->SetTPCsharingNeg(trackneg->GetTPCSharedMap());
1106 tFemtoV0->SetNdofPos(trackpos->GetTPCchi2()/trackpos->GetTPCNcls());
1107 tFemtoV0->SetNdofNeg(trackneg->GetTPCchi2()/trackneg->GetTPCNcls());
1108 tFemtoV0->SetStatusPos(trackpos->GetStatus());
1109 tFemtoV0->SetStatusNeg(trackneg->GetStatus());
1110
1111 float bfield = 5*fMagFieldSign;
1112 float globalPositionsAtRadiiPos[9][3];
1113 GetGlobalPositionAtGlobalRadiiThroughTPC(trackpos,bfield,globalPositionsAtRadiiPos);
1114 double tpcEntrancePos[3]={globalPositionsAtRadiiPos[0][0],globalPositionsAtRadiiPos[0][1],globalPositionsAtRadiiPos[0][2]};
7eb9f5d0 1115 double tpcExitPos[3]={globalPositionsAtRadiiPos[7][0],globalPositionsAtRadiiPos[7][1],globalPositionsAtRadiiPos[7][2]};
ce7b3d98 1116
1117 float globalPositionsAtRadiiNeg[9][3];
1118 GetGlobalPositionAtGlobalRadiiThroughTPC(trackneg,bfield,globalPositionsAtRadiiNeg);
1119 double tpcEntranceNeg[3]={globalPositionsAtRadiiNeg[0][0],globalPositionsAtRadiiNeg[0][1],globalPositionsAtRadiiNeg[0][2]};
7eb9f5d0 1120 double tpcExitNeg[3]={globalPositionsAtRadiiNeg[7][0],globalPositionsAtRadiiNeg[7][1],globalPositionsAtRadiiNeg[7][2]};
ce7b3d98 1121
1122 AliFemtoThreeVector tmpVec;
1123 tmpVec.SetX(tpcEntrancePos[0]); tmpVec.SetX(tpcEntrancePos[1]); tmpVec.SetX(tpcEntrancePos[2]);
ba3c23a4 1124 tFemtoV0->SetNominalTpcEntrancePointPos(tmpVec);
ba3c23a4 1125
ce7b3d98 1126 tmpVec.SetX(tpcExitPos[0]); tmpVec.SetX(tpcExitPos[1]); tmpVec.SetX(tpcExitPos[2]);
ba3c23a4 1127 tFemtoV0->SetNominalTpcExitPointPos(tmpVec);
ce7b3d98 1128
1129 tmpVec.SetX(tpcEntranceNeg[0]); tmpVec.SetX(tpcEntranceNeg[1]); tmpVec.SetX(tpcEntranceNeg[2]);
1130 tFemtoV0->SetNominalTpcEntrancePointNeg(tmpVec);
1131
1132 tmpVec.SetX(tpcExitNeg[0]); tmpVec.SetX(tpcExitNeg[1]); tmpVec.SetX(tpcExitNeg[2]);
ba3c23a4 1133 tFemtoV0->SetNominalTpcExitPointNeg(tmpVec);
1134
ce7b3d98 1135 AliFemtoThreeVector vecTpcPos[9];
1136 AliFemtoThreeVector vecTpcNeg[9];
1137 for(int i=0;i<9;i++)
1138 {
1139 vecTpcPos[i].SetX(globalPositionsAtRadiiPos[i][0]); vecTpcPos[i].SetY(globalPositionsAtRadiiPos[i][1]); vecTpcPos[i].SetZ(globalPositionsAtRadiiPos[i][2]);
1140 vecTpcNeg[i].SetX(globalPositionsAtRadiiNeg[i][0]); vecTpcNeg[i].SetY(globalPositionsAtRadiiNeg[i][1]); vecTpcNeg[i].SetZ(globalPositionsAtRadiiNeg[i][2]);
1141 }
1142 tFemtoV0->SetNominalTpcPointPos(vecTpcPos);
1143 tFemtoV0->SetNominalTpcPointNeg(vecTpcNeg);
1144
1145 tFemtoV0->SetTPCMomentumPos(trackpos->GetTPCInnerParam()->P()); //trackpos->GetTPCmomentum();
1146 tFemtoV0->SetTPCMomentumNeg(trackneg->GetTPCInnerParam()->P()); //trackneg->GetTPCmomentum();
1147
1148 tFemtoV0->SetdedxPos(trackpos->GetTPCsignal());
1149 tFemtoV0->SetdedxNeg(trackneg->GetTPCsignal());
1150
ba3c23a4 1151
ea7124dd 1152 if (fESDpid) {
ce7b3d98 1153 tFemtoV0->SetPosNSigmaTPCK(fESDpid->NumberOfSigmasTPC(trackpos,AliPID::kKaon));
1154 tFemtoV0->SetNegNSigmaTPCK(fESDpid->NumberOfSigmasTPC(trackneg,AliPID::kKaon));
1155 tFemtoV0->SetPosNSigmaTPCP(fESDpid->NumberOfSigmasTPC(trackpos,AliPID::kProton));
1156 tFemtoV0->SetNegNSigmaTPCP(fESDpid->NumberOfSigmasTPC(trackneg,AliPID::kProton));
1157 tFemtoV0->SetPosNSigmaTPCPi(fESDpid->NumberOfSigmasTPC(trackpos,AliPID::kPion));
1158 tFemtoV0->SetNegNSigmaTPCPi(fESDpid->NumberOfSigmasTPC(trackneg,AliPID::kPion));
ea7124dd 1159 }
1160 else {
ce7b3d98 1161 tFemtoV0->SetPosNSigmaTPCK(-1000);
1162 tFemtoV0->SetNegNSigmaTPCK(-1000);
1163 tFemtoV0->SetPosNSigmaTPCP(-1000);
1164 tFemtoV0->SetNegNSigmaTPCP(-1000);
1165 tFemtoV0->SetPosNSigmaTPCPi(-1000);
1166 tFemtoV0->SetNegNSigmaTPCPi(-1000);
ea7124dd 1167 }
973a91f8 1168
2952eb1d 1169 if(// (tFemtoV0->StatusPos()&AliESDtrack::kTOFpid)==0 ||
1170 (tFemtoV0->StatusPos()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusPos()&AliESDtrack::kTOFout)==0)
973a91f8 1171 {
2952eb1d 1172 if(// (tFemtoV0->StatusNeg()&AliESDtrack::kTOFpid)==0 ||
1173 (tFemtoV0->StatusNeg()&AliESDtrack::kTIME)==0 || (tFemtoV0->StatusNeg()&AliESDtrack::kTOFout)==0)
973a91f8 1174 {
ce7b3d98 1175 tFemtoV0->SetPosNSigmaTOFK(-1000);
1176 tFemtoV0->SetNegNSigmaTOFK(-1000);
1177 tFemtoV0->SetPosNSigmaTOFP(-1000);
1178 tFemtoV0->SetNegNSigmaTOFP(-1000);
1179 tFemtoV0->SetPosNSigmaTOFPi(-1000);
1180 tFemtoV0->SetNegNSigmaTOFPi(-1000);
973a91f8 1181 }
1182 }
1183 else
1184 {
ea7124dd 1185 if (fESDpid) {
ce7b3d98 1186 tFemtoV0->SetPosNSigmaTOFK(fESDpid->NumberOfSigmasTOF(trackpos,AliPID::kKaon));
1187 tFemtoV0->SetNegNSigmaTOFK(fESDpid->NumberOfSigmasTOF(trackneg,AliPID::kKaon));
1188 tFemtoV0->SetPosNSigmaTOFP(fESDpid->NumberOfSigmasTOF(trackpos,AliPID::kProton));
1189 tFemtoV0->SetNegNSigmaTOFP(fESDpid->NumberOfSigmasTOF(trackneg,AliPID::kProton));
1190 tFemtoV0->SetPosNSigmaTOFPi(fESDpid->NumberOfSigmasTOF(trackpos,AliPID::kPion));
1191 tFemtoV0->SetNegNSigmaTOFPi(fESDpid->NumberOfSigmasTOF(trackneg,AliPID::kPion));
ea7124dd 1192 }
1193 else {
ce7b3d98 1194 tFemtoV0->SetPosNSigmaTOFK(-1000);
1195 tFemtoV0->SetNegNSigmaTOFK(-1000);
1196 tFemtoV0->SetPosNSigmaTOFP(-1000);
1197 tFemtoV0->SetNegNSigmaTOFP(-1000);
1198 tFemtoV0->SetPosNSigmaTOFPi(-1000);
1199 tFemtoV0->SetNegNSigmaTOFPi(-1000);
ea7124dd 1200 }
973a91f8 1201 }
1202 }
1203 else
1204 {
1205 tFemtoV0->SetStatusPos(999);
1206 tFemtoV0->SetStatusNeg(999);
1207 }
1208 tFemtoV0->SetOnFlyStatusV0(tESDv0->GetOnFlyStatus());
1209}
1210
76ce4b5b 1211void AliFemtoEventReaderESDChain::SetReadTrackType(ReadTrackType aType)
1212{
1213 fTrackType = aType;
1214}
1215
1216//trigger
1217void AliFemtoEventReaderESDChain::SetEventTrigger(UInt_t eventtrig)
1218{
1219 fEventTrig = eventtrig;
1220}
1221
973a91f8 1222//V0 reading
1223void AliFemtoEventReaderESDChain::SetReadV0(bool a)
1224{
1225 fReadV0 = a;
1226}
76ce4b5b 1227
ce7b3d98 1228void AliFemtoEventReaderESDChain::SetMagneticFieldSign(int s)
1229{
1230 if(s>0)
1231 fMagFieldSign = 1;
1232 else if(s<0)
1233 fMagFieldSign = -1;
1234 else
1235 fMagFieldSign = 0;
1236}
1237
1238void AliFemtoEventReaderESDChain::GetGlobalPositionAtGlobalRadiiThroughTPC(AliESDtrack *track, Float_t bfield, Float_t globalPositionsAtRadii[9][3])
1239{
1240 // Gets the global position of the track at nine different radii in the TPC
1241 // track is the track you want to propagate
1242 // bfield is the magnetic field of your event
1243 // globalPositionsAtRadii is the array of global positions in the radii and xyz
1244
1245 // Initialize the array to something indicating there was no propagation
1246 for(Int_t i=0;i<9;i++){
1247 for(Int_t j=0;j<3;j++){
1248 globalPositionsAtRadii[i][j]=-9999.;
1249 }
1250 }
1251
1252 // Make a copy of the track to not change parameters of the track
1253 AliExternalTrackParam etp; etp.CopyFromVTrack(track);
1254 //printf("\nAfter CopyFromVTrack\n");
1255 //etp.Print();
1256
1257 // The global position of the the track
1258 Double_t xyz[3]={-9999.,-9999.,-9999.};
1259
1260 // Counter for which radius we want
1261 Int_t iR=0;
1262 // The radii at which we get the global positions
1263 // IROC (OROC) from 84.1 cm to 132.1 cm (134.6 cm to 246.6 cm)
1264 Float_t Rwanted[9]={85.,105.,125.,145.,165.,185.,205.,225.,245.};
1265 // The global radius we are at
1266 Float_t globalRadius=0;
1267
1268 // Propagation is done in local x of the track
1269 for (Float_t x = etp.GetX();x<247.;x+=1.){ // GetX returns local coordinates
1270 // Starts at the tracks fX and goes outwards. x = 245 is the outer radial limit
1271 // of the TPC when the track is straight, i.e. has inifinite pt and doesn't get bent.
1272 // If the track's momentum is smaller than infinite, it will develop a y-component, which
1273 // adds to the global radius
1274
1275 // Stop if the propagation was not succesful. This can happen for low pt tracks
1276 // that don't reach outer radii
1277 if(!etp.PropagateTo(x,bfield))break;
1278 etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1279 globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
1280
1281 // Roughly reached the radius we want
1282 if(globalRadius > Rwanted[iR]){
1283
1284 // Bigger loop has bad precision, we're nearly one centimeter too far, go back in small steps.
1285 while (globalRadius>Rwanted[iR]){
1286 x-=.1;
1287 // printf("propagating to x %5.2f\n",x);
1288 if(!etp.PropagateTo(x,bfield))break;
1289 etp.GetXYZ(xyz); // GetXYZ returns global coordinates
1290 globalRadius = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); //Idea to speed up: compare squared radii
1291 }
1292 //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]);
1293 globalPositionsAtRadii[iR][0]=xyz[0];
1294 globalPositionsAtRadii[iR][1]=xyz[1];
1295 globalPositionsAtRadii[iR][2]=xyz[2];
1296 // Indicate we want the next radius
1297 iR+=1;
1298 }
1299 if(iR>=8){
1300 // TPC edge reached
1301 return;
1302 }
1303 }
1304}
76ce4b5b 1305
1306
734c121a 1307void AliFemtoEventReaderESDChain::SetpA2013(Bool_t pA2013)
1308{
1309 fpA2013 = pA2013;
1310}
7dec7414 1311
1312void AliFemtoEventReaderESDChain::SetUseMVPlpSelection(Bool_t mvplp)
1313{
1314 fMVPlp = mvplp;
1315}
1316
1317void AliFemtoEventReaderESDChain::SetIsPileUpEvent(Bool_t ispileup)
1318{
1319 fisPileUp = ispileup;
1320}