SetDisplayInfo added
[u/mrichter/AliRoot.git] / HBTAN / AliHBTReaderKineTree.cxx
1 #include "AliHBTReaderKineTree.h"
2
3 #include <Riostream.h>
4 //#include <Riostream.h>
5 #include <TString.h>
6 #include <TObjString.h>
7 #include <TTree.h>
8 #include <TFile.h>
9 #include <TParticle.h>
10
11 #include <AliRun.h>
12 #include <AliStack.h>
13 #include <AliHeader.h>
14
15 #include "AliHBTRun.h"
16 #include "AliHBTEvent.h"
17 #include "AliHBTParticle.h"
18 #include "AliHBTParticleCut.h"
19
20 ClassImp(AliHBTReaderKineTree)
21 /**********************************************************/
22
23 AliHBTReaderKineTree::
24 AliHBTReaderKineTree():fFileName("galice.root")
25 {
26   fParticles = new AliHBTRun();
27   fIsRead = kFALSE;
28 }
29
30 AliHBTReaderKineTree::
31 AliHBTReaderKineTree(TString& fname):fFileName(fname)
32 {
33   fParticles = new AliHBTRun();
34   fIsRead = kFALSE;
35 }
36
37
38 /**********************************************************/
39 AliHBTReaderKineTree::
40 AliHBTReaderKineTree(TObjArray* dirs,const Char_t *filename):AliHBTReader(dirs),fFileName(filename)
41 {
42   fParticles = new AliHBTRun();
43   fIsRead = kFALSE;
44 }
45 /**********************************************************/
46
47 AliHBTEvent* AliHBTReaderKineTree::GetParticleEvent(Int_t n)
48  {
49  //returns Nth event with simulated particles
50    if (!fIsRead) 
51     if(Read(fParticles,0x0))
52      {
53        Error("GetParticleEvent","Error in reading");
54        return 0x0;
55      }
56
57    return fParticles->GetEvent(n);
58  }
59
60 /********************************************************************/
61
62 Int_t AliHBTReaderKineTree::GetNumberOfPartEvents()
63  {
64  //returns number of events of particles
65    if (!fIsRead)
66     if(Read(fParticles,0x0))
67      {
68        Error("GetNumberOfPartEvents","Error in reading");
69        return 0;
70      }
71    return fParticles->GetNumberOfEvents();
72  }
73
74
75 /**********************************************************/
76 Int_t AliHBTReaderKineTree::
77 Read(AliHBTRun* particles, AliHBTRun *tracks)
78  {
79  cout<<"AliHBTReaderKineTree::Read()"<<endl;
80  if (!particles) //check if an object is instatiated
81    {
82      Error("Read"," particles object must instatiated before passing it to the reader");
83      return 1;
84    }  
85  particles->Reset();//clear runs == delete all old events
86  
87  Int_t Ndirs;//total number of directories specified in fDirs array
88  Int_t Nevents; //number of events read in current directory
89  Int_t totalNevents = 0; //total number of read events 
90  Int_t currentdir = 0; //number of current directory name is fDirs array
91
92  if (fDirs) //if array with directories is supplied by user
93   {
94     Ndirs = fDirs->GetEntries(); //get the number if directories
95   }
96  else
97   {
98     Ndirs = 0; //if the array is not supplied read only from current directory
99   }
100
101  do  //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
102   { 
103     cout<<"________________________________________________________\n";
104     TFile* fileK = OpenFile(currentdir);
105     TTree* treeE = (TTree*)fileK->Get("TE");
106     if(treeE) 
107      {
108       Nevents = (Int_t)treeE->GetEntries();
109       cout<<"Found "<<Nevents<<" in directory "<<GetDirName(currentdir)<<endl;
110       fileK->Delete("TE");
111      }
112     else
113      {
114       Error("Read","Cannot find Header Tree (TE)");
115       Nevents = 0;
116      }
117     
118     for(Int_t currentEvent =0; currentEvent<Nevents;currentEvent++)//loop over all events
119      {
120 //      cout<<"processing event "<<currentEvent<<" in current dir."<<endl;
121       AliStack* stack = GetStack(currentEvent,fileK);
122       if (!stack)
123        {
124          Error("Read","Can not get stack for event %d",currentEvent);
125          continue;
126        }
127       Int_t npart = stack->GetNtrack();
128       Int_t nnn = 0;
129       for (Int_t i = 0;i<npart; i++)
130        {
131          
132          TParticle * p = stack->Particle(i);
133          if (p->GetFirstMother() >= 0) continue;
134          
135          if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type 
136                                               //if not take next partilce
137          
138          AliHBTParticle* part = new AliHBTParticle(*p);
139          if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
140                                                   //if it does not delete it and take next good track
141          particles->AddParticle(totalNevents,part);//put track and particle on the run
142
143          if ( (nnn%100) == 0) 
144           {
145             cout<<nnn<<"\r";
146             fflush(0);
147           }
148          nnn++;
149        }
150       cout<<"Total read "<<nnn<<endl;
151       delete stack;
152       totalNevents++;
153      }
154     delete gAlice;
155     gAlice = 0;
156     if (fileK) 
157      {
158       fileK->Close();
159       delete fileK;
160       fileK = 0;
161      }
162     
163     currentdir++;
164   }while(currentdir < Ndirs);//end of loop over directories specified in fDirs Obj Array
165  fIsRead = kTRUE;
166  return 0;
167 }
168
169 /**********************************************************/
170 AliStack* AliHBTReaderKineTree::GetStack(Int_t n, TFile* file)
171 {
172   AliHeader *header =  new AliHeader();
173   TTree* treeE = (TTree*)file->Get("TE");
174   if(treeE) 
175    {
176     treeE->SetBranchAddress("Header", &header);
177     if (!treeE->GetEntry(n)) 
178      {
179        Error("GetEvent","Cannot find event:%dn",n);
180        delete header;
181        treeE->Delete("TE");
182        return 0x0;
183      }
184    }  
185   else 
186    {
187      Error("GetStack","Cannot find Header Tree (TE)n");
188      delete header;
189      return 0x0;
190    }
191
192   AliStack* stack = header->Stack();
193   if (stack) 
194    {
195     if (!stack->GetEvent(n)) 
196      {
197        delete header;
198        treeE->Delete("TE");
199        return 0x0;
200      }
201    }
202   delete header;
203   treeE->Delete("TE");
204   return stack;
205 }
206
207 TFile* AliHBTReaderKineTree::OpenFile(Int_t n)
208 {
209 //opens file with kine tree
210
211  const TString& dirname = GetDirName(n);
212  if (dirname == "")
213   {
214    Error("OpenFiles","Can not get directory name");
215    return 0x0;
216   }
217  TString filename = dirname +"/"+ fFileName;
218  TFile *ret = TFile::Open(filename.Data()); 
219
220  if ( ret == 0x0)
221   {
222     Error("OpenFiles","Can't open file %s",filename.Data());
223     return 0x0;
224   }
225  if (!ret->IsOpen())
226   {
227     Error("OpenFiles","Can't open file  %s",filename.Data());
228     return 0x0;
229   }
230  
231  return ret;
232 }