Modifications needed for the compilation of the femtoscopy code (Adam-Mike)
[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 {
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
45 AliFemtoEventReaderESDChain::~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
63 AliFemtoString 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
71 void AliFemtoEventReaderESDChain::SetConstrained(const bool constrained)
72 {
73   fConstrained=constrained;
74 }
75 //__________________
76 // Check whether we read constrained or not constrained momentum
77 bool AliFemtoEventReaderESDChain::GetConstrained() const
78 {
79   return fConstrained;
80 }
81 //__________________
82 AliFemtoEvent* 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
97   //hbtEvent->SetEventNumber(fEvent->GetEventNumber());
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
263 void 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
270 void AliFemtoEventReaderESDChain::SetESDfriendSource(AliESDfriend *aFriend)
271 {
272   fEventFriend = aFriend;
273 }
274
275
276
277
278
279
280
281
282