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