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