HBTAN/AliHBTReaderTPC.cxx

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