Bug fix
[u/mrichter/AliRoot.git] / HBTAN / AliHBTReaderITSv1.cxx
1
2 #include "AliHBTReaderITSv1.h"
3 #include "AliHBTEvent.h"
4 #include "AliHBTRun.h"
5 #include "AliHBTParticle.h"
6 #include "AliHBTParticleCut.h"
7
8 #include <iostream.h>
9
10 #include <TROOT.h>
11 #include <TFile.h>
12 #include <TTree.h>
13 #include <TBranch.h>
14 #include <TObjArray.h>
15 #include <TParticle.h>
16 #include <TString.h>
17 #include <TObjString.h>
18
19 #include <AliRun.h>
20 #include <AliMagF.h>
21 #include <AliKalmanTrack.h>
22 #include <AliITSIOTrack.h>
23
24 ClassImp(AliHBTReaderITSv1)
25 /********************************************************************/
26
27 AliHBTReaderITSv1::
28 AliHBTReaderITSv1(const Char_t* tracksfilename,const Char_t* galicefilename):
29                  fITSTracksFileName(tracksfilename),fGAliceFileName(galicefilename)
30  {
31      fParticles = new AliHBTRun();
32      fTracks    = new AliHBTRun();
33      fIsRead = kFALSE;
34  }
35 /********************************************************************/
36
37 AliHBTReaderITSv1::
38 AliHBTReaderITSv1(TObjArray* dirs, const Char_t* tracksfilename,const Char_t* galicefilename):
39                  AliHBTReader(dirs),
40                  fITSTracksFileName(tracksfilename),fGAliceFileName(galicefilename)
41  {
42    fParticles = new AliHBTRun();
43    fTracks    = new AliHBTRun();
44    fIsRead    = kFALSE;
45  }
46 /********************************************************************/
47
48 AliHBTReaderITSv1::~AliHBTReaderITSv1()
49 {
50    delete fParticles;
51    delete fTracks;
52 }
53 /********************************************************************/
54
55 AliHBTEvent* AliHBTReaderITSv1::GetParticleEvent(Int_t n)
56  {
57  //returns Nth event with simulated particles
58    if (!fIsRead) 
59     if(Read(fParticles,fTracks))
60      {
61        Error("GetParticleEvent","Error in reading");
62        return 0x0;
63      }
64
65    return fParticles->GetEvent(n);
66  }
67 /********************************************************************/
68
69 AliHBTEvent* AliHBTReaderITSv1::GetTrackEvent(Int_t n)
70  {
71  //returns Nth event with reconstructed tracks
72    if (!fIsRead) 
73     if(Read(fParticles,fTracks))
74      {
75        Error("GetTrackEvent","Error in reading");
76        return 0x0;
77      }
78    return fTracks->GetEvent(n);
79  }
80 /********************************************************************/
81
82 Int_t AliHBTReaderITSv1::GetNumberOfPartEvents()
83  {
84  //returns number of events of particles
85    if (!fIsRead)
86     if(Read(fParticles,fTracks))
87      {
88        Error("GetNumberOfPartEvents","Error in reading");
89        return 0;
90      }
91    return fParticles->GetNumberOfEvents();
92  }
93
94 /********************************************************************/
95 Int_t AliHBTReaderITSv1::GetNumberOfTrackEvents()
96  {
97  //returns number of events of tracks
98   if (!fIsRead) 
99     if(Read(fParticles,fTracks))
100      {
101        Error("GetNumberOfTrackEvents","Error in reading");
102        return 0;
103      }
104   return fTracks->GetNumberOfEvents();
105  }
106 /********************************************************************/
107
108 Int_t AliHBTReaderITSv1::Read(AliHBTRun* particles, AliHBTRun *tracks)
109 {
110  cout<<"AliHBTReaderITSv1::Read()"<<endl;
111  Int_t Nevents = 0;
112  AliITSIOTrack *iotrack=new AliITSIOTrack;
113  Int_t currentdir = 0;
114  Int_t Ndirs;
115  Int_t totalNevents = 0;
116  
117  if (fDirs)
118   {
119     Ndirs = fDirs->GetEntries();
120   }
121  else
122   {
123     Ndirs = 0;
124   }
125  
126  do //do while is good even if 
127   {  
128    TFile* gAliceFile = OpenGAliceFile(currentdir);
129    if(gAliceFile == 0x0)
130     {
131        Error("Read","Can not open the file with gAlice");
132        delete iotrack;
133        return 1;
134     }
135    TFile *file = OpenTrackFile(currentdir);
136    if(file == 0x0)
137     {
138        Error("Read","Can not open the file with ITS tracks V1");
139        delete iotrack;
140        return 2;
141     }
142    
143    if (gAlice->TreeE())//check if tree E exists
144      {
145       Nevents = (Int_t)gAlice->TreeE()->GetEntries();//if yes get number of events in gAlice
146       cout<<"________________________________________________________\n";
147       cout<<"Found "<<Nevents<<" event(s) in directory "<<GetDirName(currentdir)<<endl;
148       cout<<"Setting Magnetic Field. Factor is "<<gAlice->Field()->Factor()<<endl;
149       AliKalmanTrack::SetConvConst(100/0.299792458/0.2/gAlice->Field()->Factor());
150      }
151     else
152      {//if not return an error
153        Error("Read","Can not find Header tree (TreeE) in gAlice");
154        delete iotrack;
155        return 4;
156      }
157     
158    Int_t naccepted = 0;
159    char tname[30];
160    
161    for (Int_t currentEvent = 0; currentEvent < Nevents; currentEvent++)
162     { 
163       cout<<"Reading Event "<<currentEvent;
164       
165       sprintf(tname,"TreeT%d",currentEvent);
166       file->cd(); 
167       TTree *tracktree=(TTree*)file->Get(tname);
168       TBranch *tbranch=tracktree->GetBranch("ITStracks");
169       tbranch->SetAddress(&iotrack);
170       
171       gAliceFile->cd();
172       gAlice->GetEvent(currentEvent);
173       gAlice->Particles();
174
175       Int_t nentr=(Int_t)tracktree->GetEntries();
176       
177       cout<<".  Found "<<nentr<<" tracks.";
178       fflush(0);
179       
180       for (Int_t i=0; i<nentr; i++) 
181        {
182
183         tracktree->GetEvent(i);
184         if(!iotrack) continue;       
185         Int_t label = iotrack->GetLabel();
186         if (label < 0) 
187          {
188            delete iotrack;
189            continue;
190          }
191         TParticle *p = (TParticle*)gAlice->Particle(label);
192         if(!p)
193          {
194            Warning("Read","Can not get particle with label &d",label);
195            continue;
196          }
197         if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type
198                                            //if not take next partilce
199
200         AliHBTParticle* part = new AliHBTParticle(*p);
201         if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
202                                                 //if it does not delete it and take next good track
203         
204         Double_t px=iotrack->GetPx();
205         Double_t py=iotrack->GetPy();
206         Double_t pz=iotrack->GetPz();
207         Double_t mass = p->GetMass();
208         Double_t tEtot = TMath::Sqrt(px*px + py*py + pz*pz + mass*mass);//total energy of the track
209  
210         Double_t x= iotrack->GetX();
211         Double_t y= iotrack->GetY();
212         Double_t z= iotrack->GetZ();
213         
214         AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(), px, py , pz, tEtot, x, y, z, 0.);
215         if(Pass(track)) { delete  track;continue;}//check if meets all criteria of any of our cuts
216                                                   //if it does not delete it and take next good track
217
218         particles->AddParticle(totalNevents,part);//put track and particle on the run
219         tracks->AddParticle(totalNevents,track);
220         naccepted++;
221        }//end loop over tracks in the event
222
223        totalNevents++;
224        cout<<"  Accepted "<<naccepted<<" tracks"<<endl;
225      }//end of loop over events in current directory
226     
227     gAliceFile->Close();
228     delete gAliceFile;
229     gAliceFile = 0;
230     
231     file->Close(); 
232     delete file;
233     file = 0;
234     currentdir++;
235    }while(currentdir < Ndirs);//end of loop over directories specified in fDirs Obj Array
236
237
238   delete iotrack;
239   fIsRead = kTRUE;
240   return 0;
241  
242  }
243 /********************************************************************/
244
245 TFile* AliHBTReaderITSv1::OpenTrackFile(Int_t ndir)
246 {
247 //opens files to be read for given directoru nomber in fDirs Array
248    const TString& dirname = GetDirName(ndir); 
249    if (dirname == "")
250     {
251       Error("OpenGAliceFile","Can not get directory name");
252       return 0x0;
253     }
254    TString filename = dirname + "/" + fITSTracksFileName;
255
256    TFile *file = TFile::Open(filename.Data());   
257    if (!file)
258     {
259       Error("Read","Can not open file %s",filename.Data());
260       return 0x0;
261     }
262    if (!file->IsOpen())
263     {
264       Error("Read","Can not open file %s",filename.Data());
265       return 0x0;
266     }
267    
268    return file;
269 }
270
271
272 /********************************************************************/
273 TFile* AliHBTReaderITSv1::OpenGAliceFile(Int_t ndir)
274 {
275   const TString& dirname = GetDirName(ndir); 
276    if (dirname == "")
277     {
278       Error("OpenGAliceFile","Can not get directory name");
279       return 0x0;
280     }
281   
282   TString filename = dirname + "/" + fGAliceFileName;
283
284   TFile* gAliceFile = TFile::Open(filename.Data());
285   if ( gAliceFile== 0x0)
286    {
287      Error("OpenFiles","Can't open file named %s",filename.Data());
288      return 0x0;
289    }
290   if (!gAliceFile->IsOpen())
291    {
292      Error("OpenFiles","Can't open file named %s",filename.Data());
293      return 0x0;
294    }
295
296   if (!(gAlice=(AliRun*)gAliceFile->Get("gAlice")))
297    {
298      Error("OpenFiles","gAlice have not been found on %s !\n",filename.Data());
299      gAliceFile->Close();
300      delete gAliceFile;
301      return 0x0;
302    }
303   
304   return gAliceFile;
305 }
306
307 /********************************************************************/
308 /********************************************************************/
309
310