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   Int_t ntracksread = 0;
 136
 137   TObjArray *tarray = new TObjArray(5000); //cotainer for tpc tracks
 138   tarray->SetOwner(); //set the ownership of the objects it contains
 139                       //when array is is deleted or cleared all objects
 140                       //that it contains are deleted
 141
 142   if (fParticlesEvent == 0x0)  fParticlesEvent = new AliHBTEvent();
 143   if (fTracksEvent == 0x0)  fTracksEvent = new AliHBTEvent();
 144
 145   fParticlesEvent->Reset();
 146   fTracksEvent->Reset();
 147
 148   do  //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
 149    {
 150
 151     if (fRunLoader == 0x0)
 152       if (OpenNextSession()) continue;//directory counter is increased inside in case of error
 153
 154     if (fCurrentEvent == fRunLoader->GetNumberOfEvents())
 155      {
 156        //read next directory
 157        delete fRunLoader;//close current session
 158        fRunLoader = 0x0;//assure pointer is null
 159        fCurrentDir++;//go to next dir
 160        continue;//directory counter is increased inside in case of error
 161      }
 162
 163     Info("ReadNext","Reading Event %d",fCurrentEvent);
 164
 165     fRunLoader->GetEvent(fCurrentEvent);
 166     TTree *tracktree = fTPCLoader->TreeT();//get the tree
 167     if (!tracktree) //check if we got the tree
 168       {//if not return with error
 169          Error("ReadNext","Can't get a tree with TPC tracks !\n");
 170          fCurrentEvent++;//go to next dir
 171          continue;
 172       }
 173
 174     TBranch *trackbranch=tracktree->GetBranch("tracks");//get the branch with tracks
 175     if (!trackbranch) ////check if we got the branch
 176       {//if not return with error
 177         Error("ReadNext","Can't get a branch with TPC tracks !\n");
 178         fCurrentEvent++;//go to next dir
 179         continue;
 180       }
 181     Int_t ntpctracks=(Int_t)tracktree->GetEntries();//get number of TPC tracks
 182     Info("ReadNext","Found %d TPC tracks.",ntpctracks);
 183     //Copy tracks to array
 184
 185     AliTPCtrack *iotrack=0;
 186
 187     for (Int_t i=0; i<ntpctracks; i++) //loop over all tpc tracks
 188      {
 189        iotrack=new AliTPCtrack;   //create new tracks
 190        trackbranch->SetAddress(&iotrack); //tell the branch ehere to put track data from tree(file)
 191        tracktree->GetEvent(i); //stream track i to the iotrack
 192        tarray->AddLast(iotrack); //put the track in the array
 193      }
 194
 195     Double_t xk;
 196     Double_t par[5];
 197     Float_t phi, lam, pt;//angles and transverse momentum
 198     Int_t label; //label of the current track
 199
 200     AliStack* stack = fRunLoader->Stack();
 201     if (stack == 0x0)
 202      {
 203        Error("ReadNext","Can not get stack for current event",fCurrentEvent);
 204        fCurrentEvent++;
 205        continue;
 206      }
 207     stack->Particles();
 208
 209     for (Int_t i=0; i<ntpctracks; i++) //loop over all good tracks
 210      {
 211        iotrack = (AliTPCtrack*)tarray->At(i);
 212        label = iotrack->GetLabel();
 213
 214        if (label < 0) continue;
 215
 216        if (CheckTrack(iotrack)) continue;
 217
 218        TParticle *p = (TParticle*)stack->Particle(label);
 219
 220        if(p == 0x0) continue; //if returned pointer is NULL
 221        if(p->GetPDG() == 0x0) continue; //if particle has crezy PDG code (not known to our database)
 222
 223        if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type
 224                                    //if not take next partilce
 225
 226        AliHBTParticle* part = new AliHBTParticle(*p,i);
 227        if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
 228                                                //if it does not delete it and take next good track
 229
 230 //       iotrack->PropagateTo(3.,0.0028,65.19);
 231 //       iotrack->PropagateToVertex(36.66,1.2e-3);//orig
 232        iotrack->PropagateToVertex(50.,0.0353);
 233
 234        iotrack->GetExternalParameters(xk,par);     //get properties of the track
 235
 236        phi=TMath::ASin(par[2]) + iotrack->GetAlpha();
 237        if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
 238        if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
 239        lam=par[3];
 240        pt=1.0/TMath::Abs(par[4]);
 241
 242        Double_t tpx = pt * TMath::Cos(phi); //track x coordinate of momentum
 243        Double_t tpy = pt * TMath::Sin(phi); //track y coordinate of momentum
 244        Double_t tpz = pt * lam; //track z coordinate of momentum
 245
 246        Double_t mass = p->GetMass();
 247        Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);//total energy of the track
 248
 249        AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(),i, tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
 250        if(Pass(track))//check if meets all criteria of any of our cuts
 251                     //if it does not delete it and take next good track
 252         {
 253           delete track;
 254           delete part;
 255           continue;
 256         }
 257        fParticlesEvent->AddParticle(part);
 258        fTracksEvent->AddParticle(track);
 259       }
 260
 261     Info("ReadNext","Read %d tracks and %d particles from event %d (event %d in dir %d).",
 262             fParticlesEvent->GetNumberOfParticles(), fTracksEvent->GetNumberOfParticles(),
 263             fNEventsRead,fCurrentEvent,fCurrentDir);
 264     fTrackCounter->Fill(fTracksEvent->GetNumberOfParticles());
 265     fCurrentEvent++;
 266     fNEventsRead++;
 267     delete tarray;
 268     return 0;
 269    }while(fCurrentDir < GetNumberOfDirs());
 270
 271   delete tarray;
 272   return 1;
 273  }
 274 /********************************************************************/
 275
 276 Int_t AliHBTReaderTPC::OpenNextSession()
 277 {
 278   TString filename = GetDirName(fCurrentDir);
 279   if (filename.IsNull())
 280    {
 281      DoOpenError("Can not get directory name");
 282      return 1;
 283    }
 284   filename = filename +"/"+ fFileName;
 285   fRunLoader = AliRunLoader::Open(filename,AliConfig::fgkDefaultEventFolderName);
 286   if( fRunLoader == 0x0)
 287    {
 288      DoOpenError("Can not open session.");
 289      return 1;
 290    }
 291
 292   fRunLoader->LoadHeader();
 293   if (fRunLoader->GetNumberOfEvents() <= 0)
 294    {
 295      DoOpenError("There is no events in this directory.");
 296      return 1;
 297    }
 298
 299   if (fRunLoader->LoadKinematics())
 300    {
 301      DoOpenError("Error occured while loading kinematics.");
 302      return 1;
 303    }
 304   fTPCLoader = (AliTPCLoader*)fRunLoader->GetLoader("TPCLoader");
 305   if ( fTPCLoader == 0x0)
 306    {
 307      DoOpenError("Exiting due to problems with opening files.");
 308      return 1;
 309    }
 310   Info("OpenNextSession","________________________________________________________");
 311   Info("OpenNextSession","Found %d event(s) in directory %s",
 312         fRunLoader->GetNumberOfEvents(),GetDirName(fCurrentDir).Data());
 313   Float_t mf;
 314   if (fUseMagFFromRun)
 315    {
 316      fRunLoader->LoadgAlice();
 317      mf = fRunLoader->GetAliRun()->Field()->SolenoidField();
 318      Info("OpenNextSession","Setting Magnetic Field from run: B=%fT",mf/10.);
 319      fRunLoader->UnloadgAlice();
 320    }
 321   else
 322    {
 323      Info("OpenNextSession","Setting Own Magnetic Field: B=%fT",fMagneticField);
 324      if (fMagneticField == 0x0)
 325       {
 326         Fatal("OpenNextSession","Magnetic field can not be 0.");
 327         return 1;//pro forma
 328       }
 329      mf = fMagneticField*10.;
 330    }
 331   AliKalmanTrack::SetConvConst(1000/0.299792458/mf);
 332
 333
 334   fRunLoader->CdGAFile();
 335   AliTPCParam *TPCParam= (AliTPCParam*)gDirectory->Get("75x40_100x60");
 336   if (!TPCParam)
 337    {
 338     TPCParam=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
 339     if (!TPCParam)
 340      {
 341        DoOpenError("TPC parameters have not been found !\n");
 342        return 1;
 343      }
 344    }
 345
 346   if (fTPCLoader->LoadTracks())
 347    {
 348      DoOpenError("Error occured while loading TPC tracks.");
 349      return 1;
 350    }
 351
 352   fCurrentEvent = 0;
 353   return 0;
 354 }
 355 /********************************************************************/
 356
 357 void AliHBTReaderTPC::DoOpenError( const char *va_(fmt), ...)
 358 {
 359   // Does error display and clean-up in case error caught on Open Next Session
 360
 361    va_list ap;
 362    va_start(ap,va_(fmt));
 363    Error("OpenNextSession", va_(fmt), ap);
 364    va_end(ap);
 365
 366    delete fRunLoader;
 367    fRunLoader = 0x0;
 368    fTPCLoader = 0x0;
 369    fCurrentDir++;
 370 }
 371
 372 /********************************************************************/
 373 Bool_t AliHBTReaderTPC::CheckTrack(AliTPCtrack* t)
 374 {
 375   //Performs check of the track
 376   if ( (t->GetNumberOfClusters() > fNClustMax) || (t->GetNumberOfClusters() < fNClustMin) ) return kTRUE;
 377
 378   Double_t cc[15];
 379   t->GetCovariance(cc);
 380   if ( (cc[9] < fC44Min) || (cc[9] > fC44Max) ) return kTRUE;
 381
 382   Float_t chisqpercl = t->GetChi2()/((Double_t)t->GetNumberOfClusters());
 383   if ( (chisqpercl < fNChi2PerClustMin) || (chisqpercl > fNChi2PerClustMax) ) return kTRUE;
 384
 385   return kFALSE;
 386
 387 }
 388 /********************************************************************/
 389
 390 void AliHBTReaderTPC::SetNClustersRange(Int_t min,Int_t max)
 391 {
 392  //sets range of Number Of Clusters that tracks have to have
 393  fNClustMin = min;
 394  fNClustMax = max;
 395 }
 396 /********************************************************************/
 397
 398 void AliHBTReaderTPC::SetChi2PerCluserRange(Float_t min, Float_t max)
 399 {
 400   //sets range of Chi2 per Cluster
 401   fNChi2PerClustMin = min;
 402   fNChi2PerClustMax = max;
 403 }
 404 /********************************************************************/
 405
 406 void AliHBTReaderTPC::SetC44Range(Float_t min, Float_t max)
 407 {
 408  //Sets range of C44 parameter of covariance matrix of the track
 409  //it defines uncertainty of the momentum
 410   fC44Min = min;
 411   fC44Max = max;
 412 }
 413
 414 /********************************************************************/
 415 /********************************************************************/
 416 /********************************************************************/
 417
 418 /*
 419 void (AliTPCtrack* iotrack, Double_t curv)
 420 {
 421   Double_t x[5];
 422   iotrack->GetExternalPara
 423   //xk is a
 424   Double_t fP4=iotrack->GetC();
 425   Double_t fP2=iotrack->GetEta();
 426
 427   Double_t x1=fX, x2=x1+(xk-x1), dx=x2-x1, y1=fP0, z1=fP1;
 428   Double_t c1=fP4*x1 - fP2, r1=sqrt(1.- c1*c1);
 429   Double_t c2=fP4*x2 - fP2, r2=sqrt(1.- c2*c2);
 430
 431   fP0 += dx*(c1+c2)/(r1+r2);
 432   fP1 += dx*(c1+c2)/(c1*r2 + c2*r1)*fP3;
 433
 434 }
 435 */
 436