Modifications needed for the compilation of the femtoscopy code (Adam-Mike)
[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemto / Reader / AliFemtoEventReaderESDChain.cxx
CommitLineData
67427ff7 1#include "AliFemtoEventReaderESDChain.h"
2
3#include "TFile.h"
4#include "TTree.h"
5#include "AliESD.h"
6#include "AliESDtrack.h"
7
b2f60a91 8#include "Infrastructure/AliFmPhysicalHelixD.h"
9#include "Infrastructure/AliFmThreeVectorF.h"
67427ff7 10
b2f60a91 11#include "Base/SystemOfUnits.h"
67427ff7 12
b2f60a91 13#include "Infrastructure/AliFemtoEvent.h"
67427ff7 14
15ClassImp(AliFemtoEventReaderESDChain)
16
17#if !(ST_NO_NAMESPACES)
18 using namespace units;
19#endif
20
21using namespace std;
22//____________________________
23//constructor with 0 parameters , look at default settings
24AliFemtoEventReaderESDChain::AliFemtoEventReaderESDChain():
25 fFileName(" "),
26 fConstrained(true),
27 fNumberofEvent(0),
28 fCurEvent(0),
29 fCurFile(0),
30 fEvent(0x0),
31 fEventFriend(0)
32{
33 fClusterPerPadrow = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
34 for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
35 fClusterPerPadrow[tPad] = new list<Int_t>();
36 }
37 fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
38 for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
39 fSharedList[tPad] = new list<Int_t>();
40 }
41}
42
43//__________________
44//Destructor
45AliFemtoEventReaderESDChain::~AliFemtoEventReaderESDChain()
46{
47 delete fEvent;
48
49 for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
50 fClusterPerPadrow[tPad]->clear();
51 delete fClusterPerPadrow[tPad];
52 }
53 delete [] fClusterPerPadrow;
54 for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
55 fSharedList[tPad]->clear();
56 delete fSharedList[tPad];
57 }
58 delete [] fSharedList;
59}
60
61//__________________
62// Simple report
63AliFemtoString AliFemtoEventReaderESDChain::Report()
64{
65 AliFemtoString temp = "\n This is the AliFemtoEventReaderESDChain\n";
66 return temp;
67}
68
69//__________________
70// Select whether to read constrained or not constrained momentum
71void AliFemtoEventReaderESDChain::SetConstrained(const bool constrained)
72{
73 fConstrained=constrained;
74}
75//__________________
76// Check whether we read constrained or not constrained momentum
77bool AliFemtoEventReaderESDChain::GetConstrained() const
78{
79 return fConstrained;
80}
81//__________________
82AliFemtoEvent* AliFemtoEventReaderESDChain::ReturnHbtEvent()
83{
84 // Get the event, read all the relevant information
85 // and fill the AliFemtoEvent class
86 // Returns a valid AliFemtoEvent
87 AliFemtoEvent *hbtEvent = 0;
88 string tFriendFileName;
89
90 // Get the friend information
91 cout<<"starting to read event "<<fCurEvent<<endl;
92 fEvent->SetESDfriend(fEventFriend);
93 vector<int> label_table;//to check labels
94
95 hbtEvent = new AliFemtoEvent;
96 //setting basic things
b2f60a91 97 //hbtEvent->SetEventNumber(fEvent->GetEventNumber());
67427ff7 98 hbtEvent->SetRunNumber(fEvent->GetRunNumber());
99 //hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
100 hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
101 hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
102 hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
103 hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
104 hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
105 hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
106 hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
107 hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
108 hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
109
110 //Vertex
111 double fV1[3];
112 fEvent->GetVertex()->GetXYZ(fV1);
113
114 AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
115 hbtEvent->SetPrimVertPos(vertex);
116
117 //starting to reading tracks
118 int nofTracks=0; //number of reconstructed tracks in event
119 nofTracks=fEvent->GetNumberOfTracks();
120 int realnofTracks=0;//number of track which we use ina analysis
121
122 // Clear the shared cluster list
123 for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
124 fClusterPerPadrow[tPad]->clear();
125 }
126 for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
127 fSharedList[tPad]->clear();
128 }
129
130
131 for (int i=0;i<nofTracks;i++) {
132 const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
133
134 list<Int_t>::iterator tClustIter;
135
136 Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster];
137 Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices);
138 for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
139 if (tTrackIndices[tNcl] >= 0) {
140 tClustIter = find(fClusterPerPadrow[tNcl]->begin(), fClusterPerPadrow[tNcl]->end(), tTrackIndices[tNcl]);
141 if (tClustIter == fClusterPerPadrow[tNcl]->end()) {
142 fClusterPerPadrow[tNcl]->push_back(tTrackIndices[tNcl]);
143 }
144 else {
145 fSharedList[tNcl]->push_back(tTrackIndices[tNcl]);
146 }
147 }
148 }
149
150 }
151
152 for (int i=0;i<nofTracks;i++)
153 {
154 bool good_momentum=true; //flaga to chcek if we can read momentum of this track
155
156 AliFemtoTrack* trackCopy = new AliFemtoTrack();
157 const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
158 const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
159
160 trackCopy->SetCharge((short)esdtrack->GetSign());
161
162 //in aliroot we have AliPID
163 //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
164 //we use only 5 first
165 double esdpid[5];
166 esdtrack->GetESDpid(esdpid);
167 trackCopy->SetPidProbElectron(esdpid[0]);
168 trackCopy->SetPidProbMuon(esdpid[1]);
169 trackCopy->SetPidProbPion(esdpid[2]);
170 trackCopy->SetPidProbKaon(esdpid[3]);
171 trackCopy->SetPidProbProton(esdpid[4]);
172
173 double pxyz[3];
174 if (fConstrained==true)
175 good_momentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
176 else
177 good_momentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
178
179 AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
180 if (v.mag() == 0) {
181 // cout << "Found 0 momentum ???? " <<endl;
182 delete trackCopy;
183 continue;
184 }
185 trackCopy->SetP(v);//setting momentum
186 trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
187 const AliFmThreeVectorD p(pxyz[0],pxyz[1],pxyz[2]);
188 const AliFmThreeVectorD origin(fV1[0],fV1[1],fV1[2]);
189 //setting helix I do not if it is ok
190 AliFmPhysicalHelixD helix(p,origin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
191 trackCopy->SetHelix(helix);
192
193 trackCopy->SetTrackId(esdtrack->GetID());
194 trackCopy->SetFlags(esdtrack->GetStatus());
195 //trackCopy->SetLabel(esdtrack->GetLabel());
196
197 //some stuff which could be useful
198 float impact[2];
199 float covimpact[3];
200 esdtrack->GetImpactParameters(impact,covimpact);
201 trackCopy->SetImpactD(impact[0]);
202 trackCopy->SetImpactZ(impact[1]);
203 trackCopy->SetCdd(covimpact[0]);
204 trackCopy->SetCdz(covimpact[1]);
205 trackCopy->SetCzz(covimpact[2]);
206 trackCopy->SetITSchi2(esdtrack->GetITSchi2());
207 trackCopy->SetITSncls(esdtrack->GetNcls(0));
208 trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());
209 trackCopy->SetTPCncls(esdtrack->GetTPCNcls());
210 trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());
211 trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class
212 trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma());
213
214 // Fill cluster per padrow information
215 Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster];
216 Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices);
217 for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
218 if (tTrackIndices[tNcl] > 0)
219 trackCopy->SetTPCcluster(tNcl, 1);
220 else
221 trackCopy->SetTPCcluster(tNcl, 0);
222 }
223
224 // Fill shared cluster information
225 list<Int_t>::iterator tClustIter;
226
227 for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
228 if (tTrackIndices[tNcl] > 0) {
229 tClustIter = find(fSharedList[tNcl]->begin(), fSharedList[tNcl]->end(), tTrackIndices[tNcl]);
230 if (tClustIter != fSharedList[tNcl]->end()) {
231 trackCopy->SetTPCshared(tNcl, 1);
232 cout << "Event next" << endl;
233 cout << "Track: " << i << endl;
234 cout << "Shared cluster: " << tNcl << " " << tTrackIndices[tNcl] << endl;
235 }
236 else {
237 trackCopy->SetTPCshared(tNcl, 0);
238 }
239 }
240 }
241
242 if (good_momentum==true)
243 {
244 hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
245 realnofTracks++;//real number of tracks
246 // delete trackCopy;
247 }
248 else
249 {
250 delete trackCopy;
251 }
252
253 }
254
255 hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
256 fCurEvent++;
257 cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
258 return hbtEvent;
259}
260//___________________
261// The chain loads the ESD for us
262// You must provide the address where it can be found
263void AliFemtoEventReaderESDChain::SetESDSource(AliESD *aESD)
264{
265 fEvent = aESD;
266}
267//___________________
268// We need the ESD tree to obtain
269// information about the friend file location
270void AliFemtoEventReaderESDChain::SetESDfriendSource(AliESDfriend *aFriend)
271{
272 fEventFriend = aFriend;
273}
274
275
276
277
278
279
280
281
282