Lookup of both parameters setup
[u/mrichter/AliRoot.git] / HBTAN / AliHBTReaderTPC.cxx
1 #include "AliHBTReaderTPC.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 #include <AliTPCtrack.h>
12 #include <AliTPCParam.h>
13 #include <AliTPCtracker.h>
14
15 #include "AliHBTRun.h"
16 #include "AliHBTEvent.h"
17 #include "AliHBTParticle.h"
18 #include "AliHBTParticleCut.h"
19
20
21 ClassImp(AliHBTReaderTPC)
22 //reader for TPC tracking
23 //needs galice.root, AliTPCtracks.root, AliTPCclusters.root, good_tracks_tpc 
24 //
25 //more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html
26 //Piotr.Skowronski@cern.ch
27
28 AliHBTReaderTPC::
29  AliHBTReaderTPC(const Char_t* trackfilename,const Char_t* clusterfilename,
30                  const Char_t* galicefilename):
31                  fTrackFileName(trackfilename),fClusterFileName(clusterfilename),
32                  fGAliceFileName(galicefilename)
33 {
34   //constructor, 
35   //Defaults:
36   //  trackfilename = "AliTPCtracks.root"
37   //  clusterfilename = "AliTPCclusters.root"
38   //  galicefilename = ""  - this means: Do not open gAlice file - 
39   //                         just leave the global pointer untached
40   
41   fParticles = new AliHBTRun();
42   fTracks    = new AliHBTRun();
43   fIsRead = kFALSE;
44 }
45 /********************************************************************/
46 AliHBTReaderTPC::
47 AliHBTReaderTPC(TObjArray* dirs,
48                   const Char_t* trackfilename, const Char_t* clusterfilename,
49                   const Char_t* galicefilename):
50                   AliHBTReader(dirs), fTrackFileName(trackfilename),
51                   fClusterFileName(clusterfilename),fGAliceFileName(galicefilename)
52
53 {
54   //constructor, 
55   //Defaults:
56   //  trackfilename = "AliTPCtracks.root"
57   //  clusterfilename = "AliTPCclusters.root"
58   //  galicefilename = ""  - this means: Do not open gAlice file - 
59   //                         just leave the global pointer untached
60   
61   fParticles = new AliHBTRun();
62   fTracks    = new AliHBTRun();
63   fIsRead = kFALSE;
64   
65 }
66 /********************************************************************/
67
68 AliHBTReaderTPC::~AliHBTReaderTPC()
69  {
70  //desctructor
71    delete fParticles;
72    delete fTracks;
73  }
74 /********************************************************************/
75
76 AliHBTEvent* AliHBTReaderTPC::GetParticleEvent(Int_t n)
77  {
78  //returns Nth event with simulated particles
79    if (!fIsRead) 
80     if(Read(fParticles,fTracks))
81      {
82        Error("GetParticleEvent","Error in reading");
83        return 0x0;
84      }
85    return fParticles->GetEvent(n);
86  }
87 /********************************************************************/
88 AliHBTEvent* AliHBTReaderTPC::GetTrackEvent(Int_t n)
89  {
90  //returns Nth event with reconstructed tracks
91    if (!fIsRead) 
92     if(Read(fParticles,fTracks))
93      {
94        Error("GetTrackEvent","Error in reading");
95        return 0x0;
96      }
97    return fTracks->GetEvent(n);
98  }
99 /********************************************************************/
100
101 Int_t AliHBTReaderTPC::GetNumberOfPartEvents()
102  {
103  //returns number of events of particles
104    if (!fIsRead) 
105     if ( Read(fParticles,fTracks))
106      {
107        Error("GetNumberOfPartEvents","Error in reading");
108        return 0;
109      }
110    return fParticles->GetNumberOfEvents();
111  }
112
113 /********************************************************************/
114 Int_t AliHBTReaderTPC::GetNumberOfTrackEvents()
115  {
116  //returns number of events of tracks
117   if (!fIsRead)
118     if(Read(fParticles,fTracks))
119      {
120        Error("GetNumberOfTrackEvents","Error in reading");
121        return 0;
122      }
123   return fTracks->GetNumberOfEvents();
124  }
125 /********************************************************************/
126
127
128 Int_t AliHBTReaderTPC::Read(AliHBTRun* particles, AliHBTRun *tracks)
129  {
130  //reads data and puts put to the particles and tracks objects
131  //reurns 0 if everything is OK
132  //
133   cout<<"AliHBTReaderTPC::Read()"<<endl;
134   Int_t i; //iterator and some temprary values
135   Int_t Nevents = 0;
136   Int_t totalNevents = 0;
137   TFile *aTracksFile;//file with tracks
138   TFile *aClustersFile;//file with clusters
139   TFile *aGAliceFile;//!ile name with galice
140
141   if (!particles) //check if an object is instatiated
142    {
143      Error("Read"," particles object must instatiated before passing it to the reader");
144    }
145   if (!tracks)  //check if an object is instatiated
146    {
147      Error("Read"," tracks object must instatiated before passing it to the reader");
148    }
149   particles->Reset();//clear runs == delete all old events
150   tracks->Reset();
151  
152   TObjArray *tarray = new TObjArray(5000); //cotainer for tpc tracks
153   tarray->SetOwner(); //set the ownership of the objects it contains
154                       //when array is is deleted or cleared all objects 
155                       //that it contains are deleted
156   Int_t currentdir = 0;
157
158   Int_t Ndirs;
159   if (fDirs) //if array with directories is supplied by user
160    {
161      Ndirs = fDirs->GetEntries(); //get the number if directories
162    }
163   else
164    {
165      Ndirs = 0; //if the array is not supplied read only from current directory
166    }
167
168   do  //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
169    {
170     
171     if( (i=OpenFiles(aTracksFile,aClustersFile,aGAliceFile,currentdir)) )
172      {
173        Error("Read","Exiting due to problems with opening files. Errorcode %d",i);
174        currentdir++;
175        continue;
176      }
177   
178     
179     if (gAlice->TreeE())//check if tree E exists
180      {
181       Nevents = (Int_t)gAlice->TreeE()->GetEntries();//if yes get number of events in gAlice
182       cout<<"________________________________________________________\n";
183       cout<<"Found "<<Nevents<<" event(s) in directory "<<GetDirName(currentdir)<<endl;
184       cout<<"Setting Magnetic Field: B="<<gAlice->Field()->SolenoidField()<<"T"<<endl;
185       AliKalmanTrack::SetConvConst(1000/0.299792458/gAlice->Field()->SolenoidField());
186      }
187     else
188      {//if not return an error
189        Error("Read","Can not find Header tree (TreeE) in gAlice");
190        currentdir++;
191        continue;
192      }
193   
194     aClustersFile->cd();//set cluster file active 
195     AliTPCParam *TPCParam= (AliTPCParam*)aClustersFile->Get("75x40_100x60_150x60");
196     if (!TPCParam) 
197       { 
198        TPCParam= (AliTPCParam*)aClustersFile->Get("75x40_100x60");
199        if (!TPCParam) 
200         { 
201           Error("Read","TPC parameters have not been found !\n");
202           currentdir++;
203           continue;
204         }
205       }
206
207   
208     for(Int_t currentEvent =0; currentEvent<Nevents;currentEvent++)//loop over all events
209      {
210        cout<<"Reading Event "<<currentEvent<<endl;
211        /**************************************/
212         /**************************************/
213          /**************************************/ 
214          
215          aTracksFile->cd();//set track file active
216           
217          Char_t  treename[100];
218          sprintf(treename,"TreeT_TPC_%d",currentEvent);//prepare name of the tree
219    
220          TTree *tracktree=0;
221          
222          tracktree=(TTree*)aTracksFile->Get(treename);//get the tree 
223          if (!tracktree) //check if we got the tree
224           {//if not return with error
225             Error("Read","Can't get a tree with TPC tracks !\n"); 
226             continue;
227           }
228    
229          TBranch *trackbranch=tracktree->GetBranch("tracks");//get the branch with tracks
230          if (!trackbranch) ////check if we got the branch
231           {//if not return with error
232             Error("Read","Can't get a branch with TPC tracks !\n"); 
233             continue;
234           }
235          Int_t NTPCtracks=(Int_t)tracktree->GetEntries();//get number of TPC tracks 
236          cout<<"Found "<<NTPCtracks<<" TPC tracks.\n";
237          //Copy tracks to array
238          
239          AliTPCtrack *iotrack=0;
240          
241          aClustersFile->cd();//set cluster file active 
242          AliTPCtracker *tracker = new AliTPCtracker(TPCParam,currentEvent);//create the tacker for this event
243          if (!tracker) //check if it has created succeffuly
244           {//if not return with error
245             Error("Read","Can't get a tracker !\n"); 
246             continue;
247           }
248          tracker->LoadInnerSectors();
249          tracker->LoadOuterSectors();
250    
251          for (i=0; i<NTPCtracks; i++) //loop over all tpc tracks
252           {
253             iotrack=new AliTPCtrack;   //create new tracks
254             trackbranch->SetAddress(&iotrack); //tell the branch ehere to put track data from tree(file)
255             tracktree->GetEvent(i); //stream track i to the iotrack
256             tracker->CookLabel(iotrack,0.1); //calculate (cook) the label of the tpc track
257                                              //which is the label of corresponding simulated particle 
258             tarray->AddLast(iotrack); //put the track in the array
259           }
260          
261          aTracksFile->Delete(treename);//delete tree from memmory (and leave untached on disk)- we do not need it any more
262          aTracksFile->Delete("tracks");//delete branch from memmory
263          delete tracker; //delete tracker
264          
265          tracker = 0x0;
266          trackbranch = 0x0;
267          tracktree = 0x0;
268    
269          Double_t xk;
270          Double_t par[5];
271          Float_t phi, lam, pt;//angles and transverse momentum
272          Int_t label; //label of the current track
273          
274          aGAliceFile->cd();
275          gAlice->GetEvent(currentEvent); 
276
277          gAlice->Particles();
278          
279          for (i=0; i<NTPCtracks; i++) //loop over all good tracks
280           { 
281             iotrack = (AliTPCtrack*)tarray->At(i);
282             label = iotrack->GetLabel();
283
284             if (label < 0) continue;
285             
286             TParticle *p = (TParticle*)gAlice->Particle(label);
287             
288             if(p == 0x0) continue; //if returned pointer is NULL
289             if(p->GetPDG() == 0x0) continue; //if particle has crezy PDG code (not known to our database)
290             
291             if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type 
292                                         //if not take next partilce
293             
294             AliHBTParticle* part = new AliHBTParticle(*p);
295             if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
296                                                     //if it does not delete it and take next good track
297          
298             iotrack->PropagateToVertex();
299
300             iotrack->GetExternalParameters(xk,par);     //get properties of the track
301             phi=TMath::ASin(par[2]) + iotrack->GetAlpha(); 
302             if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
303             if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
304             lam=par[3]; 
305             pt=1.0/TMath::Abs(par[4]);
306             
307             Double_t tpx = pt * TMath::Cos(phi); //track x coordinate of momentum
308             Double_t tpy = pt * TMath::Sin(phi); //track y coordinate of momentum
309             Double_t tpz = pt * lam; //track z coordinate of momentum
310             
311             Double_t mass = p->GetMass();
312             Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);//total energy of the track
313             
314             AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(), tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
315             if(Pass(track))//check if meets all criteria of any of our cuts
316                          //if it does not delete it and take next good track
317              { 
318                delete track;
319                delete part;
320                continue;
321              }
322             particles->AddParticle(totalNevents,part);//put track and particle on the run
323             tracks->AddParticle(totalNevents,track);
324
325           }
326          tarray->Clear(); //clear the array
327          
328         /**************************************/
329        /**************************************/
330       /**************************************/  
331      totalNevents++;
332     }
333   
334     //save environment (resouces) --
335     //clean your place after the work
336     CloseFiles(aTracksFile,aClustersFile,aGAliceFile); 
337     currentdir++;
338    }while(currentdir < Ndirs);
339
340
341   delete tarray;
342   fIsRead = kTRUE;
343   return 0;
344  }
345
346 /********************************************************************/
347 Int_t AliHBTReaderTPC::OpenFiles
348 (TFile*& aTracksFile, TFile*& aClustersFile, TFile*& agAliceFile,Int_t event)
349 {
350  //opens all the files
351    
352    
353    const TString& dirname = GetDirName(event); 
354    if (dirname == "")
355     {
356       Error("OpenFiles","Can not get directory name");
357       return 4;
358     }
359    
360    TString filename = dirname +"/"+ fTrackFileName;
361    aTracksFile = TFile::Open(filename.Data());
362    if ( aTracksFile  == 0x0 ) 
363      {
364        Error("OpenFiles","Can't open file with tacks named %s",filename.Data());
365        return 1;
366      }
367    if (!aTracksFile->IsOpen())
368      {
369        Error("OpenFiles","Can't open file with tacks named %s",filename.Data());
370        return 1;
371      }
372   
373    filename = dirname +"/"+ fClusterFileName;
374    aClustersFile = TFile::Open(filename.Data());
375    if ( aClustersFile == 0x0 )
376     {
377       Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
378       return 2;
379     }
380    if (!aClustersFile->IsOpen())
381     {
382       Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
383       return 2;
384     }
385
386    filename = dirname +"/"+ fGAliceFileName;
387    agAliceFile = TFile::Open(filename.Data());
388    if ( agAliceFile== 0x0)
389     {
390       Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
391       return 3;
392     }
393    if (!agAliceFile->IsOpen())
394     {
395       Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
396       return 3;
397     } 
398    
399    if (!(gAlice=(AliRun*)agAliceFile->Get("gAlice"))) 
400     {
401       Error("OpenFiles","gAlice have not been found on %s !\n",filename.Data());
402       return 5;
403     }
404
405    return 0; 
406 }
407 /********************************************************************/
408
409 /********************************************************************/
410   
411 void AliHBTReaderTPC::CloseFiles(TFile*& tracksFile, TFile*& clustersFile, TFile*& gAliceFile)
412 {
413   //closes the files
414   tracksFile->Close();
415   delete tracksFile;
416   tracksFile = 0x0;
417   clustersFile->Close();
418   delete clustersFile;
419   clustersFile = 0x0;
420
421   delete gAlice;
422   gAlice = 0;
423
424   if (gAliceFile) 
425    {
426      gAliceFile->Close();
427      delete gAliceFile;
428      gAliceFile = 0x0;
429    }
430 }
431
432 /********************************************************************/
433
434 /********************************************************************/
435 /********************************************************************/
436 /********************************************************************/
437