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