Initial check-in of the model classes
[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),
0215f606 31 fEventFriend(0),
32 fSharedList(0x0),
33 fClusterPerPadrow(0x0)
67427ff7 34{
35 fClusterPerPadrow = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
36 for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
37 fClusterPerPadrow[tPad] = new list<Int_t>();
38 }
39 fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
40 for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
41 fSharedList[tPad] = new list<Int_t>();
42 }
43}
44
0215f606 45//__________________
46// Copy constructor
47AliFemtoEventReaderESDChain::AliFemtoEventReaderESDChain(const AliFemtoEventReaderESDChain& aReader):
48 fFileName(" "),
49 fConstrained(true),
50 fNumberofEvent(0),
51 fCurEvent(0),
52 fCurFile(0),
53 fEvent(0x0),
54 fEventFriend(0),
55 fSharedList(0x0),
56 fClusterPerPadrow(0x0)
57{
58 fConstrained = aReader.fConstrained;
59 fNumberofEvent = aReader.fNumberofEvent;
60 fCurEvent = aReader.fCurEvent;
61 fCurFile = aReader.fCurFile;
62 fEvent = new AliESD(*aReader.fEvent);
63 fEventFriend = aReader.fEventFriend;
64 fClusterPerPadrow = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
65 for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
66 fClusterPerPadrow[tPad] = new list<Int_t>();
67 list<Int_t>::iterator iter;
68 for (iter=aReader.fClusterPerPadrow[tPad]->begin(); iter!=aReader.fClusterPerPadrow[tPad]->end(); iter++) {
69 fClusterPerPadrow[tPad]->push_back(*iter);
70 }
71 }
72 fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
73 for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
74 fSharedList[tPad] = new list<Int_t>();
75 list<Int_t>::iterator iter;
76 for (iter=aReader.fSharedList[tPad]->begin(); iter!=aReader.fSharedList[tPad]->end(); iter++) {
77 fSharedList[tPad]->push_back(*iter);
78 }
79 }
80}
67427ff7 81//__________________
82//Destructor
83AliFemtoEventReaderESDChain::~AliFemtoEventReaderESDChain()
84{
85 delete fEvent;
86
87 for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
88 fClusterPerPadrow[tPad]->clear();
89 delete fClusterPerPadrow[tPad];
90 }
91 delete [] fClusterPerPadrow;
92 for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
93 fSharedList[tPad]->clear();
94 delete fSharedList[tPad];
95 }
96 delete [] fSharedList;
97}
98
0215f606 99//__________________
100// Assignment operator
101AliFemtoEventReaderESDChain& AliFemtoEventReaderESDChain::operator=(const AliFemtoEventReaderESDChain& aReader)
102{
103 if (this == &aReader)
104 return *this;
105
106 fConstrained = aReader.fConstrained;
107 fNumberofEvent = aReader.fNumberofEvent;
108 fCurEvent = aReader.fCurEvent;
109 fCurFile = aReader.fCurFile;
110 if (fEvent) delete fEvent;
111 fEvent = new AliESD(*aReader.fEvent);
112
113 fEventFriend = aReader.fEventFriend;
114
115 if (fClusterPerPadrow) {
116 for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
117 fClusterPerPadrow[tPad]->clear();
118 delete fClusterPerPadrow[tPad];
119 }
120 delete [] fClusterPerPadrow;
121 }
122
123 if (fSharedList) {
124 for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
125 fSharedList[tPad]->clear();
126 delete fSharedList[tPad];
127 }
128 delete [] fSharedList;
129 }
130
131 fClusterPerPadrow = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
132 for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
133 fClusterPerPadrow[tPad] = new list<Int_t>();
134 list<Int_t>::iterator iter;
135 for (iter=aReader.fClusterPerPadrow[tPad]->begin(); iter!=aReader.fClusterPerPadrow[tPad]->end(); iter++) {
136 fClusterPerPadrow[tPad]->push_back(*iter);
137 }
138 }
139 fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
140 for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
141 fSharedList[tPad] = new list<Int_t>();
142 list<Int_t>::iterator iter;
143 for (iter=aReader.fSharedList[tPad]->begin(); iter!=aReader.fSharedList[tPad]->end(); iter++) {
144 fSharedList[tPad]->push_back(*iter);
145 }
146 }
147
148 return *this;
149}
67427ff7 150//__________________
151// Simple report
152AliFemtoString AliFemtoEventReaderESDChain::Report()
153{
154 AliFemtoString temp = "\n This is the AliFemtoEventReaderESDChain\n";
155 return temp;
156}
157
158//__________________
159// Select whether to read constrained or not constrained momentum
160void AliFemtoEventReaderESDChain::SetConstrained(const bool constrained)
161{
162 fConstrained=constrained;
163}
164//__________________
165// Check whether we read constrained or not constrained momentum
166bool AliFemtoEventReaderESDChain::GetConstrained() const
167{
168 return fConstrained;
169}
170//__________________
171AliFemtoEvent* AliFemtoEventReaderESDChain::ReturnHbtEvent()
172{
173 // Get the event, read all the relevant information
174 // and fill the AliFemtoEvent class
175 // Returns a valid AliFemtoEvent
176 AliFemtoEvent *hbtEvent = 0;
177 string tFriendFileName;
178
179 // Get the friend information
180 cout<<"starting to read event "<<fCurEvent<<endl;
181 fEvent->SetESDfriend(fEventFriend);
182 vector<int> label_table;//to check labels
183
184 hbtEvent = new AliFemtoEvent;
185 //setting basic things
a531c6fc 186 // hbtEvent->SetEventNumber(fEvent->GetEventNumber());
67427ff7 187 hbtEvent->SetRunNumber(fEvent->GetRunNumber());
188 //hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
189 hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
190 hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
191 hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
192 hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
193 hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
194 hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
195 hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
196 hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
197 hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
198
199 //Vertex
200 double fV1[3];
201 fEvent->GetVertex()->GetXYZ(fV1);
202
203 AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
204 hbtEvent->SetPrimVertPos(vertex);
205
206 //starting to reading tracks
207 int nofTracks=0; //number of reconstructed tracks in event
208 nofTracks=fEvent->GetNumberOfTracks();
209 int realnofTracks=0;//number of track which we use ina analysis
210
211 // Clear the shared cluster list
212 for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
213 fClusterPerPadrow[tPad]->clear();
214 }
215 for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
216 fSharedList[tPad]->clear();
217 }
218
219
220 for (int i=0;i<nofTracks;i++) {
221 const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
222
223 list<Int_t>::iterator tClustIter;
224
225 Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster];
226 Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices);
227 for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
228 if (tTrackIndices[tNcl] >= 0) {
229 tClustIter = find(fClusterPerPadrow[tNcl]->begin(), fClusterPerPadrow[tNcl]->end(), tTrackIndices[tNcl]);
230 if (tClustIter == fClusterPerPadrow[tNcl]->end()) {
231 fClusterPerPadrow[tNcl]->push_back(tTrackIndices[tNcl]);
232 }
233 else {
234 fSharedList[tNcl]->push_back(tTrackIndices[tNcl]);
235 }
236 }
237 }
238
239 }
240
241 for (int i=0;i<nofTracks;i++)
242 {
243 bool good_momentum=true; //flaga to chcek if we can read momentum of this track
244
245 AliFemtoTrack* trackCopy = new AliFemtoTrack();
246 const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
0215f606 247 // const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
67427ff7 248
249 trackCopy->SetCharge((short)esdtrack->GetSign());
250
251 //in aliroot we have AliPID
252 //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon
253 //we use only 5 first
254 double esdpid[5];
255 esdtrack->GetESDpid(esdpid);
256 trackCopy->SetPidProbElectron(esdpid[0]);
257 trackCopy->SetPidProbMuon(esdpid[1]);
258 trackCopy->SetPidProbPion(esdpid[2]);
259 trackCopy->SetPidProbKaon(esdpid[3]);
260 trackCopy->SetPidProbProton(esdpid[4]);
261
262 double pxyz[3];
263 if (fConstrained==true)
264 good_momentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
265 else
266 good_momentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
267
268 AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
269 if (v.mag() == 0) {
270 // cout << "Found 0 momentum ???? " <<endl;
271 delete trackCopy;
272 continue;
273 }
274 trackCopy->SetP(v);//setting momentum
275 trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
276 const AliFmThreeVectorD p(pxyz[0],pxyz[1],pxyz[2]);
277 const AliFmThreeVectorD origin(fV1[0],fV1[1],fV1[2]);
278 //setting helix I do not if it is ok
279 AliFmPhysicalHelixD helix(p,origin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge()));
280 trackCopy->SetHelix(helix);
281
282 trackCopy->SetTrackId(esdtrack->GetID());
283 trackCopy->SetFlags(esdtrack->GetStatus());
284 //trackCopy->SetLabel(esdtrack->GetLabel());
285
286 //some stuff which could be useful
287 float impact[2];
288 float covimpact[3];
289 esdtrack->GetImpactParameters(impact,covimpact);
290 trackCopy->SetImpactD(impact[0]);
291 trackCopy->SetImpactZ(impact[1]);
292 trackCopy->SetCdd(covimpact[0]);
293 trackCopy->SetCdz(covimpact[1]);
294 trackCopy->SetCzz(covimpact[2]);
295 trackCopy->SetITSchi2(esdtrack->GetITSchi2());
296 trackCopy->SetITSncls(esdtrack->GetNcls(0));
297 trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());
298 trackCopy->SetTPCncls(esdtrack->GetTPCNcls());
299 trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());
300 trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class
301 trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma());
302
303 // Fill cluster per padrow information
304 Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster];
305 Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices);
306 for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
307 if (tTrackIndices[tNcl] > 0)
308 trackCopy->SetTPCcluster(tNcl, 1);
309 else
310 trackCopy->SetTPCcluster(tNcl, 0);
311 }
312
313 // Fill shared cluster information
314 list<Int_t>::iterator tClustIter;
315
316 for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
317 if (tTrackIndices[tNcl] > 0) {
318 tClustIter = find(fSharedList[tNcl]->begin(), fSharedList[tNcl]->end(), tTrackIndices[tNcl]);
319 if (tClustIter != fSharedList[tNcl]->end()) {
320 trackCopy->SetTPCshared(tNcl, 1);
321 cout << "Event next" << endl;
322 cout << "Track: " << i << endl;
323 cout << "Shared cluster: " << tNcl << " " << tTrackIndices[tNcl] << endl;
324 }
325 else {
326 trackCopy->SetTPCshared(tNcl, 0);
327 }
328 }
329 }
330
331 if (good_momentum==true)
332 {
333 hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
334 realnofTracks++;//real number of tracks
335 // delete trackCopy;
336 }
337 else
338 {
339 delete trackCopy;
340 }
341
342 }
343
344 hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
345 fCurEvent++;
346 cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
347 return hbtEvent;
348}
349//___________________
350// The chain loads the ESD for us
351// You must provide the address where it can be found
352void AliFemtoEventReaderESDChain::SetESDSource(AliESD *aESD)
353{
354 fEvent = aESD;
355}
356//___________________
357// We need the ESD tree to obtain
358// information about the friend file location
359void AliFemtoEventReaderESDChain::SetESDfriendSource(AliESDfriend *aFriend)
360{
361 fEventFriend = aFriend;
362}
363
364
365
366
367
368
369
370
371