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