ReaderInternal: Internal data format implementation
[u/mrichter/AliRoot.git] / HBTAN / AliHBTReaderITSv2.cxx
1
2 #include "AliHBTReaderITSv2.h"
3
4 #include <iostream.h>
5 #include <fstream.h>
6 #include <TString.h>
7 #include <TObjString.h>
8 #include <TTree.h>
9 #include <TFile.h>
10 #include <TParticle.h>
11
12 #include <AliRun.h>
13 #include <AliMagF.h>
14 #include <AliITStrackV2.h>
15 //#include <AliITSParam.h>
16 #include <AliITStrackerV2.h>
17 #include <AliITSgeom.h>
18
19 #include "AliHBTRun.h"
20 #include "AliHBTEvent.h"
21 #include "AliHBTParticle.h"
22 #include "AliHBTParticleCut.h"
23
24
25 ClassImp(AliHBTReaderITSv2)
26
27 AliHBTReaderITSv2::
28  AliHBTReaderITSv2(const Char_t* trackfilename, const Char_t* clusterfilename,
29         const Char_t* galicefilename)
30                   :fTrackFileName(trackfilename),fClusterFileName(clusterfilename),
31                     fGAliceFileName(galicefilename)
32 {
33   //constructor, 
34   //Defaults:
35   //  trackfilename = "AliITStracksV2.root"
36   //  clusterfilename = "AliITSclustersV2.root"
37   //  galicefilename = "galice.root"
38   fParticles = new AliHBTRun();
39   fTracks    = new AliHBTRun();
40   fIsRead = kFALSE;
41 }
42 /********************************************************************/
43
44 AliHBTReaderITSv2::
45  AliHBTReaderITSv2(TObjArray* dirs, const Char_t* trackfilename, 
46                    const Char_t* clusterfilename, const Char_t* galicefilename)
47                   : AliHBTReader(dirs),
48                     fTrackFileName(trackfilename),fClusterFileName(clusterfilename),
49                     fGAliceFileName(galicefilename)
50 {
51   //constructor, 
52   //Defaults:
53   //  trackfilename = "AliITStracksV2.root"
54   //  clusterfilename = "AliITSclustersV2.root"
55   //  galicefilename = "galice.root"
56   
57   fParticles = new AliHBTRun();
58   fTracks    = new AliHBTRun();
59   fIsRead = kFALSE;
60 }
61 /********************************************************************/
62
63 AliHBTReaderITSv2::~AliHBTReaderITSv2()
64  {
65    if (fParticles) delete fParticles;
66    if (fTracks) delete fTracks;
67  }
68 /********************************************************************/
69 /********************************************************************/
70
71 AliHBTEvent* AliHBTReaderITSv2::GetParticleEvent(Int_t n)
72  {
73  //returns Nth event with simulated particles
74    if (!fIsRead) 
75     if(Read(fParticles,fTracks))
76      {
77        Error("GetParticleEvent","Error in reading");
78        return 0x0;
79      }
80
81    return fParticles->GetEvent(n);
82  }
83 /********************************************************************/
84
85 AliHBTEvent* AliHBTReaderITSv2::GetTrackEvent(Int_t n)
86  {
87  //returns Nth event with reconstructed tracks
88    if (!fIsRead) 
89     if(Read(fParticles,fTracks))
90      {
91        Error("GetTrackEvent","Error in reading");
92        return 0x0;
93      }
94    return fTracks->GetEvent(n);
95  }
96 /********************************************************************/
97
98 Int_t AliHBTReaderITSv2::GetNumberOfPartEvents()
99  {
100  //returns number of events of particles
101    if (!fIsRead)
102     if(Read(fParticles,fTracks))
103      {
104        Error("GetNumberOfPartEvents","Error in reading");
105        return 0;
106      }
107    return fParticles->GetNumberOfEvents();
108  }
109
110 /********************************************************************/
111 Int_t AliHBTReaderITSv2::GetNumberOfTrackEvents()
112  {
113  //returns number of events of tracks
114   if (!fIsRead) 
115     if(Read(fParticles,fTracks))
116      {
117        Error("GetNumberOfTrackEvents","Error in reading");
118        return 0;
119      }
120   return fTracks->GetNumberOfEvents();
121  }
122
123
124 /********************************************************************/
125 /********************************************************************/
126 Int_t AliHBTReaderITSv2::Read(AliHBTRun* particles, AliHBTRun *tracks)
127 {
128  Int_t Nevents = 0; //number of events found in given directory
129  Int_t Ndirs; //number of the directories to be read
130  Int_t Ntracks; //number of tracks in current event
131  Int_t currentdir = 0; //number of events in the current directory 
132  Int_t totalNevents = 0; //total number of events read from all directories up to now
133  register Int_t i; //iterator
134  
135  TFile *aTracksFile;//file with tracks
136  TFile *aClustersFile;//file with clusters
137  TFile *aGAliceFile;//file name with galice
138
139 // AliITStrackerV2 *tracker; // ITS tracker - used for cooking labels
140  TTree *tracktree; // tree for tracks
141  
142  Double_t xk;
143  Double_t par[5]; //Kalman track parameters
144  Float_t phi, lam, pt;//angles and transverse momentum
145  Int_t label; //label of the current track
146
147  char tname[100]; //buffer for tree name
148  AliITStrackV2 *iotrack= new AliITStrackV2(); //buffer track for reading data from tree
149
150  if (!particles) //check if an object is instatiated
151   {
152     Error("Read"," particles object must instatiated before passing it to the reader");
153   }
154  if (!tracks)  //check if an object is instatiated
155   {
156     Error("Read"," tracks object must instatiated before passing it to the reader");
157   }
158  particles->Reset();//clear runs == delete all old events
159  tracks->Reset();
160
161  if (fDirs) //if array with directories is supplied by user
162   {
163     Ndirs = fDirs->GetEntries(); //get the number if directories
164   }
165  else
166   {
167     Ndirs = 0; //if the array is not supplied read only from current directory
168   }
169  
170 // cout<<"Found "<<Ndirs<<" directory entries"<<endl;
171  
172  do //do while is good even if Ndirs==0 (than read from current directory)
173    {
174     if( (i=OpenFiles(aTracksFile,aClustersFile,aGAliceFile,currentdir)) )
175      {
176        Error("Read","Exiting due to problems with opening files. Errorcode %d",i);
177        delete iotrack;
178        return i;
179      }
180     
181     if (gAlice->TreeE())//check if tree E exists
182      {
183       Nevents = (Int_t)gAlice->TreeE()->GetEntries();//if yes get number of events in gAlice
184       cout<<"________________________________________________________\n";
185       cout<<"Found "<<Nevents<<" event(s) in directory "<<GetDirName(currentdir)<<endl;
186       cout<<"Setting Magnetic Field. Factor is "<<gAlice->Field()->Factor()<<endl;
187       AliKalmanTrack::SetConvConst(100/0.299792458/0.2/gAlice->Field()->Factor());  
188      }
189     else
190      {//if not return an error
191        Error("Read","Can not find Header tree (TreeE) in gAlice");
192        delete iotrack;
193        return 1;
194      }
195     
196     AliITSgeom *geom=(AliITSgeom*)aClustersFile->Get("AliITSgeom");
197     if (!geom) 
198      { 
199        Error("Read","Can't get the ITS geometry!"); 
200        delete iotrack;
201        return 2; 
202      }
203
204     for(Int_t currentEvent =0; currentEvent<Nevents;currentEvent++)//loop over all events
205      {
206        cout<<"Reading Event "<<currentEvent<<endl;
207        
208        aGAliceFile->cd();
209        gAlice->GetEvent(currentEvent);
210
211 //       TParticle * part = gAlice->Particle(0);
212       
213        aClustersFile->cd();
214 //       Double_t orz=part->Vz();
215 //       tracker = new AliITStrackerV2(geom,currentEvent,orz); //<---- this is for Massimo version
216 //       tracker = new AliITStrackerV2(geom,currentEvent);
217        sprintf(tname,"TreeT_ITS_%d",currentEvent);
218        
219        tracktree=(TTree*)aTracksFile->Get(tname);
220        if (!tracktree) 
221          {
222            Error("Read","Can't get a tree with ITS tracks"); 
223            delete iotrack;
224  //          delete tracker;
225            return 4;
226          }
227        TBranch *tbranch=tracktree->GetBranch("tracks");
228       
229        Ntracks=(Int_t)tracktree->GetEntries();
230
231        Int_t accepted = 0;
232        Int_t tpcfault = 0;
233        Int_t itsfault = 0;
234        for (i=0; i<Ntracks; i++) //loop over all tpc tracks
235         { 
236           if(i%100 == 0)cout<<"all: "<<i<<"   accepted: "<<accepted<<"   tpc faults: "<<tpcfault<<"\r";
237           
238           tbranch->SetAddress(&iotrack);
239           tracktree->GetEvent(i);
240
241           label=iotrack->GetLabel();
242           if (label < 0) 
243            {
244              tpcfault++;
245              continue;
246            }
247
248           TParticle *p = (TParticle*)gAlice->Particle(label);
249           if(p == 0x0) continue; //if returned pointer is NULL
250           if(p->GetPDG() == 0x0) continue; //if particle has crezy PDG code (not known to our database)
251
252           if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type 
253                                               //if not take next partilce
254             
255           AliHBTParticle* part = new AliHBTParticle(*p);
256           if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
257                                                   //if it does not delete it and take next good track
258
259           
260           iotrack->PropagateTo(3.,0.0028,65.19);
261           iotrack->PropagateToVertex();
262  
263           iotrack->GetExternalParameters(xk,par);     //get properties of the track
264           phi=TMath::ASin(par[2]) + iotrack->GetAlpha(); 
265           if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
266           if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
267           lam=par[3]; 
268           pt=1.0/TMath::Abs(par[4]);
269             
270           Double_t tpx = pt * TMath::Cos(phi); //track x coordinate of momentum
271           Double_t tpy = pt * TMath::Sin(phi); //track y coordinate of momentum
272           Double_t tpz = pt * lam; //track z coordinate of momentum
273            
274           Double_t mass = p->GetMass();
275           Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);//total energy of the track
276             
277           AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(), tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
278           if(Pass(track))//check if meets all criteria of any of our cuts
279                          //if it does not delete it and take next good track
280            { 
281             delete track;
282             delete part;
283             continue;
284            }
285           particles->AddParticle(totalNevents,part);//put track and particle on the run
286           tracks->AddParticle(totalNevents,track);
287           accepted++;
288         }//end of loop over tracks in the event
289         
290        aTracksFile->Delete(tname);
291        aTracksFile->Delete("tracks");
292 //       delete tracker;
293        
294        totalNevents++;
295        cout<<"all: "<<i<<"   accepted: "<<accepted<<"   tpc faults: "<<tpcfault<<"   its faults: "<<itsfault<<endl;
296      
297      }//end of loop over events in current directory
298     CloseFiles(aTracksFile,aClustersFile,aGAliceFile);     
299     currentdir++;
300    }while(currentdir < Ndirs);//end of loop over directories specified in fDirs Obj Array
301
302  delete iotrack;
303  fIsRead = kTRUE;
304  return 0;
305 }
306
307 /********************************************************************/
308 Int_t AliHBTReaderITSv2::OpenFiles
309 (TFile*& aTracksFile, TFile*& aClustersFile, TFile*& agAliceFile,Int_t event)
310 {
311  //opens all the files
312    
313    
314    const TString& dirname = GetDirName(event); 
315    if (dirname == "")
316     {
317       Error("OpenFiles","Can not get directory name");
318       return 4;
319     }
320    
321    TString filename = dirname +"/"+ fTrackFileName;
322    aTracksFile = TFile::Open(filename.Data());
323    if ( aTracksFile  == 0x0 ) 
324      {
325        Error("OpenFiles","Can't open file with tacks named %s",filename.Data());
326        return 1;
327      }
328    if (!aTracksFile->IsOpen())
329      {
330        Error("OpenFiles","Can't open file with tacks named %s",filename.Data());
331        return 1;
332      }
333   
334    filename = dirname +"/"+ fClusterFileName;
335    aClustersFile = TFile::Open(filename.Data());
336    if ( aClustersFile == 0x0 )
337     {
338       Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
339       return 2;
340     }
341    if (!aClustersFile->IsOpen())
342     {
343       Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
344       return 2;
345     }
346
347    filename = dirname +"/"+ fGAliceFileName;
348    agAliceFile = TFile::Open(filename.Data());
349    if ( agAliceFile== 0x0)
350     {
351       Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
352       return 3;
353     }
354    if (!agAliceFile->IsOpen())
355     {
356       Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
357       return 3;
358     } 
359    
360    if (!(gAlice=(AliRun*)agAliceFile->Get("gAlice"))) 
361     {
362       Error("OpenFiles","gAlice have not been found on %s !\n",filename.Data());
363       return 5;
364     }
365
366    return 0; 
367 }
368 /********************************************************************/
369
370 /********************************************************************/
371   
372 void AliHBTReaderITSv2::CloseFiles(TFile*& tracksFile, TFile*& clustersFile, TFile*& gAliceFile)
373 {
374   //closes the files
375   tracksFile->Close();
376   delete tracksFile;
377   tracksFile = 0x0;
378   
379   clustersFile->Close();
380   delete clustersFile;
381   clustersFile = 0x0;
382   
383   delete gAlice;
384   gAlice = 0;
385
386   if (gAliceFile) 
387    {
388      gAliceFile->Close();
389      delete gAliceFile;
390      gAliceFile = 0x0;
391    }
392 }
393
394 /********************************************************************/
395
396