HBTAN/AliHBTReaderTPC.cxx

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