HBTAN/AliHBTReaderTPC.cxx

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