Retrofeeded recent developement in from AliHBTReaderESD (ITS track points, additional...
[u/mrichter/AliRoot.git] / HBTAN / AliHBTReaderITSv2.cxx
1 #include "AliHBTReaderITSv2.h"
2
3 #include <Varargs.h>
4
5 #include <TString.h>
6 #include <TObjString.h>
7 #include <TTree.h>
8 #include <TParticle.h>
9
10 #include <AliRun.h>
11 #include <AliRunLoader.h>
12 #include <AliLoader.h>
13 #include <AliStack.h>
14 #include <AliMagF.h>
15 #include <AliITStrackV2.h>
16
17 #include "AliHBTRun.h"
18 #include "AliHBTEvent.h"
19 #include "AliHBTParticle.h"
20 #include "AliHBTParticleCut.h"
21
22
23 ClassImp(AliHBTReaderITSv2)
24
25 AliHBTReaderITSv2::AliHBTReaderITSv2():
26  fFileName("galice.root"),
27  fRunLoader(0x0),
28  fITSLoader(0x0),
29  fMagneticField(0.0),
30  fUseMagFFromRun(kTRUE)
31 {
32   //constructor, 
33   //Defaults:
34   //  galicefilename = "galice.root"
35 }
36 /********************************************************************/
37
38 AliHBTReaderITSv2::AliHBTReaderITSv2(const Char_t* galicefilename):
39  fFileName(galicefilename),
40  fRunLoader(0x0),
41  fITSLoader(0x0),
42  fMagneticField(0.0),
43  fUseMagFFromRun(kTRUE)
44 {
45   //constructor, 
46   //Defaults:
47   //  galicefilename = "galice.root"
48 }
49 /********************************************************************/
50
51 AliHBTReaderITSv2::AliHBTReaderITSv2(TObjArray* dirs, const Char_t* galicefilename): 
52  AliHBTReader(dirs),
53  fFileName(galicefilename),
54  fRunLoader(0x0),
55  fITSLoader(0x0),
56  fMagneticField(0.0),
57  fUseMagFFromRun(kTRUE)
58 {
59   //constructor, 
60   //Defaults:
61   //  galicefilename = "galice.root"
62   
63 }
64 /********************************************************************/
65  
66 void AliHBTReaderITSv2::Rewind()
67 {
68   //rewinds reading
69   delete fRunLoader;
70   fRunLoader = 0x0;
71   fCurrentDir = 0;
72   fNEventsRead= 0;
73 }
74 /********************************************************************/
75
76 AliHBTReaderITSv2::~AliHBTReaderITSv2()
77 {
78   //dtor
79   delete fRunLoader;
80 }
81 /********************************************************************/
82  
83 Int_t AliHBTReaderITSv2::ReadNext()
84 {
85 //reads data from next event
86   
87  register Int_t i = 0; //iterator
88  
89 // AliITStrackerV2 *tracker; // ITS tracker - used for cooking labels
90  TTree *tracktree = 0x0; // tree for tracks
91  AliITStrackV2 *iotrack = 0x0;
92  
93  Double_t xk;
94  Double_t par[5]; //Kalman track parameters
95  Float_t phi, lam, pt;//angles and transverse momentum
96  Int_t label; //label of the current track
97  
98  
99  if (fParticlesEvent == 0x0)  fParticlesEvent = new AliHBTEvent();
100  if (fTracksEvent == 0x0)  fTracksEvent = new AliHBTEvent();
101
102  fParticlesEvent->Reset();
103  fTracksEvent->Reset();
104  do //do while is good even if Ndirs==0 (than read from current directory)
105    {
106     if (fRunLoader == 0x0) 
107       if (OpenNextFile()) continue;//directory counter is increased inside in case of error
108
109     if (fCurrentEvent == fRunLoader->GetNumberOfEvents())
110      {
111        //read next directory
112        delete fRunLoader;//close current session
113        fRunLoader = 0x0;//assure pointer is null
114        fCurrentDir++;//go to next dir
115        continue;//directory counter is increased inside in case of error
116      }
117      
118     Info("ReadNext","Reading Event %d",fCurrentEvent);
119      
120     fRunLoader->GetEvent(fCurrentEvent);
121     
122     tracktree=fITSLoader->TreeT();
123     if (!tracktree) 
124      {
125        Error("ReadNext","Can't get a tree with ITS tracks"); 
126        fCurrentEvent++;
127        continue;
128      }
129       
130     TBranch *tbranch=tracktree->GetBranch("tracks");
131     if (!tbranch) 
132      {
133        Error("ReadNext","Can't get a branch with ITS tracks"); 
134        fCurrentEvent++;
135        continue;
136      }
137
138     AliStack* stack = fRunLoader->Stack();
139     if (stack == 0x0)
140      {
141        Error("ReadNext","Can not get stack for current event",fCurrentEvent);
142        fCurrentEvent++;
143        continue;
144      }
145      
146     //must be here because on the beginning conv. const. is not set yet 
147     if (iotrack == 0x0) iotrack = new AliITStrackV2(); //buffer track for reading data from tree
148     
149     Int_t ntr = (Int_t)tracktree->GetEntries();
150     
151     for (i=0; i < ntr; i++) //loop over all tpc tracks
152      { 
153        tbranch->SetAddress(&iotrack);
154        tracktree->GetEvent(i);
155
156        label=iotrack->GetLabel();
157        if (label < 0) 
158         {
159           continue;
160         }
161
162        TParticle *p = stack->Particle(label);
163        if(p == 0x0) continue; //if returned pointer is NULL
164        if(p->GetPDG() == 0x0) continue; //if particle has crezy PDG code (not known to our database)
165
166        if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type 
167                                            //if not take next partilce
168
169        AliHBTParticle* part = new AliHBTParticle(*p,i);
170        if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
171                                                //if it does not delete it and take next good track
172
173        iotrack->PropagateTo(3.,0.0028,65.19);
174        iotrack->PropagateToVertex();
175
176        iotrack->GetExternalParameters(xk,par);     //get properties of the track
177        phi=TMath::ASin(par[2]) + iotrack->GetAlpha(); 
178        if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
179        if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
180        lam=par[3]; 
181        pt=1.0/TMath::Abs(par[4]);
182
183        Double_t tpx = pt * TMath::Cos(phi); //track x coordinate of momentum
184        Double_t tpy = pt * TMath::Sin(phi); //track y coordinate of momentum
185        Double_t tpz = pt * lam; //track z coordinate of momentum
186
187        Double_t mass = p->GetMass();
188        Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);//total energy of the track
189
190        AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(), i, tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
191        if(Pass(track))//check if meets all criteria of any of our cuts
192                       //if it does not delete it and take next good track
193         { 
194          delete track;
195          delete part;
196          continue;
197         }
198         
199        fParticlesEvent->AddParticle(part);
200        fTracksEvent->AddParticle(track);
201      }//end of loop over tracks in the event
202        
203     Info("ReadNext","Read %d tracks and %d particles from event %d (event %d in dir %d).",
204             fParticlesEvent->GetNumberOfParticles(), fTracksEvent->GetNumberOfParticles(),
205             fNEventsRead,fCurrentEvent,fCurrentDir);
206      
207     fCurrentEvent++;
208     fNEventsRead++;
209     delete iotrack;
210     return 0;
211    }while(fCurrentDir < GetNumberOfDirs());//end of loop over directories specified in fDirs Obj Array
212
213  delete iotrack;
214  return 1;
215 }
216
217 /********************************************************************/
218 Int_t AliHBTReaderITSv2::OpenNextFile()
219 {
220   //opens next file
221   TString filename = GetDirName(fCurrentDir);
222   if (filename.IsNull())
223    {
224      DoOpenError("Can not get directory name");
225      return 1;
226    }
227   filename = filename +"/"+ fFileName;
228   fRunLoader = AliRunLoader::Open(filename,AliConfig::GetDefaultEventFolderName());
229   if( fRunLoader == 0x0)
230    {
231      DoOpenError("Can not open session.");
232      return 1;
233    }
234
235   if (fRunLoader->GetNumberOfEvents() <= 0)
236    {
237      DoOpenError("There is no events in this directory.");
238      return 1;
239    }
240
241   if (fRunLoader->LoadKinematics())
242    {
243      DoOpenError("Error occured while loading kinematics.");
244      return 1;
245    }
246   fITSLoader = fRunLoader->GetLoader("ITSLoader");
247   if ( fITSLoader == 0x0)
248    {
249      DoOpenError("Exiting due to problems with opening files.");
250      return 1;
251    }
252    
253   Info("OpenNextSession","________________________________________________________");
254   Info("OpenNextSession","Found %d event(s) in directory %s",
255         fRunLoader->GetNumberOfEvents(),GetDirName(fCurrentDir).Data());
256   Float_t mf;
257   if (fUseMagFFromRun)
258    {
259      if (fRunLoader->LoadgAlice())
260       {
261         DoOpenError("Error occured while loading AliRun.");
262         return 1;
263       }
264      mf = fRunLoader->GetAliRun()->Field()->SolenoidField();
265      Info("OpenNextSession","Setting Magnetic Field from run: B=%fT",mf/10.);
266      fRunLoader->UnloadgAlice();
267    }
268   else
269    {
270      Info("OpenNextSession","Setting Own Magnetic Field: B=%fT",fMagneticField);
271      if (fMagneticField == 0x0)
272       {
273         Fatal("OpenNextSession","Magnetic field can not be 0.");
274         return 1;//pro forma
275       }
276      mf = fMagneticField*10.;
277    }
278   AliKalmanTrack::SetConvConst(1000/0.299792458/mf);
279
280   if (fITSLoader->LoadTracks())
281    {
282      DoOpenError("Error occured while loading TPC tracks.");
283      return 1;
284    }
285   
286   fCurrentEvent = 0;
287   return 0;
288 }
289 /********************************************************************/
290
291 void AliHBTReaderITSv2::DoOpenError( const char *va_(fmt), ...)
292 {
293   // Does error display and clean-up in case error caught on Open Next Session
294
295    va_list ap;
296    va_start(ap,va_(fmt));
297    Error("OpenNextFile", va_(fmt), ap);
298    va_end(ap);
299    
300    delete fRunLoader;
301    fRunLoader = 0x0;
302    fITSLoader = 0x0;
303    fCurrentDir++;
304 }
305
306 /********************************************************************/