HBTAN/AliHBTReaderTPC.cxx

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