Comilation Warning removal
[u/mrichter/AliRoot.git] / HBTAN / AliHBTReaderTPC.cxx
1 #include "AliHBTReaderTPC.h"
2 //______________________________________________
3 //
4 // class AliHBTReaderTPC
5 //
6 // reader for TPC tracking
7 // needs galice.root, AliTPCtracks.root, AliTPCclusters.root
8 //
9 // more info: http://aliweb.cern.ch/people/skowron/analyzer/index.html
10 // Piotr.Skowronski@cern.ch
11 //
12 ///////////////////////////////////////////////////////////////////////////
13 #include <TTree.h>
14 #include <TFile.h>
15 #include <TParticle.h>
16 #include <TH1.h>
17
18 #include <Varargs.h>
19
20 #include <AliRun.h>
21 #include <AliLoader.h>
22 #include <AliStack.h>
23 #include <AliMagF.h>
24 #include <AliTPCtrack.h>
25 #include <AliTPCParam.h>
26 #include <AliTPCtracker.h>
27 #include <AliTPCLoader.h>
28
29 #include "AliHBTRun.h"
30 #include "AliHBTEvent.h"
31 #include "AliHBTParticle.h"
32 #include "AliHBTParticleCut.h"
33
34
35
36 ClassImp(AliHBTReaderTPC)
37
38 AliHBTReaderTPC::AliHBTReaderTPC():
39  fFileName("galice.root"),
40  fRunLoader(0x0),
41  fTPCLoader(0x0),
42  fMagneticField(0.0),
43  fUseMagFFromRun(kTRUE),
44  fNClustMin(0),
45  fNClustMax(190),
46  fNChi2PerClustMin(0.0),
47  fNChi2PerClustMax(10e5),
48  fC44Min(0.0),
49  fC44Max(10e5)
50 {
51   //constructor, 
52   //Defaults:
53   //  galicefilename = ""  - this means: Do not open gAlice file - 
54   //                         just leave the global pointer untouched
55   
56 }
57
58 AliHBTReaderTPC::AliHBTReaderTPC(const Char_t* galicefilename):
59  fFileName(galicefilename),
60  fRunLoader(0x0),
61  fTPCLoader(0x0),
62  fMagneticField(0.0),
63  fUseMagFFromRun(kTRUE),
64  fNClustMin(0),
65  fNClustMax(190),
66  fNChi2PerClustMin(0.0),
67  fNChi2PerClustMax(10e5),
68  fC44Min(0.0),
69  fC44Max(10e5)
70 {
71   //constructor, 
72   //Defaults:
73   //  galicefilename = ""  - this means: Do not open gAlice file - 
74   //                         just leave the global pointer untouched
75   
76 }
77 /********************************************************************/
78 AliHBTReaderTPC::AliHBTReaderTPC(TObjArray* dirs, const Char_t* galicefilename):
79  AliHBTReader(dirs), 
80  fFileName(galicefilename),
81  fRunLoader(0x0),
82  fTPCLoader(0x0),
83  fMagneticField(0.0),
84  fUseMagFFromRun(kTRUE),
85  fNClustMin(0),
86  fNClustMax(190),
87  fNChi2PerClustMin(0.0),
88  fNChi2PerClustMax(10e5),
89  fC44Min(0.0),
90  fC44Max(10e5)
91 {
92   //constructor, 
93   //Defaults:
94   //  galicefilename = ""  - this means: Do not open gAlice file - 
95   //                         just leave the global pointer untached
96   
97 }
98 /********************************************************************/
99
100 AliHBTReaderTPC::~AliHBTReaderTPC()
101 {
102  //desctructor
103    if (AliHBTParticle::GetDebug()) 
104     {
105       Info("~AliHBTReaderTPC","deleting Run Loader");
106       AliLoader::SetDebug(AliHBTParticle::GetDebug());
107     }
108    
109    delete fRunLoader;
110    
111    if (AliHBTParticle::GetDebug()) 
112     {
113       Info("~AliHBTReaderTPC","deleting Run Loader Done");
114     }
115 }
116 /********************************************************************/
117  
118 void AliHBTReaderTPC::Rewind()
119 {
120   delete fRunLoader;
121   fRunLoader = 0x0;
122   fCurrentDir = 0;
123   fNEventsRead= 0;
124   if (fTrackCounter) fTrackCounter->Reset();
125 }
126 /********************************************************************/
127
128 Int_t AliHBTReaderTPC::ReadNext()
129  {
130  //reads data and puts put to the particles and tracks objects
131  //reurns 0 if everything is OK
132  //
133   Info("Read","");
134   
135   TObjArray *tarray = new TObjArray(5000); //cotainer for tpc tracks
136   tarray->SetOwner(); //set the ownership of the objects it contains
137                       //when array is is deleted or cleared all objects 
138                       //that it contains are deleted
139
140   if (fParticlesEvent == 0x0)  fParticlesEvent = new AliHBTEvent();
141   if (fTracksEvent == 0x0)  fTracksEvent = new AliHBTEvent();
142
143   fParticlesEvent->Reset();
144   fTracksEvent->Reset();
145
146   do  //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
147    {
148       
149     if (fRunLoader == 0x0) 
150       if (OpenNextSession()) continue;//directory counter is increased inside in case of error
151
152     if (fCurrentEvent == fRunLoader->GetNumberOfEvents())
153      {
154        //read next directory
155        delete fRunLoader;//close current session
156        fRunLoader = 0x0;//assure pointer is null
157        fCurrentDir++;//go to next dir
158        continue;//directory counter is increased inside in case of error
159      }
160      
161     Info("ReadNext","Reading Event %d",fCurrentEvent);
162     
163     fRunLoader->GetEvent(fCurrentEvent);
164     TTree *tracktree = fTPCLoader->TreeT();//get the tree 
165     if (!tracktree) //check if we got the tree
166       {//if not return with error
167          Error("ReadNext","Can't get a tree with TPC tracks !\n"); 
168          fCurrentEvent++;//go to next dir
169          continue;
170       }
171    
172     TBranch *trackbranch=tracktree->GetBranch("tracks");//get the branch with tracks
173     if (!trackbranch) ////check if we got the branch
174       {//if not return with error
175         Error("ReadNext","Can't get a branch with TPC tracks !\n"); 
176         fCurrentEvent++;//go to next dir
177         continue;
178       }
179     Int_t ntpctracks=(Int_t)tracktree->GetEntries();//get number of TPC tracks 
180     Info("ReadNext","Found %d TPC tracks.",ntpctracks);
181     //Copy tracks to array
182
183     AliTPCtrack *iotrack=0;
184
185     for (Int_t i=0; i<ntpctracks; i++) //loop over all tpc tracks
186      {
187        iotrack=new AliTPCtrack;   //create new tracks
188        trackbranch->SetAddress(&iotrack); //tell the branch ehere to put track data from tree(file)
189        tracktree->GetEvent(i); //stream track i to the iotrack
190        tarray->AddLast(iotrack); //put the track in the array
191      }
192
193     Double_t xk;
194     Double_t par[5];
195     Float_t phi, lam, pt;//angles and transverse momentum
196     Int_t label; //label of the current track
197
198     AliStack* stack = fRunLoader->Stack();
199     if (stack == 0x0)
200      {
201        Error("ReadNext","Can not get stack for current event",fCurrentEvent);
202        fCurrentEvent++;
203        continue;
204      }
205     stack->Particles();
206
207     for (Int_t i=0; i<ntpctracks; i++) //loop over all good tracks
208      { 
209        iotrack = (AliTPCtrack*)tarray->At(i);
210        label = iotrack->GetLabel();
211
212        if (label < 0) continue;
213        
214        if (CheckTrack(iotrack)) continue;
215        
216        TParticle *p = (TParticle*)stack->Particle(label);
217
218        if(p == 0x0) continue; //if returned pointer is NULL
219        if(p->GetPDG() == 0x0) continue; //if particle has crezy PDG code (not known to our database)
220
221        if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type 
222                                    //if not take next partilce
223
224        AliHBTParticle* part = new AliHBTParticle(*p,i);
225        if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
226                                                //if it does not delete it and take next good track
227
228 //       iotrack->PropagateTo(3.,0.0028,65.19);
229 //       iotrack->PropagateToVertex(36.66,1.2e-3);//orig
230        iotrack->PropagateToVertex(50.,0.0353);
231        
232        iotrack->GetExternalParameters(xk,par);     //get properties of the track
233        
234        phi=TMath::ASin(par[2]) + iotrack->GetAlpha(); 
235        if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
236        if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
237        lam=par[3]; 
238        pt=1.0/TMath::Abs(par[4]);
239
240        Double_t tpx = pt * TMath::Cos(phi); //track x coordinate of momentum
241        Double_t tpy = pt * TMath::Sin(phi); //track y coordinate of momentum
242        Double_t tpz = pt * lam; //track z coordinate of momentum
243
244        Double_t mass = p->GetMass();
245        Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);//total energy of the track
246
247        AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(),i, tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
248        if(Pass(track))//check if meets all criteria of any of our cuts
249                     //if it does not delete it and take next good track
250         { 
251           delete track;
252           delete part;
253           continue;
254         }
255        fParticlesEvent->AddParticle(part);
256        fTracksEvent->AddParticle(track);
257       }
258       
259     Info("ReadNext","Read %d tracks and %d particles from event %d (event %d in dir %d).",
260             fParticlesEvent->GetNumberOfParticles(), fTracksEvent->GetNumberOfParticles(),
261             fNEventsRead,fCurrentEvent,fCurrentDir);
262     fTrackCounter->Fill(fTracksEvent->GetNumberOfParticles());
263     fCurrentEvent++;
264     fNEventsRead++;
265     delete tarray;
266     return 0;
267    }while(fCurrentDir < GetNumberOfDirs());
268
269   delete tarray;
270   return 1;
271  }
272 /********************************************************************/
273
274 Int_t AliHBTReaderTPC::OpenNextSession()
275 {
276   TString filename = GetDirName(fCurrentDir);
277   if (filename.IsNull())
278    {
279      DoOpenError("Can not get directory name");
280      return 1;
281    }
282   filename = filename +"/"+ fFileName;
283   fRunLoader = AliRunLoader::Open(filename,AliConfig::fgkDefaultEventFolderName);
284   if( fRunLoader == 0x0)
285    {
286      DoOpenError("Can not open session.");
287      return 1;
288    }
289
290   fRunLoader->LoadHeader();
291   if (fRunLoader->GetNumberOfEvents() <= 0)
292    {
293      DoOpenError("There is no events in this directory.");
294      return 1;
295    }
296
297   if (fRunLoader->LoadKinematics())
298    {
299      DoOpenError("Error occured while loading kinematics.");
300      return 1;
301    }
302   fTPCLoader = (AliTPCLoader*)fRunLoader->GetLoader("TPCLoader");
303   if ( fTPCLoader == 0x0)
304    {
305      DoOpenError("Exiting due to problems with opening files.");
306      return 1;
307    }
308   Info("OpenNextSession","________________________________________________________");
309   Info("OpenNextSession","Found %d event(s) in directory %s",
310         fRunLoader->GetNumberOfEvents(),GetDirName(fCurrentDir).Data());
311   Float_t mf;
312   if (fUseMagFFromRun)
313    {
314      fRunLoader->LoadgAlice();
315      mf = fRunLoader->GetAliRun()->Field()->SolenoidField();
316      Info("OpenNextSession","Setting Magnetic Field from run: B=%fT",mf/10.);
317      fRunLoader->UnloadgAlice();
318    }
319   else
320    {
321      Info("OpenNextSession","Setting Own Magnetic Field: B=%fT",fMagneticField);
322      if (fMagneticField == 0x0)
323       {
324         Fatal("OpenNextSession","Magnetic field can not be 0.");
325         return 1;//pro forma
326       }
327      mf = fMagneticField*10.;
328    }
329   AliKalmanTrack::SetConvConst(1000/0.299792458/mf);
330
331
332   fRunLoader->CdGAFile();
333   AliTPCParam *TPCParam= (AliTPCParam*)gDirectory->Get("75x40_100x60");
334   if (!TPCParam) 
335    {
336     TPCParam=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
337     if (!TPCParam) 
338      { 
339        DoOpenError("TPC parameters have not been found !\n");
340        return 1;
341      }
342    }
343
344   if (fTPCLoader->LoadTracks())
345    {
346      DoOpenError("Error occured while loading TPC tracks.");
347      return 1;
348    }
349   
350   fCurrentEvent = 0;
351   return 0;
352 }
353 /********************************************************************/
354
355 void AliHBTReaderTPC::DoOpenError( const char *va_(fmt), ...)
356 {
357   // Does error display and clean-up in case error caught on Open Next Session
358
359    va_list ap;
360    va_start(ap,va_(fmt));
361    Error("OpenNextSession", va_(fmt), ap);
362    va_end(ap);
363    
364    delete fRunLoader;
365    fRunLoader = 0x0;
366    fTPCLoader = 0x0;
367    fCurrentDir++;
368 }
369
370 /********************************************************************/
371 Bool_t AliHBTReaderTPC::CheckTrack(AliTPCtrack* t)
372 {
373   //Performs check of the track
374   if ( (t->GetNumberOfClusters() > fNClustMax) || (t->GetNumberOfClusters() < fNClustMin) ) return kTRUE;
375
376   Double_t cc[15];
377   t->GetCovariance(cc);
378   if ( (cc[9] < fC44Min) || (cc[9] > fC44Max) ) return kTRUE;
379   
380   Float_t chisqpercl = t->GetChi2()/((Double_t)t->GetNumberOfClusters());
381   if ( (chisqpercl < fNChi2PerClustMin) || (chisqpercl > fNChi2PerClustMax) ) return kTRUE;
382   
383   return kFALSE;
384   
385 }
386 /********************************************************************/
387
388 void AliHBTReaderTPC::SetNClustersRange(Int_t min,Int_t max)
389 {
390  //sets range of Number Of Clusters that tracks have to have
391  fNClustMin = min;
392  fNClustMax = max;
393 }
394 /********************************************************************/
395
396 void AliHBTReaderTPC::SetChi2PerCluserRange(Float_t min, Float_t max)
397 {
398   //sets range of Chi2 per Cluster
399   fNChi2PerClustMin = min;
400   fNChi2PerClustMax = max;
401 }
402 /********************************************************************/
403
404 void AliHBTReaderTPC::SetC44Range(Float_t min, Float_t max)
405 {
406  //Sets range of C44 parameter of covariance matrix of the track
407  //it defines uncertainty of the momentum
408   fC44Min = min;
409   fC44Max = max;
410 }
411
412 /********************************************************************/
413 /********************************************************************/
414 /********************************************************************/
415
416 /*
417 void (AliTPCtrack* iotrack, Double_t curv)
418 {
419   Double_t x[5];
420   iotrack->GetExternalPara
421   //xk is a 
422   Double_t fP4=iotrack->GetC();
423   Double_t fP2=iotrack->GetEta();
424   
425   Double_t x1=fX, x2=x1+(xk-x1), dx=x2-x1, y1=fP0, z1=fP1;
426   Double_t c1=fP4*x1 - fP2, r1=sqrt(1.- c1*c1);
427   Double_t c2=fP4*x2 - fP2, r2=sqrt(1.- c2*c2);
428   
429   fP0 += dx*(c1+c2)/(r1+r2);
430   fP1 += dx*(c1+c2)/(c1*r2 + c2*r1)*fP3;
431
432 }
433 */
434