be4b8c94b439757686c5410f37934d1f41d4d84d
[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   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 }
81 //__________________
82 //Destructor
83 AliFemtoEventReaderESDChain::~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
99 //__________________
100 // Assignment operator
101 AliFemtoEventReaderESDChain& 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 }
150 //__________________
151 // Simple report
152 AliFemtoString 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
160 void AliFemtoEventReaderESDChain::SetConstrained(const bool constrained)
161 {
162   fConstrained=constrained;
163 }
164 //__________________
165 // Check whether we read constrained or not constrained momentum
166 bool AliFemtoEventReaderESDChain::GetConstrained() const
167 {
168   return fConstrained;
169 }
170 //__________________
171 AliFemtoEvent* 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
186   //  hbtEvent->SetEventNumber(fEvent->GetEventNumber());
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
247       //      const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
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
352 void 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
359 void AliFemtoEventReaderESDChain::SetESDfriendSource(AliESDfriend *aFriend)
360 {
361   fEventFriend = aFriend;
362 }
363
364
365
366
367
368
369
370
371