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