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