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