903bfc615b4347a6159eae2c347ab7fed7eef6b1
[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemto / Reader / AliFemtoEventReaderESDChain.cxx
1 #include "AliFemtoEventReaderESDChain.h"
2
3 #include "TFile.h"
4 #include "TTree.h"
5 #include "AliESD.h"
6 #include "AliESDtrack.h"
7
8 #include "Infrastructure/AliFmPhysicalHelixD.h"
9 #include "Infrastructure/AliFmThreeVectorF.h"
10
11 #include "Base/SystemOfUnits.h"
12
13 #include "Infrastructure/AliFemtoEvent.h"
14
15 ClassImp(AliFemtoEventReaderESDChain)
16
17 #if !(ST_NO_NAMESPACES)
18   using namespace units;
19 #endif
20
21 using namespace std;
22 //____________________________
23 //constructor with 0 parameters , look at default settings 
24 AliFemtoEventReaderESDChain::AliFemtoEventReaderESDChain():
25   fFileName(" "),
26   fConstrained(true),
27   fNumberofEvent(0),
28   fCurEvent(0),
29   fCurFile(0),
30   fEvent(0x0),
31   fEventFriend(0),
32   fSharedList(0x0),
33   fClusterPerPadrow(0x0)
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
45 //__________________
46 // Copy constructor
47 AliFemtoEventReaderESDChain::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   fEvent = new AliESD();
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 }
82 //__________________
83 //Destructor
84 AliFemtoEventReaderESDChain::~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
100 //__________________
101 // Assignment operator
102 AliFemtoEventReaderESDChain& 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;
112   //fEvent = new AliESD(*aReader.fEvent);
113   fEvent = new AliESD();
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 }
152 //__________________
153 // Simple report
154 AliFemtoString 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
162 void AliFemtoEventReaderESDChain::SetConstrained(const bool constrained)
163 {
164   fConstrained=constrained;
165 }
166 //__________________
167 // Check whether we read constrained or not constrained momentum
168 bool AliFemtoEventReaderESDChain::GetConstrained() const
169 {
170   return fConstrained;
171 }
172 //__________________
173 AliFemtoEvent* 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
188   //  hbtEvent->SetEventNumber(fEvent->GetEventNumber());
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
249       //      const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
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
354 void 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
361 void AliFemtoEventReaderESDChain::SetESDfriendSource(AliESDfriend *aFriend)
362 {
363   fEventFriend = aFriend;
364 }
365
366
367
368
369
370
371
372
373