HBTAN/AliHBTReaderTPC.cxx

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