HBTAN/AliHBTReaderTPC.cxx

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