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