fdeea3d33c401a05d6216e8d237860df0aef5a7b
[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemto / Reader / AliFemtoEventReaderESD.cxx
1 /*
2  *$Id$
3  *$Log$
4  *Revision 1.5  2007-04-03 16:00:08  mchojnacki
5  *Changes to iprove memory managing
6  *
7  *Revision 1.4  2007/03/13 15:30:03  mchojnacki
8  *adding reader for simulated data
9  *
10  *Revision 1.3  2007/03/08 14:58:03  mchojnacki
11  *adding some alice stuff
12  *
13  *Revision 1.2  2007/03/07 13:36:16  mchojnacki
14  *Add some comments
15  *
16  *Revision 1.1.1.1  2007/03/07 10:14:49  mchojnacki
17  *First version on CVS
18  *
19  */
20
21 #include "AliFemtoEventReaderESD.h"
22
23 #include "TFile.h"
24 #include "TTree.h"
25 #include "AliESD.h"
26 #include "AliESDtrack.h"
27
28 //#include "TSystem.h"
29
30 #include "AliFmPhysicalHelixD.h"
31 #include "AliFmThreeVectorF.h"
32
33 #include "SystemOfUnits.h"
34
35 #include "AliFemtoEvent.h"
36
37 ClassImp(AliFemtoEventReaderESD)
38
39 #if !(ST_NO_NAMESPACES)
40   using namespace units;
41 #endif
42
43 using namespace std;
44 //____________________________
45 //constructor with 0 parameters , look at default settings 
46 AliFemtoEventReaderESD::AliFemtoEventReaderESD():
47   fInputFile(" "),
48   fFileName(" "),
49   fConstrained(true),
50   fNumberofEvent(0),
51   fCurEvent(0),
52   fCurFile(0),
53   fTree(0x0),
54   fEvent(0x0),
55   fEsdFile(0x0),
56   fEventFriend(0)
57 {
58   fClusterPerPadrow = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
59   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
60     fClusterPerPadrow[tPad] = new list<Int_t>();
61   }
62   fSharedList = (list<Int_t> **) malloc(sizeof(list<Int_t> *) * AliESDfriendTrack::kMaxTPCcluster);
63   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
64     fSharedList[tPad] = new list<Int_t>();
65   }
66 }
67
68 //__________________
69 //Destructor
70 AliFemtoEventReaderESD::~AliFemtoEventReaderESD()
71 {
72   //delete fListOfFiles;
73   delete fTree;
74   delete fEvent;
75   delete fEsdFile;
76
77   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
78     fClusterPerPadrow[tPad]->clear();
79     delete fClusterPerPadrow[tPad];
80   }
81   delete [] fClusterPerPadrow;
82   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
83     fSharedList[tPad]->clear();
84     delete fSharedList[tPad];
85   }
86   delete [] fSharedList;
87 }
88
89 //__________________
90 AliFemtoString AliFemtoEventReaderESD::Report()
91 {
92   AliFemtoString temp = "\n This is the AliFemtoEventReaderESD\n";
93   return temp;
94 }
95
96 //__________________
97 //setting the name of file where names of ESD file are written 
98 //it takes only this files which have good trees
99 void AliFemtoEventReaderESD::SetInputFile(const char* inputFile)
100 {
101   char buffer[256];
102   fInputFile=string(inputFile);
103   cout<<"Input File set on "<<fInputFile<<endl;
104   ifstream infile(inputFile);
105   if(infile.good()==true)
106     { 
107       //checking if all give files have good tree inside
108       while (infile.eof()==false)
109         {
110           infile.getline(buffer,256);
111           //ifstream test_file(buffer);
112           TFile *esdFile=TFile::Open(buffer,"READ");
113           if (esdFile!=0x0)
114             {   
115               TTree* tree = (TTree*) esdFile->Get("esdTree");
116               if (tree!=0x0)
117                 {
118                   cout<<"putting file  "<<string(buffer)<<" into analysis"<<endl;
119                   fListOfFiles.push_back(string(buffer));
120                   delete tree;
121                 }
122               esdFile->Close(); 
123             }
124           delete esdFile;
125         }
126     }
127 }
128
129 //setting the next file to read 
130 bool AliFemtoEventReaderESD::GetNextFile()
131 {       
132   if (fCurFile>=fListOfFiles.size())
133     return false;
134   fFileName=fListOfFiles.at(fCurFile);  
135   cout<<"FileName set on "<<fFileName<<" "<<fCurFile<<endl;
136
137   fCurFile++;
138   return true;
139 }
140 void AliFemtoEventReaderESD::SetConstrained(const bool constrained)
141 {
142   fConstrained=constrained;
143 }
144
145 bool AliFemtoEventReaderESD::GetConstrained() const
146 {
147   return fConstrained;
148 }
149
150 AliFemtoEvent* AliFemtoEventReaderESD::ReturnHbtEvent()
151 {
152   AliFemtoEvent *hbtEvent = 0;
153   string tFriendFileName;
154
155   if (fCurEvent==fNumberofEvent)//open next file  
156     {
157       cout<<"next file"<<endl;
158       if(GetNextFile()) 
159         {
160           delete fEventFriend;
161           fEventFriend = 0;
162           delete fEvent;//added 1.04.2007
163           fEvent=new AliESD();
164           delete fTree;
165           //fTree=0;
166           delete fEsdFile;
167                 
168           //ESD data
169           fEsdFile=TFile::Open(fFileName.c_str(),"READ");
170           fTree = (TTree*) fEsdFile->Get("esdTree");                    
171           fTree->SetBranchAddress("ESD", &fEvent);                      
172
173           // Attach the friend tree with additional information
174           tFriendFileName = fFileName;
175           tFriendFileName.insert(tFriendFileName.find("s.root"),"friend");
176           cout << "Reading friend " << tFriendFileName.c_str() << endl;;
177           fTree->AddFriend("esdFriendTree",tFriendFileName.c_str());
178           fTree->SetBranchAddress("ESDfriend",&fEventFriend);
179
180           fNumberofEvent=fTree->GetEntries();
181           cout<<"Number of Entries in file "<<fNumberofEvent<<endl;
182           fCurEvent=0;
183           //sim data
184         }
185       else //no more data to read
186         {
187           cout<<"no more files "<<hbtEvent<<endl;
188           fReaderStatus=1;
189           return hbtEvent; 
190         }
191     }           
192   cout<<"starting to read event "<<fCurEvent<<endl;
193   fTree->GetEvent(fCurEvent);//getting next event
194   fEvent->SetESDfriend(fEventFriend);
195   vector<int> label_table;//to check labels
196         
197   hbtEvent = new AliFemtoEvent;
198   //setting basic things
199   hbtEvent->SetEventNumber(fEvent->GetEventNumber());
200   hbtEvent->SetRunNumber(fEvent->GetRunNumber());
201   //hbtEvent->SetNumberOfTracks(fEvent->GetNumberOfTracks());
202   hbtEvent->SetMagneticField(fEvent->GetMagneticField()*kilogauss);//to check if here is ok
203   hbtEvent->SetZDCN1Energy(fEvent->GetZDCN1Energy());
204   hbtEvent->SetZDCP1Energy(fEvent->GetZDCP1Energy());
205   hbtEvent->SetZDCN2Energy(fEvent->GetZDCN2Energy());
206   hbtEvent->SetZDCP2Energy(fEvent->GetZDCP2Energy());
207   hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
208   hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
209   hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
210   hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
211         
212   //Vertex
213   double fV1[3];
214   fEvent->GetVertex()->GetXYZ(fV1);
215
216   AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
217   hbtEvent->SetPrimVertPos(vertex);
218         
219   //starting to reading tracks
220   int nofTracks=0;  //number of reconstructed tracks in event
221   nofTracks=fEvent->GetNumberOfTracks();
222   int realnofTracks=0;//number of track which we use ina analysis
223
224   // Clear the shared cluster list
225   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
226     fClusterPerPadrow[tPad]->clear();
227   }
228   for (int tPad=0; tPad<AliESDfriendTrack::kMaxTPCcluster; tPad++) {
229     fSharedList[tPad]->clear();
230   }
231
232
233   for (int i=0;i<nofTracks;i++) {
234     const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
235
236     list<Int_t>::iterator tClustIter;
237
238     Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster];
239     Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices);
240     for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
241       if (tTrackIndices[tNcl] >= 0) {
242         tClustIter = find(fClusterPerPadrow[tNcl]->begin(), fClusterPerPadrow[tNcl]->end(), tTrackIndices[tNcl]);
243         if (tClustIter == fClusterPerPadrow[tNcl]->end()) {
244           fClusterPerPadrow[tNcl]->push_back(tTrackIndices[tNcl]);
245         }
246         else {
247           fSharedList[tNcl]->push_back(tTrackIndices[tNcl]);
248         }
249       }
250     }
251       
252   }
253
254   for (int i=0;i<nofTracks;i++)
255     {
256       bool  good_momentum=true; //flaga to chcek if we can read momentum of this track
257                 
258       AliFemtoTrack* trackCopy = new AliFemtoTrack();   
259       const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
260       const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
261
262       trackCopy->SetCharge((short)esdtrack->GetSign());
263
264       //in aliroot we have AliPID 
265       //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon   
266       //we use only 5 first
267       double esdpid[5];
268       esdtrack->GetESDpid(esdpid);
269       trackCopy->SetPidProbElectron(esdpid[0]);
270       trackCopy->SetPidProbMuon(esdpid[1]);
271       trackCopy->SetPidProbPion(esdpid[2]);
272       trackCopy->SetPidProbKaon(esdpid[3]);
273       trackCopy->SetPidProbProton(esdpid[4]);
274                                                 
275       double pxyz[3];
276       if (fConstrained==true)               
277         good_momentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
278       else
279         good_momentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
280       AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
281       trackCopy->SetP(v);//setting momentum
282       trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
283       const AliFmThreeVectorD p(pxyz[0],pxyz[1],pxyz[2]);
284       if (p.mag() == 0) {
285         delete trackCopy;
286         continue;
287       }
288       const AliFmThreeVectorD origin(fV1[0],fV1[1],fV1[2]);
289       //setting helix I do not if it is ok
290       AliFmPhysicalHelixD helix(p,origin,(double)(fEvent->GetMagneticField())*kilogauss,(double)(trackCopy->Charge())); 
291       trackCopy->SetHelix(helix);
292                 
293       trackCopy->SetTrackId(esdtrack->GetID());
294       trackCopy->SetFlags(esdtrack->GetStatus());
295       //trackCopy->SetLabel(esdtrack->GetLabel());
296                 
297       //some stuff which could be useful 
298       float impact[2];
299       float covimpact[3];
300       esdtrack->GetImpactParameters(impact,covimpact);
301       trackCopy->SetImpactD(impact[0]);
302       trackCopy->SetImpactZ(impact[1]);
303       trackCopy->SetCdd(covimpact[0]);
304       trackCopy->SetCdz(covimpact[1]);
305       trackCopy->SetCzz(covimpact[2]);
306       trackCopy->SetITSchi2(esdtrack->GetITSchi2());    
307       trackCopy->SetITSncls(esdtrack->GetNcls(0));     
308       trackCopy->SetTPCchi2(esdtrack->GetTPCchi2());       
309       trackCopy->SetTPCncls(esdtrack->GetTPCNcls());       
310       trackCopy->SetTPCnclsF(esdtrack->GetTPCNclsF());      
311       trackCopy->SetTPCsignalN((short)esdtrack->GetTPCsignalN()); //due to bug in aliesdtrack class   
312       trackCopy->SetTPCsignalS(esdtrack->GetTPCsignalSigma()); 
313
314       // Fill cluster per padrow information
315       Int_t tTrackIndices[AliESDfriendTrack::kMaxTPCcluster];
316       Int_t tNClusters = esdtrack->GetTPCclusters(tTrackIndices);
317       for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
318         if (tTrackIndices[tNcl] > 0)
319           trackCopy->SetTPCcluster(tNcl, 1);
320         else
321           trackCopy->SetTPCcluster(tNcl, 0);
322       }
323       
324       // Fill shared cluster information
325       list<Int_t>::iterator tClustIter;
326
327       for (int tNcl=0; tNcl<AliESDfriendTrack::kMaxTPCcluster; tNcl++) {
328         if (tTrackIndices[tNcl] > 0) {
329           tClustIter = find(fSharedList[tNcl]->begin(), fSharedList[tNcl]->end(), tTrackIndices[tNcl]);
330           if (tClustIter != fSharedList[tNcl]->end()) {
331             trackCopy->SetTPCshared(tNcl, 1);
332             cout << "Event next" <<  endl;
333             cout << "Track: " << i << endl;
334             cout << "Shared cluster: " << tNcl << " " << tTrackIndices[tNcl] << endl;
335           }
336           else {
337             trackCopy->SetTPCshared(tNcl, 0);
338           }
339         }
340       }
341
342       //decision if we want this track
343       //if we using diffrent labels we want that this label was use for first time 
344       //if we use hidden info we want to have match between sim data and ESD
345       if (good_momentum==true)
346         {
347           hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
348           realnofTracks++;//real number of tracks
349         }
350       else
351         {
352           delete  trackCopy;
353         }
354                 
355     }
356
357   hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event   
358   fCurEvent++;  
359   cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
360   if (fCurEvent== fNumberofEvent)//if end of current file close all
361     {   
362       fTree->Reset(); 
363       delete fTree;
364       fEsdFile->Close();
365     }
366   return hbtEvent; 
367 }
368
369
370
371
372
373
374
375
376