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