Monitor functions implemented by Z.Chajecki.
[u/mrichter/AliRoot.git] / HBTAN / AliHBTReaderInternal.cxx
1 #include "AliHBTReaderInternal.h"
2
3 #include <iostream.h>
4 //#include <fstream.h>
5 #include <TTree.h>
6 #include <TFile.h>
7 #include <TParticle.h>
8
9 #include <AliRun.h>
10 #include <AliMagF.h>
11
12 #include "AliHBTRun.h"
13 #include "AliHBTEvent.h"
14 #include "AliHBTParticle.h"
15 #include "AliHBTParticleCut.h"
16
17
18 ClassImp(AliHBTReaderInternal)
19 /********************************************************************/
20
21 AliHBTReaderInternal::AliHBTReaderInternal()
22 {
23   fParticles = 0x0; 
24   fTracks = 0x0;
25   fIsRead = kFALSE;
26 }
27 /********************************************************************/
28
29 AliHBTReaderInternal::AliHBTReaderInternal(const char *filename):fFileName(filename)
30
31  fParticles = new AliHBTRun();
32  fTracks    = new AliHBTRun();
33  fIsRead = kFALSE;
34 }
35 /********************************************************************/
36 AliHBTReaderInternal::AliHBTReaderInternal(TObjArray* dirs, const char *filename):
37   AliHBTReader(dirs),fFileName(filename)
38
39  fParticles = new AliHBTRun();
40  fTracks    = new AliHBTRun();
41  fIsRead = kFALSE;
42 }
43 AliHBTReaderInternal::~AliHBTReaderInternal()
44  {
45  //desctructor
46    delete fParticles;
47    delete fTracks;
48  }
49 /********************************************************************/
50 AliHBTEvent* AliHBTReaderInternal::GetParticleEvent(Int_t n)
51  {
52  //returns Nth event with simulated particles
53    if (!fIsRead) 
54     if(Read(fParticles,fTracks))
55      {
56        Error("GetParticleEvent","Error in reading");
57        return 0x0;
58      }
59    return fParticles->GetEvent(n);
60  }
61 /********************************************************************/
62 AliHBTEvent* AliHBTReaderInternal::GetTrackEvent(Int_t n)
63  {
64  //returns Nth event with reconstructed tracks
65    if (!fIsRead) 
66     if(Read(fParticles,fTracks))
67      {
68        Error("GetTrackEvent","Error in reading");
69        return 0x0;
70      }
71    return fTracks->GetEvent(n);
72  }
73 /********************************************************************/
74
75 Int_t AliHBTReaderInternal::GetNumberOfPartEvents()
76  {
77  //returns number of events of particles
78    if (!fIsRead) 
79     if ( Read(fParticles,fTracks))
80      {
81        Error("GetNumberOfPartEvents","Error in reading");
82        return 0;
83      }
84    return fParticles->GetNumberOfEvents();
85  }
86
87 /********************************************************************/
88 Int_t AliHBTReaderInternal::GetNumberOfTrackEvents()
89  {
90  //returns number of events of tracks
91   if (!fIsRead)
92     if(Read(fParticles,fTracks))
93      {
94        Error("GetNumberOfTrackEvents","Error in reading");
95        return 0;
96      }
97   return fTracks->GetNumberOfEvents();
98  }
99 /********************************************************************/
100
101
102 Int_t AliHBTReaderInternal::Read(AliHBTRun* particles, AliHBTRun *tracks)
103  {
104  //reads data and puts put to the particles and tracks objects
105  //reurns 0 if everything is OK
106  //
107   cout<<"AliHBTReaderInternal::Read()"<<endl;
108   Int_t i; //iterator and some temprary values
109   Int_t totalNevents = 0; //total number of read events
110   Int_t Nevents = 0;
111   Int_t currentdir = 0;
112   Int_t Ndirs;
113   
114   TFile *aFile;//file with tracks
115   AliHBTParticle* tpart = 0x0, *ttrack = 0x0;
116   
117   if (!particles) //check if an object is instatiated
118    {
119      Error("Read"," particles object must instatiated before passing it to the reader");
120    }
121   if (!tracks)  //check if an object is instatiated
122    {
123      Error("Read"," tracks object must instatiated before passing it to the reader");
124    }
125   particles->Reset();//clear runs == delete all old events
126   tracks->Reset();
127  
128   if (fDirs) //if array with directories is supplied by user
129    {
130      Ndirs = fDirs->GetEntries(); //get the number if directories
131    }
132   else
133    {
134      Ndirs = 0; //if the array is not supplied read only from current directory
135    }
136
137   TClonesArray* pbuffer = new TClonesArray("AliHBTParticle",15000);
138   TClonesArray* tbuffer = new TClonesArray("AliHBTParticle",15000);
139
140   do  //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
141    {
142     
143     if( (i=OpenFile(aFile, currentdir)) )
144      {
145        Error("Read","Exiting due to problems with opening files. Errorcode %d",i);
146        return i;
147      }
148    /***************************/
149    /***************************/
150    /***************************/
151     
152      TTree* tree = (TTree*)aFile->Get("data");
153      if (tree == 0x0)
154       {
155        Error("Read","Can not get the tree");
156        return 1;
157       }
158      
159      TBranch *trackbranch=tree->GetBranch("tracks");//get the branch with tracks
160      if (trackbranch == 0x0) ////check if we got the branch
161        {//if not return with error
162          Warning("Read","Can't find a branch with tracks !\n"); 
163        }
164      else
165       {
166         trackbranch->SetAddress(&tbuffer);
167       }
168
169      TBranch *partbranch=tree->GetBranch("particles");//get the branch with particles
170      if (partbranch == 0x0) ////check if we got the branch
171        {//if not return with error
172          Warning("Read","Can't find a branch with particles !\n"); 
173        }
174      else
175       {
176         partbranch->SetAddress(&pbuffer);
177       }
178      
179      Nevents = (Int_t)tree->GetEntries();
180      cout<<"________________________________________________________\n";
181      cout<<"Found "<<Nevents<<" event(s) in directory "<<GetDirName(currentdir)<<endl;
182      
183      for (Int_t currentEvent =0; currentEvent<Nevents;currentEvent++)
184       {
185        cout<<"Event "<<currentEvent<<"\n";
186        tree->GetEvent(currentEvent);
187          
188        if (partbranch && trackbranch)
189         {
190            for(i = 0; i < pbuffer->GetEntries(); i++)
191              {
192                cout<<i<<"\r";
193                tpart = dynamic_cast<AliHBTParticle*>(pbuffer->At(i));
194                ttrack =  dynamic_cast<AliHBTParticle*>(tbuffer->At(i));
195
196                if( (tpart == 0x0) || (ttrack == 0x0) ) continue; //if returned pointer is NULL
197                if( (tpart->GetPDG()==0x0)||(ttrack->GetPDG()==0x0) ) continue; //if particle has crezy PDG code (not known to our database)
198                if((Pass(tpart))||(Pass(ttrack))) continue; //check if we are intersted with particles of this type 
199                                                           //if not take next partilce
200                AliHBTParticle* part = new AliHBTParticle(*tpart);
201                AliHBTParticle* track = new AliHBTParticle(*ttrack);
202                particles->AddParticle(totalNevents,part);//put track and particle on the run
203                tracks->AddParticle(totalNevents,track);
204              }
205             cout<<"Read: "<<particles->GetNumberOfParticlesInEvent(totalNevents)<<" particles and tracks \n";
206         }
207        else
208         { 
209          if (partbranch)
210           {
211             for(i = 0; i < pbuffer->GetEntries(); i++)
212              {
213                cout<<i<<"\r";
214                tpart = dynamic_cast<AliHBTParticle*>(pbuffer->At(i));
215                if(tpart == 0x0) continue; //if returned pointer is NULL
216                if(tpart->GetPDG() == 0x0) continue; //if particle has crezy PDG code (not known to our database)
217                if(Pass(tpart)) continue; //check if we are intersted with particles of this type 
218                                          //if not take next partilce
219  
220                AliHBTParticle* part = new AliHBTParticle(*tpart);
221                particles->AddParticle(totalNevents,part);//put track and particle on the run
222              }
223             cout<<"\nRead: "<<particles->GetNumberOfParticlesInEvent(totalNevents)<<" particles  \n";
224           }
225          else cout<<"Read: 0 particles  \n";
226          
227          if (trackbranch)
228           {
229             for(i = 0; i < tbuffer->GetEntries(); i++)
230              {
231                cout<<i<<"\r";
232                tpart = dynamic_cast<AliHBTParticle*>(tbuffer->At(i));
233                if(tpart == 0x0) continue; //if returned pointer is NULL
234                if(tpart->GetPDG() == 0x0) continue; //if particle has crezy PDG code (not known to our database)
235                if(Pass(tpart)) continue; //check if we are intersted with particles of this type 
236                                                    //if not take next partilce
237                AliHBTParticle* part = new AliHBTParticle(*tpart);
238                tracks->AddParticle(totalNevents,part);//put track and particle on the run
239              }
240             cout<<"\nRead: "<<tracks->GetNumberOfParticlesInEvent(totalNevents)<<" tracks"<<endl;
241           }
242          else cout<<" 0 tracks"<<endl;
243         }
244         totalNevents++;
245          
246        }
247
248    /***************************/
249    /***************************/
250    /***************************/
251    currentdir++;
252    aFile->Close();
253    aFile = 0x0;
254
255    }while(currentdir < Ndirs);
256
257   fIsRead = kTRUE;
258   return 0;
259 }
260
261 /********************************************************************/
262
263 Int_t AliHBTReaderInternal::OpenFile(TFile*& aFile,Int_t event)
264 {
265
266    const TString& dirname = GetDirName(event); 
267    if (dirname == "")
268     {
269       Error("OpenFile","Can not get directory name");
270       return 4;
271     }
272    
273    TString filename = dirname +"/"+ fFileName;
274    aFile = TFile::Open(filename.Data());
275    if ( aFile  == 0x0 ) 
276      {
277        Error("OpenFiles","Can't open file with tacks named %s",filename.Data());
278        return 1;
279      }
280    if (!aFile->IsOpen())
281      {
282        Error("OpenFiles","Can't open file with tacks named %s",filename.Data());
283        return 1;
284      }
285    return 0; 
286 }
287
288
289
290
291 Int_t AliHBTReaderInternal::Write(AliHBTReader* reader,const char* outfile)
292  {
293   //reads tracks from runs and writes them to file
294     Int_t i,j;
295   
296   
297   cout<<"________________________________________________________\n";
298   cout<<"________________________________________________________\n";
299   cout<<"________________________________________________________\n";
300
301   TFile *histoOutput = TFile::Open(outfile,"recreate");
302   
303   if (!histoOutput->IsOpen())
304    {
305      cout<<"File is not opened"<<endl;
306      return 1;
307    }
308     
309
310   TTree *tracktree = new TTree("data","Tree with tracks");
311
312   TClonesArray* pbuffer = new TClonesArray("AliHBTParticle",15000);
313   TClonesArray* tbuffer = new TClonesArray("AliHBTParticle",15000);
314
315   TClonesArray &particles = *pbuffer;
316   TClonesArray &tracks = *tbuffer;
317     
318   TString name("Tracks");
319  
320   
321   Int_t NT = reader->GetNumberOfTrackEvents();
322   Int_t NP = reader->GetNumberOfPartEvents();
323   
324   Bool_t trck = (NT > 0) ? kTRUE : kFALSE;
325   Bool_t part = (NP > 0) ? kTRUE : kFALSE;
326
327   TBranch *trackbranch = 0x0, *partbranch = 0x0;
328   
329   
330   if (trck) trackbranch = tracktree->Branch("tracks","TClonesArray",&tbuffer);
331   if (part) partbranch = tracktree->Branch("particles","TClonesArray",&pbuffer);
332
333
334   
335   if ( (trck) && (part) && (NP != NT))
336    {
337      cerr<<"Warning number of track and particle events is different"<<endl;
338    }
339   
340   Int_t N;
341   if (NT >= NP ) N = NT; else N = NP;
342
343   for ( i =0;i< N; i++)
344     {
345       cout<<"Event "<<i+1<<"\n";
346       if (trck && (i<=NT))
347        {
348          cout<<"Tracks: \n";
349          AliHBTEvent* trackev = reader->GetTrackEvent(i);
350          for ( j = 0; j< trackev->GetNumberOfParticles();j++)
351           {
352             cout<<j<<"\r";
353             const AliHBTParticle& t= *(trackev->GetParticle(j));
354             new (tracks[j]) AliHBTParticle(t);
355           }
356
357        }else cout<<"NO TRACKS";
358       
359       cout<<endl;
360       
361       if (part && (i<=NP))
362        {
363         cout<<"Particles: \n";
364         AliHBTEvent* partev = reader->GetParticleEvent(i);
365         for ( j = 0; j< partev->GetNumberOfParticles();j++)
366          {
367            const AliHBTParticle& part= *(partev->GetParticle(j));
368            cout<<j<<"\r";
369            new (particles[j]) AliHBTParticle(part);
370          }
371        }else cout<<"NO PARTICLES";
372       cout<<endl;
373
374       histoOutput->cd();
375       tracktree->Fill();
376       tracktree->AutoSave();
377       tbuffer->Delete();
378       pbuffer->Delete();
379      }
380
381   histoOutput->cd();
382   tracktree->Write(0,TObject::kOverwrite);
383   delete tracktree;
384
385   tbuffer->SetOwner();
386   pbuffer->SetOwner();
387   delete pbuffer;
388   delete tbuffer;
389
390   histoOutput->Close();
391   return 0;
392  }