HBTAN/AliHBTReaderPPprod.cxx

   1 #include "AliHBTReaderPPprod.h"
   2
   3 #include <iostream.h>
   4 #include <fstream.h>
   5 #include <TString.h>
   6 #include <TTree.h>
   7 #include <TFile.h>
   8
   9 #include <AliTPCtrack.h>
  10 #include <AliTPCParam.h>
  11 #include <AliTPCtracker.h>
  12
  13 #include "AliRun.h"
  14 #include "AliHBTRun.h"
  15 #include "AliHBTEvent.h"
  16 #include "AliHBTParticle.h"
  17 #include "AliHBTParticleCut.h"
  18
  19
  20 ClassImp(AliHBTReaderPPprod)
  21
  22 AliHBTReaderPPprod::
  23  AliHBTReaderPPprod(const Char_t* trackfilename,const Char_t* clusterfilename,
  24                  const Char_t* goodtracksfilename,const Char_t* galicefilename):
  25                  fTrackFileName(trackfilename),fClusterFileName(clusterfilename),
  26                  fGAliceFileName(galicefilename),
  27                  fGoodTPCTracksFileName(goodtracksfilename)
  28 {
  29   //constructor, only file names are set
  30   //Defaults:
  31   //  trackfilename = "AliTPCtracks.root"
  32   //  clusterfilename = "AliTPCclusters.root"
  33   //  goodtracksfilename = "good_tracks_tpc"
  34   //  galicefilename = ""  - this means: Do not open gAlice file -
  35   //                         just leave the global pointer untached
  36
  37   fParticles = new AliHBTRun();
  38   fTracks    = new AliHBTRun();
  39
  40   fTracksFile   = 0x0;  //files are opened during reading only
  41   fClustersFile = 0x0;
  42
  43   fIsRead = kFALSE;
  44 }
  45 /********************************************************************/
  46
  47 AliHBTReaderPPprod::~AliHBTReaderPPprod()
  48  {
  49    delete fParticles;
  50    delete fTracks;
  51  }
  52 /********************************************************************/
  53
  54 AliHBTEvent* AliHBTReaderPPprod::GetParticleEvent(Int_t n)
  55  {
  56  //returns Nth event with simulated particles
  57    if (!fIsRead) Read(fParticles,fTracks);
  58    return fParticles->GetEvent(n);
  59  }
  60 /********************************************************************/
  61 AliHBTEvent* AliHBTReaderPPprod::GetTrackEvent(Int_t n)
  62  {
  63  //returns Nth event with reconstructed tracks
  64    if (!fIsRead) Read(fParticles,fTracks);
  65    return fTracks->GetEvent(n);
  66  }
  67 /********************************************************************/
  68
  69 Int_t AliHBTReaderPPprod::GetNumberOfPartEvents()
  70  {
  71  //returns number of events of particles
  72    if (!fIsRead) Read(fParticles,fTracks);
  73    return fParticles->GetNumberOfEvents();
  74  }
  75
  76 /********************************************************************/
  77 Int_t AliHBTReaderPPprod::GetNumberOfTrackEvents()
  78  {
  79  //returns number of events of tracks
  80   if (!fIsRead) Read(fParticles,fTracks);
  81   return fTracks->GetNumberOfEvents();
  82  }
  83 /********************************************************************/
  84
  85
  86 Int_t AliHBTReaderPPprod::Read(AliHBTRun* particles, AliHBTRun *tracks)
  87  {
  88  //reads data and puts put to the particles and tracks objects
  89  //reurns 0 if everything is OK
  90  //
  91   Int_t i; //iterator and some temprary values
  92   Int_t Nevents;
  93   if (!particles) //check if an object is instatiated
  94    {
  95      Error("AliHBTReaderPPprod::Read"," particles object must instatiated before passing it to the reader");
  96    }
  97   if (!tracks)  //check if an object is instatiated
  98    {
  99      Error("AliHBTReaderPPprod::Read"," tracks object must instatiated before passing it to the reader");
 100    }
 101   particles->Reset();//clear runs == delete all old events
 102   tracks->Reset();
 103
 104   if( (i=OpenFiles()) )
 105    {
 106      Error("AliHBTReaderPPprod::Read","Exiting due to problems with opening files. Errorcode %d",i);
 107      return i;
 108    }
 109
 110   AliGoodTracksPP *goodTPCTracks = new AliGoodTracksPP(fGoodTPCTracksFileName);
 111   if (!goodTPCTracks)
 112    {
 113      Error("AliHBTReaderPPprod::Read","Exiting due to problems with opening files. Errorcode %d",i);
 114      return 1;
 115    }
 116
 117   Nevents = 100;
 118
 119   fClustersFile->cd();
 120   AliTPCParam *TPCParam= (AliTPCParam*)fClustersFile->Get("75x40_100x60");
 121   if (!TPCParam)
 122     {
 123      Error("AliHBTReaderPPprod::Read","TPC parameters have not been found !\n");
 124      return 1;
 125     }
 126
 127   TObjArray *tarray = new TObjArray(5000);
 128   tarray->SetOwner();//  this causes memory leak, but in some cases deleting is infinite loop
 129
 130
 131   for(Int_t currentEvent =0; currentEvent<Nevents;currentEvent++)
 132    {
 133      cout<<"Reading Event "<<currentEvent<<endl;
 134      /**************************************/
 135       /**************************************/
 136        /**************************************/
 137          fTracksFile->cd();
 138          Char_t  treename[100];
 139          sprintf(treename,"TreeT_TPC_%d",currentEvent);
 140
 141          TTree *tracktree=0;
 142
 143          tracktree=(TTree*)fTracksFile->Get(treename);
 144          if (!tracktree)
 145           {
 146             Error("AliHBTReaderPPprod::Read","Can't get a tree with TPC tracks !\n");
 147             return 1;
 148           }
 149
 150          TBranch *trackbranch=tracktree->GetBranch("tracks");
 151          if (!trackbranch)
 152           {
 153             Error("AliHBTReaderPPprod::Read","Can't get a branch with TPC tracks !\n");
 154             return 2;
 155           }
 156          Int_t NTPCtracks=(Int_t)tracktree->GetEntries();
 157          cout<<"Found "<<NTPCtracks<<" TPC tracks.\n";
 158          //Copy tracks to array
 159
 160          AliTPCtrack *iotrack=0;
 161
 162          fClustersFile->cd();
 163          AliTPCtracker *tracker = new AliTPCtracker(TPCParam,currentEvent);
 164          if (!tracker)
 165           {
 166             Error("AliHBTReaderPPprod::Read","Can't get a tracker !\n");
 167             return 3;
 168           }
 169          tracker->LoadInnerSectors();
 170          tracker->LoadOuterSectors();
 171
 172          for (i=0; i<NTPCtracks; i++)
 173           {
 174             iotrack=new AliTPCtrack;
 175             trackbranch->SetAddress(&iotrack);
 176             tracktree->GetEvent(i);
 177             tracker->CookLabel(iotrack,0.1);
 178             tarray->AddLast(iotrack);
 179           }
 180
 181
 182          fTracksFile->Delete(treename);//delete tree from memmory (and leave untached on disk)- we do not need it any more
 183          fTracksFile->Delete("tracks");
 184
 185          delete tracker;
 186
 187          tracker = 0x0;
 188          trackbranch = 0x0;
 189          tracktree = 0x0;
 190
 191          Int_t & ngood = goodTPCTracks->fGoodInEvent[currentEvent];
 192
 193          Double_t xk;
 194          Double_t par[5];
 195          Float_t phi, lam, pt;
 196          Int_t label;
 197          Bool_t found;
 198
 199          for (i=0; i<ngood; i++)
 200           {
 201             const struct GoodTrack & gt = goodTPCTracks->GetTrack(currentEvent,i);
 202
 203             if(Pass(gt.code)) continue;
 204
 205             label = gt.lab;
 206             found = kFALSE; //guard in case we don't find track with such a label
 207             for (Int_t j=0;j<NTPCtracks;j++)
 208               {
 209                 iotrack = (AliTPCtrack*)tarray->At(j);
 210                 if (iotrack->GetLabel() == label)
 211                   {
 212                     found = kTRUE;
 213                     break;
 214                   }
 215               }
 216             if(!found)
 217               {
 218                 Warning("Read",
 219                 "Sth is going wrong with tracks - there is no TPC track corresponding to goodtrack.\nGood tack label %d",label);
 220                 continue; //put comunicate on the screen and continue loop
 221               }
 222
 223             Double_t mass = TDatabasePDG::Instance()->GetParticle(gt.code)->Mass();
 224             Double_t pEtot = TMath::Sqrt(gt.px*gt.px + gt.py*gt.py + gt.pz*gt.pz + mass*mass);
 225
 226             AliHBTParticle* part = new AliHBTParticle(gt.code, gt.px, gt.py, gt.pz, pEtot, gt.x, gt.y, gt.z, 0.0);
 227             if(Pass(part)) continue;
 228
 229
 230             iotrack->PropagateTo(gt.x);
 231             iotrack->GetExternalParameters(xk,par);
 232             phi=TMath::ASin(par[2]) + iotrack->GetAlpha();
 233             if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
 234             if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
 235             lam=par[3];
 236             pt=1.0/TMath::Abs(par[4]);
 237
 238             Double_t tpx = pt * TMath::Cos(phi); //track x coordinate of momentum
 239             Double_t tpy = pt * TMath::Sin(phi); //track x coordinate of momentum
 240             Double_t tpz = pt * lam;
 241
 242             Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);
 243
 244             AliHBTParticle* track = new AliHBTParticle(gt.code, tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
 245             if(Pass(track)) continue;
 246
 247             particles->AddParticle(currentEvent,part);
 248             tracks->AddParticle(currentEvent,track);
 249
 250           }
 251          tarray->Clear();
 252        /**************************************/
 253       /**************************************/
 254      /**************************************/
 255    }
 256
 257   CloseFiles();
 258   delete tarray;
 259   delete goodTPCTracks;
 260   fIsRead = kTRUE;
 261   return 0;
 262  }
 263
 264 /********************************************************************/
 265 Int_t AliHBTReaderPPprod::OpenFiles()
 266 {
 267
 268    fTracksFile = 0;
 269    fTracksFile=TFile::Open("AliTPCtracks.root");
 270    if (!fTracksFile->IsOpen())
 271      {
 272        Error("AliHBTReaderPPprod::OpenFiles","Can't open AliTPCtracks.root");
 273        return 1;
 274      }
 275
 276    fClustersFile = 0;
 277
 278    fClustersFile=TFile::Open("AliTPCclusters.root");
 279    if (!fClustersFile->IsOpen())
 280     {
 281       Error("AliHBTReaderPPprod::OpenFiles","Can't open AliTPCclusters.root");
 282       return 2;
 283     }
 284
 285
 286
 287  return 0;
 288 }
 289
 290
 291
 292 /********************************************************************/
 293
 294 void AliHBTReaderPPprod::CloseFiles()
 295 {
 296   fTracksFile->Close();
 297   fClustersFile->Close();
 298 }
 299
 300 /********************************************************************/
 301
 302 /********************************************************************/
 303 /********************************************************************/
 304 /********************************************************************/
 305
 306
 307 AliGoodTracksPP::~AliGoodTracksPP()
 308 {
 309  delete [] fGoodInEvent;
 310  for (Int_t i = 0;i<fNevents;i++)
 311    delete [] fData[i];
 312  delete [] fData;
 313 }
 314 /********************************************************************/
 315 AliGoodTracksPP::AliGoodTracksPP(const TString& infilename)
 316 {
 317
 318   fNevents = 100;
 319   ifstream in(infilename.Data());
 320
 321   if(!in)
 322     {
 323       cerr<<"Can not open file with Good TPC Tracks named:"<<infilename.Data()<<endl;
 324       delete this;
 325       return;
 326     }
 327
 328
 329   fGoodInEvent = new Int_t[fNevents];
 330   fData = new struct GoodTrack* [fNevents];
 331
 332   Int_t i;
 333   for( i = 0;i<fNevents;i++)
 334    {
 335     fGoodInEvent[i] =0;
 336     fData[i] = new struct GoodTrack[500];
 337    }
 338   Float_t tmp;
 339   Int_t evno;
 340   while(in>>evno)
 341    {
 342     if(fGoodInEvent[evno]>=500)
 343      {
 344       cerr<<"AliGoodTracksPP::AliGoodTracksPP() : Not enough place in the array\n";
 345       continue;
 346      }
 347     in>>fData[evno][fGoodInEvent[evno]].lab;
 348     in>>fData[evno][fGoodInEvent[evno]].code;
 349     in>>fData[evno][fGoodInEvent[evno]].px;
 350     in>>fData[evno][fGoodInEvent[evno]].py;
 351     in>>fData[evno][fGoodInEvent[evno]].pz;
 352     in>>fData[evno][fGoodInEvent[evno]].x;
 353     in>>fData[evno][fGoodInEvent[evno]].y;
 354     in>>fData[evno][fGoodInEvent[evno]].z;
 355
 356     in>>tmp;
 357     in>>tmp;
 358     in>>tmp;
 359
 360  /* cout<<evno<<" ";
 361   cout<<fData[evno][fGoodInEvent[evno]].lab;
 362   cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].code;
 363   cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].px;
 364   cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].py;
 365   cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].pz;
 366   cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].x;
 367   cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].y;
 368   cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].z;
 369   cout<<"\n";
 370  */
 371   fGoodInEvent[evno]++;
 372  }
 373  in.close();
 374  cout<<"AliGoodTracksPP::AliGoodTracksPP()  ....  Done\n";
 375 }
 376
 377
 378
 379 const GoodTrack& AliGoodTracksPP::GetTrack(Int_t event, Int_t n) const
 380  {
 381
 382   if( (event>fNevents) || (event<0))
 383    {
 384      gAlice->Fatal("AliGoodTracksPP::GetTrack","No such Event %d",event);
 385    }
 386   if( (n>fGoodInEvent[event]) || (n<0))
 387    {
 388      gAlice->Fatal("AliGoodTracksPP::GetTrack","No such Good TPC Track %d",n);
 389    }
 390   return fData[event][n];
 391
 392  }