1 ////////////////////////////////////////////////////////////////////////////////
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 ///
8 ////////////////////////////////////////////////////////////////////////////////
9 #include "AliFemtoEventReaderESDChain.h"
13 #include "AliESDEvent.h"
14 #include "AliESDtrack.h"
15 #include "AliESDVertex.h"
17 #include "AliFmPhysicalHelixD.h"
18 #include "AliFmThreeVectorF.h"
20 #include "SystemOfUnits.h"
22 #include "AliFemtoEvent.h"
23 #include "AliFemtoModelHiddenInfo.h"
25 ClassImp(AliFemtoEventReaderESDChain)
27 #if !(ST_NO_NAMESPACES)
28 using namespace units;
32 //____________________________
33 AliFemtoEventReaderESDChain::AliFemtoEventReaderESDChain():
42 fUsePhysicsSel(kFALSE),
46 //constructor with 0 parameters , look at default settings
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>();
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>();
58 AliFemtoEventReaderESDChain::AliFemtoEventReaderESDChain(const AliFemtoEventReaderESDChain& aReader):
59 AliFemtoEventReader(aReader),
68 fUsePhysicsSel(kFALSE),
73 fConstrained = aReader.fConstrained;
74 fReadInner = aReader.fReadInner;
75 fUseTPCOnly = aReader.fUseTPCOnly;
76 fNumberofEvent = aReader.fNumberofEvent;
77 fCurEvent = aReader.fCurEvent;
78 fCurFile = aReader.fCurFile;
79 // fEvent = new AliESD(*aReader.fEvent);
80 fEvent = new AliESDEvent();
81 fUsePhysicsSel = aReader.fUsePhysicsSel;
82 if (aReader.fUsePhysicsSel)
83 fSelect = new AliPhysicsSelection();
84 fTrackType = aReader.fTrackType;
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);
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);
104 AliFemtoEventReaderESDChain::~AliFemtoEventReaderESDChain()
109 // for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
110 // fClusterPerPadrow[tPad]->clear();
111 // delete fClusterPerPadrow[tPad];
113 // delete [] fClusterPerPadrow;
114 // for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
115 // fSharedList[tPad]->clear();
116 // delete fSharedList[tPad];
118 // delete [] fSharedList;
119 if (fSelect) delete fSelect;
123 AliFemtoEventReaderESDChain& AliFemtoEventReaderESDChain::operator=(const AliFemtoEventReaderESDChain& aReader)
125 // Assignment operator
126 if (this == &aReader)
129 fConstrained = aReader.fConstrained;
130 fReadInner = aReader.fReadInner;
131 fUseTPCOnly = aReader.fUseTPCOnly;
132 fNumberofEvent = aReader.fNumberofEvent;
133 fCurEvent = aReader.fCurEvent;
134 fCurFile = aReader.fCurFile;
135 if (fEvent) delete fEvent;
136 fEvent = new AliESDEvent();
137 fTrackType = aReader.fTrackType;
139 fUsePhysicsSel = aReader.fUsePhysicsSel;
140 if (aReader.fUsePhysicsSel)
141 fSelect = new AliPhysicsSelection();
142 // fEventFriend = aReader.fEventFriend;
144 // if (fClusterPerPadrow) {
145 // for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
146 // fClusterPerPadrow[tPad]->clear();
147 // delete fClusterPerPadrow[tPad];
149 // delete [] fClusterPerPadrow;
152 // if (fSharedList) {
153 // for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
154 // fSharedList[tPad]->clear();
155 // delete fSharedList[tPad];
157 // delete [] fSharedList;
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);
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);
181 AliFemtoString AliFemtoEventReaderESDChain::Report()
183 AliFemtoString temp = "\n This is the AliFemtoEventReaderESDChain\n";
188 void AliFemtoEventReaderESDChain::SetConstrained(const bool constrained)
190 // Select whether to read constrained or not constrained momentum
191 fConstrained=constrained;
194 bool AliFemtoEventReaderESDChain::GetConstrained() const
196 // Check whether we read constrained or not constrained momentum
200 void AliFemtoEventReaderESDChain::SetReadTPCInner(const bool readinner)
202 fReadInner=readinner;
205 bool AliFemtoEventReaderESDChain::GetReadTPCInner() const
211 void AliFemtoEventReaderESDChain::SetUseTPCOnly(const bool usetpconly)
213 fUseTPCOnly=usetpconly;
216 bool AliFemtoEventReaderESDChain::GetUseTPCOnly() const
221 void AliFemtoEventReaderESDChain::SetUsePhysicsSelection(const bool usephysics)
223 fUsePhysicsSel = usephysics;
224 if (!fSelect) fSelect = new AliPhysicsSelection();
227 AliFemtoEvent* AliFemtoEventReaderESDChain::ReturnHbtEvent()
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;
235 // Get the friend information
236 cout<<"starting to read event "<<fCurEvent<<endl;
237 // fEvent->SetESDfriend(fEventFriend);
238 if(fEvent->GetAliESDOld())fEvent->CopyFromOldESD();
240 hbtEvent = new AliFemtoEvent;
242 if (fUsePhysicsSel) {
243 hbtEvent->SetIsCollisionCandidate(fSelect->IsCollisionCandidate(fEvent));
244 if (!(fSelect->IsCollisionCandidate(fEvent)))
245 printf("Event not a collision candidate\n");
248 hbtEvent->SetIsCollisionCandidate(kTRUE);
250 //setting basic things
251 // hbtEvent->SetEventNumber(fEvent->GetEventNumber());
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());
268 fEvent->GetPrimaryVertexTPC()->GetXYZ(fV1);
269 fEvent->GetPrimaryVertexTPC()->GetCovMatrix(fVCov);
270 if (!fEvent->GetPrimaryVertexTPC()->GetStatus())
274 fEvent->GetPrimaryVertex()->GetXYZ(fV1);
275 fEvent->GetPrimaryVertex()->GetCovMatrix(fVCov);
276 if (!fEvent->GetPrimaryVertex()->GetStatus())
280 AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
281 hbtEvent->SetPrimVertPos(vertex);
282 hbtEvent->SetPrimVertCov(fVCov);
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
289 // // Clear the shared cluster list
290 // for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
291 // fClusterPerPadrow[tPad]->clear();
293 // for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
294 // fSharedList[tPad]->clear();
298 // for (int i=0;i<nofTracks;i++) {
299 // const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
301 // list<Int_t>::iterator tClustIter;
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]);
312 // fSharedList[tNcl]->push_back(tTrackIndices[tNcl]);
320 for (int i=0;i<nofTracks;i++)
322 bool tGoodMomentum=true; //flaga to chcek if we can read momentum of this track
324 const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
325 // const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
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)
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();
340 for (int iter=0; iter<6; iter++) if (iclm&(1<<iter)) incls++;
342 cout << "Rejecting track with " << incls << " clusters" << endl;
347 AliFemtoTrack* trackCopy = new AliFemtoTrack();
348 trackCopy->SetCharge((short)esdtrack->GetSign());
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
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]);
367 if (!esdtrack->GetTPCInnerParam()) {
373 AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
375 param->PropagateToDCA(fEvent->GetPrimaryVertexTPC(), (fEvent->GetMagneticField()), 10000, impact, covimpact);
376 param->GetPxPyPz(pxyz);//reading noconstarined momentum
378 if (fReadInner == true) {
379 AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
380 tInfo->SetPDGPid(211);
381 tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]);
382 tInfo->SetMass(0.13957);
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);
386 trackCopy->SetHiddenInfo(tInfo);
389 AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
390 if (v.Mag() < 0.0001) {
391 // cout << "Found 0 momentum ???? " <<endl;
395 trackCopy->SetP(v);//setting momentum
396 trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
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);
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));
415 if (fReadInner == true) {
417 if (esdtrack->GetTPCInnerParam()) {
418 AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
420 param->PropagateToDCA(fEvent->GetPrimaryVertex(), (fEvent->GetMagneticField()), 10000);
421 param->GetPxPyPz(pxyz);//reading noconstarined momentum
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);
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);
431 trackCopy->SetHiddenInfo(tInfo);
434 if (fConstrained==true)
435 tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
437 tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
439 AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
440 if (v.Mag() < 0.0001) {
441 // cout << "Found 0 momentum ???? " <<endl;
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);
453 //some stuff which could be useful
456 esdtrack->GetImpactParameters(imp,cim);
460 covimpact[0] = cim[0];
461 covimpact[1] = cim[1];
462 covimpact[2] = cim[2];
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));
472 trackCopy->SetTrackId(esdtrack->GetID());
473 trackCopy->SetFlags(esdtrack->GetStatus());
474 trackCopy->SetLabel(esdtrack->GetLabel());
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());
484 trackCopy->SetTPCClusterMap(esdtrack->GetTPCClusterMap());
485 trackCopy->SetTPCSharedMap(esdtrack->GetTPCSharedMap());
488 esdtrack->GetInnerXYZ(xtpc);
490 trackCopy->SetNominalTPCEntrancePoint(xtpc);
492 esdtrack->GetOuterXYZ(xtpc);
494 trackCopy->SetNominalTPCExitPoint(xtpc);
497 for (int ik=0; ik<3; ik++) {
498 indexes[ik] = esdtrack->GetKinkIndex(ik);
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
504 if (tGoodMomentum==true)
506 hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
507 realnofTracks++;//real number of tracks
517 hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
518 hbtEvent->SetNormalizedMult(tNormMult);
520 cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
523 //___________________
524 void AliFemtoEventReaderESDChain::SetESDSource(AliESDEvent *aESD)
526 // The chain loads the ESD for us
527 // You must provide the address where it can be found
530 //___________________
531 // void AliFemtoEventReaderESDChain::SetESDfriendSource(AliESDfriend *aFriend)
533 // // We need the ESD tree to obtain
534 // // information about the friend file location
535 // fEventFriend = aFriend;
538 //____________________________________________________________________
539 Float_t AliFemtoEventReaderESDChain::GetSigmaToVertex(double *impact, double *covar)
541 // Calculates the number of sigma to the vertex.
553 bRes[0] = TMath::Sqrt(bCov[0]);
554 bRes[1] = TMath::Sqrt(bCov[2]);
556 // -----------------------------------
557 // How to get to a n-sigma cut?
559 // The accumulated statistics from 0 to d is
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)
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?
567 if (bRes[0] == 0 || bRes[1] ==0)
570 Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
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)
577 d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
581 void AliFemtoEventReaderESDChain::SetReadTrackType(ReadTrackType aType)