HBTAN/AliHBTReaderPPprod.cxx

   1 #include "AliHBTReaderPPprod.h"
   2
   3 #include <Riostream.h>
   4 #include <Riostream.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   cout<<"New I/O not implemented\n";
  92   return 0;
  93
  94   Int_t i; //iterator and some temprary values
  95   Int_t Nevents;
  96   if (!particles) //check if an object is instatiated
  97    {
  98      Error("AliHBTReaderPPprod::Read"," particles object must instatiated before passing it to the reader");
  99    }
 100   if (!tracks)  //check if an object is instatiated
 101    {
 102      Error("AliHBTReaderPPprod::Read"," tracks object must instatiated before passing it to the reader");
 103    }
 104   particles->Reset();//clear runs == delete all old events
 105   tracks->Reset();
 106
 107   if( (i=OpenFiles()) )
 108    {
 109      Error("AliHBTReaderPPprod::Read","Exiting due to problems with opening files. Errorcode %d",i);
 110      return i;
 111    }
 112
 113   AliGoodTracksPP *goodTPCTracks = new AliGoodTracksPP(fGoodTPCTracksFileName);
 114   if (!goodTPCTracks)
 115    {
 116      Error("AliHBTReaderPPprod::Read","Exiting due to problems with opening files. Errorcode %d",i);
 117      return 1;
 118    }
 119
 120   Nevents = 100;
 121
 122   fClustersFile->cd();
 123   AliTPCParam *TPCParam= (AliTPCParam*)fClustersFile->Get("75x40_100x60");
 124   if (!TPCParam)
 125     {
 126      Error("AliHBTReaderPPprod::Read","TPC parameters have not been found !\n");
 127      return 1;
 128     }
 129
 130   TObjArray *tarray = new TObjArray(5000);
 131   tarray->SetOwner();//  this causes memory leak, but in some cases deleting is infinite loop
 132
 133
 134   for(Int_t currentEvent =0; currentEvent<Nevents;currentEvent++)
 135    {
 136      cout<<"Reading Event "<<currentEvent<<endl;
 137      /**************************************/
 138       /**************************************/
 139        /**************************************/
 140          fTracksFile->cd();
 141          Char_t  treename[100];
 142          sprintf(treename,"TreeT_TPC_%d",currentEvent);
 143
 144          TTree *tracktree=0;
 145
 146          tracktree=(TTree*)fTracksFile->Get(treename);
 147          if (!tracktree)
 148           {
 149             Error("AliHBTReaderPPprod::Read","Can't get a tree with TPC tracks !\n");
 150             return 1;
 151           }
 152
 153          TBranch *trackbranch=tracktree->GetBranch("tracks");
 154          if (!trackbranch)
 155           {
 156             Error("AliHBTReaderPPprod::Read","Can't get a branch with TPC tracks !\n");
 157             return 2;
 158           }
 159          Int_t NTPCtracks=(Int_t)tracktree->GetEntries();
 160          cout<<"Found "<<NTPCtracks<<" TPC tracks.\n";
 161          //Copy tracks to array
 162
 163          AliTPCtrack *iotrack=0;
 164
 165          fClustersFile->cd();
 166 printf("This method is not converted to the NewIO !\n"); //I.B.
 167 return 1; //I.B.
 168          AliTPCtracker *tracker = new AliTPCtracker(TPCParam); //I.B.
 169          //AliTPCtracker *tracker = new AliTPCtracker(TPCParam,currentEvent,"");
 170          if (!tracker)
 171           {
 172             Error("AliHBTReaderPPprod::Read","Can't get a tracker !\n");
 173             return 3;
 174           }
 175 tracker->LoadClusters(0);//I.Belikov, "0" should be a pointer to a tree
 176
 177          for (i=0; i<NTPCtracks; i++)
 178           {
 179             iotrack=new AliTPCtrack;
 180             trackbranch->SetAddress(&iotrack);
 181             tracktree->GetEvent(i);
 182             tracker->CookLabel(iotrack,0.1);
 183             tarray->AddLast(iotrack);
 184           }
 185
 186
 187          fTracksFile->Delete(treename);//delete tree from memmory (and leave untached on disk)- we do not need it any more
 188          fTracksFile->Delete("tracks");
 189
 190          delete tracker;
 191
 192          tracker = 0x0;
 193          trackbranch = 0x0;
 194          tracktree = 0x0;
 195
 196          Int_t & ngood = goodTPCTracks->fGoodInEvent[currentEvent];
 197
 198          Double_t xk;
 199          Double_t par[5];
 200          Float_t phi, lam, pt;
 201          Int_t label;
 202          Bool_t found;
 203
 204          for (i=0; i<ngood; i++)
 205           {
 206             const struct GoodTrack & gt = goodTPCTracks->GetTrack(currentEvent,i);
 207
 208             if(Pass(gt.code)) continue;
 209
 210             label = gt.lab;
 211             found = kFALSE; //guard in case we don't find track with such a label
 212             for (Int_t j=0;j<NTPCtracks;j++)
 213               {
 214                 iotrack = (AliTPCtrack*)tarray->At(j);
 215                 if (iotrack->GetLabel() == label)
 216                   {
 217                     found = kTRUE;
 218                     break;
 219                   }
 220               }
 221             if(!found)
 222               {
 223                 Warning("Read",
 224                 "Sth is going wrong with tracks - there is no TPC track corresponding to goodtrack.\nGood tack label %d",label);
 225                 continue; //put comunicate on the screen and continue loop
 226               }
 227
 228             Double_t mass = TDatabasePDG::Instance()->GetParticle(gt.code)->Mass();
 229             Double_t pEtot = TMath::Sqrt(gt.px*gt.px + gt.py*gt.py + gt.pz*gt.pz + mass*mass);
 230
 231             AliHBTParticle* part = new AliHBTParticle(gt.code,i, gt.px, gt.py, gt.pz, pEtot, gt.x, gt.y, gt.z, 0.0);
 232             if(Pass(part)) continue;
 233
 234
 235             iotrack->PropagateTo(gt.x);
 236             iotrack->GetExternalParameters(xk,par);
 237             phi=TMath::ASin(par[2]) + iotrack->GetAlpha();
 238             if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
 239             if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
 240             lam=par[3];
 241             pt=1.0/TMath::Abs(par[4]);
 242
 243             Double_t tpx = pt * TMath::Cos(phi); //track x coordinate of momentum
 244             Double_t tpy = pt * TMath::Sin(phi); //track x coordinate of momentum
 245             Double_t tpz = pt * lam;
 246
 247             Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);
 248
 249             AliHBTParticle* track = new AliHBTParticle(gt.code,i, tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
 250             if(Pass(track)) continue;
 251
 252             particles->AddParticle(currentEvent,part);
 253             tracks->AddParticle(currentEvent,track);
 254
 255           }
 256          tarray->Clear();
 257        /**************************************/
 258       /**************************************/
 259      /**************************************/
 260    }
 261
 262   CloseFiles();
 263   delete tarray;
 264   delete goodTPCTracks;
 265   fIsRead = kTRUE;
 266   return 0;
 267  }
 268
 269 /********************************************************************/
 270 Int_t AliHBTReaderPPprod::OpenFiles()
 271 {
 272
 273    fTracksFile = 0;
 274    fTracksFile=TFile::Open("AliTPCtracks.root");
 275    if (!fTracksFile->IsOpen())
 276      {
 277        Error("AliHBTReaderPPprod::OpenFiles","Can't open AliTPCtracks.root");
 278        return 1;
 279      }
 280
 281    fClustersFile = 0;
 282
 283    fClustersFile=TFile::Open("AliTPCclusters.root");
 284    if (!fClustersFile->IsOpen())
 285     {
 286       Error("AliHBTReaderPPprod::OpenFiles","Can't open AliTPCclusters.root");
 287       return 2;
 288     }
 289
 290
 291
 292  return 0;
 293 }
 294
 295
 296
 297 /********************************************************************/
 298
 299 void AliHBTReaderPPprod::CloseFiles()
 300 {
 301   fTracksFile->Close();
 302   fClustersFile->Close();
 303 }
 304
 305 /********************************************************************/
 306
 307 /********************************************************************/
 308 /********************************************************************/
 309 /********************************************************************/
 310
 311
 312 AliGoodTracksPP::~AliGoodTracksPP()
 313 {
 314  delete [] fGoodInEvent;
 315  for (Int_t i = 0;i<fNevents;i++)
 316    delete [] fData[i];
 317  delete [] fData;
 318 }
 319 /********************************************************************/
 320 AliGoodTracksPP::AliGoodTracksPP(const TString& infilename)
 321 {
 322
 323   fNevents = 100;
 324   ifstream in(infilename.Data());
 325
 326   if(!in)
 327     {
 328       cerr<<"Can not open file with Good TPC Tracks named:"<<infilename.Data()<<endl;
 329       delete this;
 330       return;
 331     }
 332
 333
 334   fGoodInEvent = new Int_t[fNevents];
 335   fData = new struct GoodTrack* [fNevents];
 336
 337   Int_t i;
 338   for( i = 0;i<fNevents;i++)
 339    {
 340     fGoodInEvent[i] =0;
 341     fData[i] = new struct GoodTrack[500];
 342    }
 343   Float_t tmp;
 344   Int_t evno;
 345   while(in>>evno)
 346    {
 347     if(fGoodInEvent[evno]>=500)
 348      {
 349       cerr<<"AliGoodTracksPP::AliGoodTracksPP() : Not enough place in the array\n";
 350       continue;
 351      }
 352     in>>fData[evno][fGoodInEvent[evno]].lab;
 353     in>>fData[evno][fGoodInEvent[evno]].code;
 354     in>>fData[evno][fGoodInEvent[evno]].px;
 355     in>>fData[evno][fGoodInEvent[evno]].py;
 356     in>>fData[evno][fGoodInEvent[evno]].pz;
 357     in>>fData[evno][fGoodInEvent[evno]].x;
 358     in>>fData[evno][fGoodInEvent[evno]].y;
 359     in>>fData[evno][fGoodInEvent[evno]].z;
 360
 361     in>>tmp;
 362     in>>tmp;
 363     in>>tmp;
 364
 365  /* cout<<evno<<" ";
 366   cout<<fData[evno][fGoodInEvent[evno]].lab;
 367   cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].code;
 368   cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].px;
 369   cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].py;
 370   cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].pz;
 371   cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].x;
 372   cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].y;
 373   cout<<" ";cout<<fData[evno][fGoodInEvent[evno]].z;
 374   cout<<"\n";
 375  */
 376   fGoodInEvent[evno]++;
 377  }
 378  in.close();
 379  cout<<"AliGoodTracksPP::AliGoodTracksPP()  ....  Done\n";
 380 }
 381
 382
 383
 384 const GoodTrack& AliGoodTracksPP::GetTrack(Int_t event, Int_t n) const
 385  {
 386
 387   if( (event>fNevents) || (event<0))
 388    {
 389      gAlice->Fatal("AliGoodTracksPP::GetTrack","No such Event %d",event);
 390    }
 391   if( (n>fGoodInEvent[event]) || (n<0))
 392    {
 393      gAlice->Fatal("AliGoodTracksPP::GetTrack","No such Good TPC Track %d",n);
 394    }
 395   return fData[event][n];
 396
 397  }