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