]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/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
Fix Gamma calculation^CliFemtoModelGausLCMSFreezeOutGenerator.cxx
[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemto / AliFemtoEventReaderESDChain.cxx
CommitLineData
d0e92d9a 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////////////////////////////////////////////////////////////////////////////////
67427ff7 9#include "AliFemtoEventReaderESDChain.h"
10
11#include "TFile.h"
12#include "TTree.h"
ea77036b 13#include "AliESDEvent.h"
67427ff7 14#include "AliESDtrack.h"
ea77036b 15#include "AliESDVertex.h"
67427ff7 16
d0e92d9a 17#include "AliFmPhysicalHelixD.h"
18#include "AliFmThreeVectorF.h"
67427ff7 19
d0e92d9a 20#include "SystemOfUnits.h"
67427ff7 21
d0e92d9a 22#include "AliFemtoEvent.h"
d8f11913 23#include "AliFemtoModelHiddenInfo.h"
67427ff7 24
25ClassImp(AliFemtoEventReaderESDChain)
26
27#if !(ST_NO_NAMESPACES)
28 using namespace units;
29#endif
30
31using namespace std;
32//____________________________
67427ff7 33AliFemtoEventReaderESDChain::AliFemtoEventReaderESDChain():
34 fFileName(" "),
35 fConstrained(true),
d8f11913 36 fReadInner(false),
537bd404 37 fUseTPCOnly(false),
67427ff7 38 fNumberofEvent(0),
39 fCurEvent(0),
40 fCurFile(0),
73695088 41 fEvent(0x0),
42 fUsePhysicsSel(kFALSE),
43 fSelect(0x0),
44 fTrackType(kGlobal)
67427ff7 45{
d0e92d9a 46 //constructor with 0 parameters , look at default settings
ea77036b 47// fClusterPerPadrow = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
48// for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
49// fClusterPerPadrow[tPad] = new list<Int_t>();
50// }
51// fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
52// for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
53// fSharedList[tPad] = new list<Int_t>();
54// }
67427ff7 55}
56
0215f606 57//__________________
0215f606 58AliFemtoEventReaderESDChain::AliFemtoEventReaderESDChain(const AliFemtoEventReaderESDChain& aReader):
fcda1d4e 59 AliFemtoEventReader(aReader),
0215f606 60 fFileName(" "),
61 fConstrained(true),
d8f11913 62 fReadInner(false),
537bd404 63 fUseTPCOnly(false),
0215f606 64 fNumberofEvent(0),
65 fCurEvent(0),
66 fCurFile(0),
73695088 67 fEvent(0x0),
68 fUsePhysicsSel(kFALSE),
69 fSelect(0x0),
70 fTrackType(kGlobal)
0215f606 71{
d0e92d9a 72 // Copy constructor
0215f606 73 fConstrained = aReader.fConstrained;
d8f11913 74 fReadInner = aReader.fReadInner;
537bd404 75 fUseTPCOnly = aReader.fUseTPCOnly;
0215f606 76 fNumberofEvent = aReader.fNumberofEvent;
77 fCurEvent = aReader.fCurEvent;
78 fCurFile = aReader.fCurFile;
d0e92d9a 79 // fEvent = new AliESD(*aReader.fEvent);
ea77036b 80 fEvent = new AliESDEvent();
73695088 81 fUsePhysicsSel = aReader.fUsePhysicsSel;
82 if (aReader.fUsePhysicsSel)
83 fSelect = new AliPhysicsSelection();
84 fTrackType = aReader.fTrackType;
ea77036b 85// fEventFriend = aReader.fEventFriend;
86// fClusterPerPadrow = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
87// for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
88// fClusterPerPadrow[tPad] = new list<Int_t>();
89// list<Int_t>::iterator iter;
90// for (iter=aReader.fClusterPerPadrow[tPad]->begin(); iter!=aReader.fClusterPerPadrow[tPad]->end(); iter++) {
91// fClusterPerPadrow[tPad]->push_back(*iter);
92// }
93// }
94// fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
95// for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
96// fSharedList[tPad] = new list<Int_t>();
97// list<Int_t>::iterator iter;
98// for (iter=aReader.fSharedList[tPad]->begin(); iter!=aReader.fSharedList[tPad]->end(); iter++) {
99// fSharedList[tPad]->push_back(*iter);
100// }
101// }
0215f606 102}
67427ff7 103//__________________
67427ff7 104AliFemtoEventReaderESDChain::~AliFemtoEventReaderESDChain()
105{
d0e92d9a 106 //Destructor
67427ff7 107 delete fEvent;
108
ea77036b 109// for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
110// fClusterPerPadrow[tPad]->clear();
111// delete fClusterPerPadrow[tPad];
112// }
113// delete [] fClusterPerPadrow;
114// for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
115// fSharedList[tPad]->clear();
116// delete fSharedList[tPad];
117// }
118// delete [] fSharedList;
73695088 119 if (fSelect) delete fSelect;
67427ff7 120}
121
0215f606 122//__________________
0215f606 123AliFemtoEventReaderESDChain& AliFemtoEventReaderESDChain::operator=(const AliFemtoEventReaderESDChain& aReader)
124{
d0e92d9a 125 // Assignment operator
0215f606 126 if (this == &aReader)
127 return *this;
128
129 fConstrained = aReader.fConstrained;
d8f11913 130 fReadInner = aReader.fReadInner;
537bd404 131 fUseTPCOnly = aReader.fUseTPCOnly;
0215f606 132 fNumberofEvent = aReader.fNumberofEvent;
133 fCurEvent = aReader.fCurEvent;
134 fCurFile = aReader.fCurFile;
135 if (fEvent) delete fEvent;
ea77036b 136 fEvent = new AliESDEvent();
73695088 137 fTrackType = aReader.fTrackType;
0215f606 138
73695088 139 fUsePhysicsSel = aReader.fUsePhysicsSel;
140 if (aReader.fUsePhysicsSel)
141 fSelect = new AliPhysicsSelection();
ea77036b 142 // fEventFriend = aReader.fEventFriend;
0215f606 143
ea77036b 144// if (fClusterPerPadrow) {
145// for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
146// fClusterPerPadrow[tPad]->clear();
147// delete fClusterPerPadrow[tPad];
148// }
149// delete [] fClusterPerPadrow;
150// }
0215f606 151
ea77036b 152// if (fSharedList) {
153// for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
154// fSharedList[tPad]->clear();
155// delete fSharedList[tPad];
156// }
157// delete [] fSharedList;
158// }
159
160// fClusterPerPadrow = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
161// for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
162// fClusterPerPadrow[tPad] = new list<Int_t>();
163// list<Int_t>::iterator iter;
164// for (iter=aReader.fClusterPerPadrow[tPad]->begin(); iter!=aReader.fClusterPerPadrow[tPad]->end(); iter++) {
165// fClusterPerPadrow[tPad]->push_back(*iter);
166// }
167// }
168// fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
169// for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
170// fSharedList[tPad] = new list<Int_t>();
171// list<Int_t>::iterator iter;
172// for (iter=aReader.fSharedList[tPad]->begin(); iter!=aReader.fSharedList[tPad]->end(); iter++) {
173// fSharedList[tPad]->push_back(*iter);
174// }
175// }
0215f606 176
177 return *this;
178}
67427ff7 179//__________________
180// Simple report
181AliFemtoString AliFemtoEventReaderESDChain::Report()
182{
183 AliFemtoString temp = "\n This is the AliFemtoEventReaderESDChain\n";
184 return temp;
185}
186
187//__________________
67427ff7 188void AliFemtoEventReaderESDChain::SetConstrained(const bool constrained)
189{
d0e92d9a 190 // Select whether to read constrained or not constrained momentum
67427ff7 191 fConstrained=constrained;
192}
193//__________________
67427ff7 194bool AliFemtoEventReaderESDChain::GetConstrained() const
195{
d0e92d9a 196 // Check whether we read constrained or not constrained momentum
67427ff7 197 return fConstrained;
198}
199//__________________
d8f11913 200void AliFemtoEventReaderESDChain::SetReadTPCInner(const bool readinner)
201{
202 fReadInner=readinner;
203}
204
205bool AliFemtoEventReaderESDChain::GetReadTPCInner() const
206{
207 return fReadInner;
208}
209
537bd404 210//__________________
211void AliFemtoEventReaderESDChain::SetUseTPCOnly(const bool usetpconly)
212{
213 fUseTPCOnly=usetpconly;
214}
215
216bool AliFemtoEventReaderESDChain::GetUseTPCOnly() const
217{
218 return fUseTPCOnly;
219}
220
73695088 221void AliFemtoEventReaderESDChain::SetUsePhysicsSelection(const bool usephysics)
222{
223 fUsePhysicsSel = usephysics;
224 if (!fSelect) fSelect = new AliPhysicsSelection();
225}
226
67427ff7 227AliFemtoEvent* AliFemtoEventReaderESDChain::ReturnHbtEvent()
228{
229 // Get the event, read all the relevant information
230 // and fill the AliFemtoEvent class
231 // Returns a valid AliFemtoEvent
232 AliFemtoEvent *hbtEvent = 0;
233 string tFriendFileName;
234
235 // Get the friend information
236 cout<<"starting to read event "<<fCurEvent<<endl;
ea77036b 237 // fEvent->SetESDfriend(fEventFriend);
376f40b1 238 if(fEvent->GetAliESDOld())fEvent->CopyFromOldESD();
73695088 239
67427ff7 240 hbtEvent = new AliFemtoEvent;
73695088 241
242 if (fUsePhysicsSel) {
243 hbtEvent->SetIsCollisionCandidate(fSelect->IsCollisionCandidate(fEvent));
244 if (!(fSelect->IsCollisionCandidate(fEvent)))
245 printf("Event not a collision candidate\n");
246 }
247 else
248 hbtEvent->SetIsCollisionCandidate(kTRUE);
249
67427ff7 250 //setting basic things
a531c6fc 251 // hbtEvent->SetEventNumber(fEvent->GetEventNumber());
67427ff7 252 hbtEvent->SetRunNumber(fEvent->GetRunNumber());
253 //hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
254 hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
255 hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
256 hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
257 hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
258 hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
259 hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
260 hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
261 hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
262 hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
263
264 //Vertex
265 double fV1[3];
63a5982a 266 double fVCov[6];
537bd404 267 if (fUseTPCOnly) {
268 fEvent->GetPrimaryVertexTPC()->GetXYZ(fV1);
63a5982a 269 fEvent->GetPrimaryVertexTPC()->GetCovMatrix(fVCov);
270 if (!fEvent->GetPrimaryVertexTPC()->GetStatus())
271 fVCov[4] = -1001.0;
537bd404 272 }
273 else {
274 fEvent->GetPrimaryVertex()->GetXYZ(fV1);
63a5982a 275 fEvent->GetPrimaryVertex()->GetCovMatrix(fVCov);
276 if (!fEvent->GetPrimaryVertex()->GetStatus())
277 fVCov[4] = -1001.0;
537bd404 278 }
67427ff7 279
280 AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
281 hbtEvent->SetPrimVertPos(vertex);
63a5982a 282 hbtEvent->SetPrimVertCov(fVCov);
67427ff7 283
284 //starting to reading tracks
285 int nofTracks=0; //number of reconstructed tracks in event
286 nofTracks=fEvent->GetNumberOfTracks();
287 int realnofTracks=0;//number of track which we use ina analysis
288
ea77036b 289// // Clear the shared cluster list
290// for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
291// fClusterPerPadrow[tPad]->clear();
292// }
293// for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
294// fSharedList[tPad]->clear();
295// }
296
297
298// for (int i=0;i<nofTracks;i++) {
299// const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
300
301// list<Int_t>::iterator tClustIter;
302
303// Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster];
304// Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices);
305// for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
306// if (tTrackIndices[tNcl] >= 0) {
307// tClustIter = find(fClusterPerPadrow[tNcl]->begin(), fClusterPerPadrow[tNcl]->end(), tTrackIndices[tNcl]);
308// if (tClustIter == fClusterPerPadrow[tNcl]->end()) {
309// fClusterPerPadrow[tNcl]->push_back(tTrackIndices[tNcl]);
310// }
311// else {
312// fSharedList[tNcl]->push_back(tTrackIndices[tNcl]);
313// }
314// }
315// }
67427ff7 316
ea77036b 317// }
67427ff7 318
03decc29 319 int tNormMult = 0;
67427ff7 320 for (int i=0;i<nofTracks;i++)
321 {
d0e92d9a 322 bool tGoodMomentum=true; //flaga to chcek if we can read momentum of this track
03decc29 323
67427ff7 324 const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
0215f606 325 // const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
03decc29 326 if (esdtrack->GetStatus() & AliESDtrack::kTPCrefit)
327 if (esdtrack->GetTPCNcls() > 80)
328 if (esdtrack->GetTPCchi2()/esdtrack->GetTPCNcls() < 6.0)
329 if (esdtrack->GetConstrainedParam())
330 if (esdtrack->GetConstrainedParam()->Eta() < 0.9)
331 tNormMult++;
67427ff7 332
73695088 333 // If reading ITS-only tracks, reject all with TPC
334 if (fTrackType == kITSOnly) {
335 if (esdtrack->GetStatus() & AliESDtrack::kTPCrefit) continue;
336 if (!(esdtrack->GetStatus() & AliESDtrack::kITSrefit)) continue;
337 if (esdtrack->GetStatus() & AliESDtrack::kTPCin) continue;
338 UChar_t iclm = esdtrack->GetITSClusterMap();
339 Int_t incls = 0;
340 for (int iter=0; iter<6; iter++) if (iclm&(1<<iter)) incls++;
341 if (incls<=3) {
342 cout << "Rejecting track with " << incls << " clusters" << endl;
343 continue;
344 }
345 }
346
347 AliFemtoTrack* trackCopy = new AliFemtoTrack();
67427ff7 348 trackCopy->SetCharge((short)esdtrack->GetSign());
349
350 //in aliroot we have AliPID
351 //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
352 //we use only 5 first
353 double esdpid[5];
354 esdtrack->GetESDpid(esdpid);
355 trackCopy->SetPidProbElectron(esdpid[0]);
356 trackCopy->SetPidProbMuon(esdpid[1]);
357 trackCopy->SetPidProbPion(esdpid[2]);
358 trackCopy->SetPidProbKaon(esdpid[3]);
359 trackCopy->SetPidProbProton(esdpid[4]);
360
361 double pxyz[3];
537bd404 362 double rxyz[3];
363 double impact[2];
364 double covimpact[3];
365
366 if (fUseTPCOnly) {
367 if (!esdtrack->GetTPCInnerParam()) {
368 delete trackCopy;
369 continue;
370 }
371
372
373 AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
374 param->GetXYZ(rxyz);
375 param->PropagateToDCA(fEvent->GetPrimaryVertexTPC(), (fEvent->GetMagneticField()), 10000, impact, covimpact);
376 param->GetPxPyPz(pxyz);//reading noconstarined momentum
d8f11913 377
537bd404 378 if (fReadInner == true) {
d8f11913 379 AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
380 tInfo->SetPDGPid(211);
381 tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]);
382 tInfo->SetMass(0.13957);
537bd404 383 // tInfo->SetEmissionPoint(rxyz[0], rxyz[1], rxyz[2], 0.0);
384 // tInfo->SetEmissionPoint(fV1[0], fV1[1], fV1[2], 0.0);
385 tInfo->SetEmissionPoint(rxyz[0]-fV1[0], rxyz[1]-fV1[1], rxyz[2]-fV1[2], 0.0);
d8f11913 386 trackCopy->SetHiddenInfo(tInfo);
387 }
537bd404 388
389 AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
390 if (v.mag() < 0.0001) {
391 // cout << "Found 0 momentum ???? " <<endl;
392 delete trackCopy;
393 continue;
394 }
395 trackCopy->SetP(v);//setting momentum
396 trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
397
398 const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
399 const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
400 //setting helix I do not if it is ok
401 AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
402 trackCopy->SetHelix(helix);
403
404 //some stuff which could be useful
405 trackCopy->SetImpactD(impact[0]);
406 trackCopy->SetImpactZ(impact[1]);
407 trackCopy->SetCdd(covimpact[0]);
408 trackCopy->SetCdz(covimpact[1]);
409 trackCopy->SetCzz(covimpact[2]);
410 trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact, covimpact));
411
412 delete param;
d8f11913 413 }
537bd404 414 else {
415 if (fReadInner == true) {
416
417 if (esdtrack->GetTPCInnerParam()) {
418 AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
419 param->GetXYZ(rxyz);
420 param->PropagateToDCA(fEvent->GetPrimaryVertex(), (fEvent->GetMagneticField()), 10000);
421 param->GetPxPyPz(pxyz);//reading noconstarined momentum
422 delete param;
423
424 AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
425 tInfo->SetPDGPid(211);
426 tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]);
427 tInfo->SetMass(0.13957);
428 // tInfo->SetEmissionPoint(rxyz[0], rxyz[1], rxyz[2], 0.0);
8597d9a1 429 //tInfo->SetEmissionPoint(fV1[0], fV1[1], fV1[2], 0.0);
430 tInfo->SetEmissionPoint(rxyz[0]-fV1[0], rxyz[1]-fV1[1], rxyz[2]-fV1[2], 0.0);
537bd404 431 trackCopy->SetHiddenInfo(tInfo);
432 }
433 }
434 if (fConstrained==true)
435 tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
436 else
437 tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
438
439 AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
440 if (v.mag() < 0.0001) {
441 // cout << "Found 0 momentum ???? " <<endl;
442 delete trackCopy;
443 continue;
444 }
445 trackCopy->SetP(v);//setting momentum
446 trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
447 const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
448 const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
449 //setting helix I do not if it is ok
450 AliFmPhysicalHelixD helix(kP,kOrigin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
451 trackCopy->SetHelix(helix);
452
453 //some stuff which could be useful
454 float imp[2];
455 float cim[3];
456 esdtrack->GetImpactParameters(imp,cim);
457
458 impact[0] = imp[0];
459 impact[1] = imp[1];
460 covimpact[0] = cim[0];
461 covimpact[1] = cim[1];
462 covimpact[2] = cim[2];
463
464 trackCopy->SetImpactD(impact[0]);
465 trackCopy->SetImpactZ(impact[1]);
466 trackCopy->SetCdd(covimpact[0]);
467 trackCopy->SetCdz(covimpact[1]);
468 trackCopy->SetCzz(covimpact[2]);
469 trackCopy->SetSigmaToVertex(GetSigmaToVertex(impact,covimpact));
67427ff7 470 }
537bd404 471
67427ff7 472 trackCopy->SetTrackId(esdtrack->GetID());
473 trackCopy->SetFlags(esdtrack->GetStatus());
d9241b1b 474 trackCopy->SetLabel(esdtrack->GetLabel());
67427ff7 475
67427ff7 476 trackCopy->SetITSchi2(esdtrack->GetITSchi2());
477 trackCopy->SetITSncls(esdtrack->GetNcls(0));
478 trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());
479 trackCopy->SetTPCncls(esdtrack->GetTPCNcls());
480 trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());
481 trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class
482 trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma());
ea77036b 483
484 trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap());
485 trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap());
67427ff7 486
0b3bd1ac 487 double xtpc[3];
488 esdtrack->GetInnerXYZ(xtpc);
537bd404 489 xtpc[2] -= fV1[2];
0b3bd1ac 490 trackCopy->SetNominalTPCEntrancePoint(xtpc);
491
492 esdtrack->GetOuterXYZ(xtpc);
537bd404 493 xtpc[2] -= fV1[2];
0b3bd1ac 494 trackCopy->SetNominalTPCExitPoint(xtpc);
495
496 int indexes[3];
497 for (int ik=0; ik<3; ik++) {
498 indexes[ik] = esdtrack->GetKinkIndex(ik);
499 }
500 trackCopy->SetKinkIndexes(indexes);
501 //decision if we want this track
502 //if we using diffrent labels we want that this label was use for first time
503 //if we use hidden info we want to have match between sim data and ESD
d0e92d9a 504 if (tGoodMomentum==true)
67427ff7 505 {
506 hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
507 realnofTracks++;//real number of tracks
508 // delete trackCopy;
509 }
510 else
511 {
512 delete trackCopy;
513 }
514
515 }
516
517 hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
03decc29 518 hbtEvent->SetNormalizedMult(tNormMult);
67427ff7 519 fCurEvent++;
520 cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
521 return hbtEvent;
522}
523//___________________
ea77036b 524void AliFemtoEventReaderESDChain::SetESDSource(AliESDEvent *aESD)
67427ff7 525{
d0e92d9a 526 // The chain loads the ESD for us
527 // You must provide the address where it can be found
67427ff7 528 fEvent = aESD;
529}
530//___________________
ea77036b 531// void AliFemtoEventReaderESDChain::SetESDfriendSource(AliESDfriend *aFriend)
532// {
533// // We need the ESD tree to obtain
534// // information about the friend file location
535// fEventFriend = aFriend;
536// }
67427ff7 537
0b3bd1ac 538//____________________________________________________________________
537bd404 539Float_t AliFemtoEventReaderESDChain::GetSigmaToVertex(double *impact, double *covar)
0b3bd1ac 540{
541 // Calculates the number of sigma to the vertex.
542
543 Float_t b[2];
544 Float_t bRes[2];
545 Float_t bCov[3];
537bd404 546
547 b[0] = impact[0];
548 b[1] = impact[1];
549 bCov[0] = covar[0];
550 bCov[1] = covar[1];
551 bCov[2] = covar[2];
552
0b3bd1ac 553 bRes[0] = TMath::Sqrt(bCov[0]);
554 bRes[1] = TMath::Sqrt(bCov[2]);
555
556 // -----------------------------------
557 // How to get to a n-sigma cut?
558 //
559 // The accumulated statistics from 0 to d is
560 //
561 // -> Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
562 // -> 1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
563 //
564 // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
565 // Can this be expressed in a different way?
566
567 if (bRes[0] == 0 || bRes[1] ==0)
568 return -1;
569
570 Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
571
572 // stupid rounding problem screws up everything:
573 // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
574 if (TMath::Exp(-d * d / 2) < 1e-10)
575 return 1000;
576
577 d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
578 return d;
579}
67427ff7 580
73695088 581void AliFemtoEventReaderESDChain::SetReadTrackType(ReadTrackType aType)
582{
583 fTrackType = aType;
584}
67427ff7 585
586
587
588
589
590
591