HBTAN/AliHBTReaderITSv2.cxx

   1
   2 #include "AliHBTReaderITSv2.h"
   3
   4 #include <Riostream.h>
   5 #include <Riostream.h>
   6 #include <TString.h>
   7 #include <TObjString.h>
   8 #include <TTree.h>
   9 #include <TFile.h>
  10 #include <TParticle.h>
  11
  12 #include <AliRun.h>
  13 #include <AliRunLoader.h>
  14 #include <AliLoader.h>
  15 #include <AliMagF.h>
  16 #include <AliITS.h>
  17 #include <AliITStrackV2.h>
  18 #include <AliITStrackerV2.h>
  19 #include <AliITSgeom.h>
  20
  21 #include "AliHBTRun.h"
  22 #include "AliHBTEvent.h"
  23 #include "AliHBTParticle.h"
  24 #include "AliHBTParticleCut.h"
  25
  26
  27 ClassImp(AliHBTReaderITSv2)
  28
  29 AliHBTReaderITSv2::AliHBTReaderITSv2():fFileName("galice.root")
  30 {
  31   //constructor,
  32   //Defaults:
  33   //  galicefilename = "galice.root"
  34   fParticles = 0x0;
  35   fTracks    = 0x0;
  36   fIsRead = kFALSE;
  37 }
  38 /********************************************************************/
  39
  40 AliHBTReaderITSv2::AliHBTReaderITSv2(const Char_t* galicefilename):fFileName(galicefilename)
  41 {
  42   //constructor,
  43   //Defaults:
  44   //  galicefilename = "galice.root"
  45   fParticles = new AliHBTRun();
  46   fTracks    = new AliHBTRun();
  47   fIsRead = kFALSE;
  48 }
  49 /********************************************************************/
  50
  51 AliHBTReaderITSv2::AliHBTReaderITSv2(TObjArray* dirs, const Char_t* galicefilename):
  52        AliHBTReader(dirs), fFileName(galicefilename)
  53 {
  54   //constructor,
  55   //Defaults:
  56   //  galicefilename = "galice.root"
  57
  58   fParticles = new AliHBTRun();
  59   fTracks    = new AliHBTRun();
  60   fIsRead = kFALSE;
  61 }
  62 /********************************************************************/
  63
  64 AliHBTReaderITSv2::~AliHBTReaderITSv2()
  65  {
  66    if (fParticles) delete fParticles;
  67    if (fTracks) delete fTracks;
  68  }
  69 /********************************************************************/
  70 /********************************************************************/
  71
  72 AliHBTEvent* AliHBTReaderITSv2::GetParticleEvent(Int_t n)
  73  {
  74  //returns Nth event with simulated particles
  75  if (!fIsRead)
  76   {
  77     if (fParticles == 0x0) fParticles = new AliHBTRun();
  78     if (fTracks == 0x0) fTracks    = new AliHBTRun();
  79     if(Read(fParticles,fTracks))
  80      {
  81        Error("GetParticleEvent","Error in reading");
  82        return 0x0;
  83      }
  84   }
  85  return (fParticles)?fParticles->GetEvent(n):0x0;
  86 }
  87 /********************************************************************/
  88
  89 AliHBTEvent* AliHBTReaderITSv2::GetTrackEvent(Int_t n)
  90  {
  91  //returns Nth event with reconstructed tracks
  92  if (!fIsRead)
  93   {
  94     if (fParticles == 0x0) fParticles = new AliHBTRun();
  95     if (fTracks == 0x0) fTracks    = new AliHBTRun();
  96     if(Read(fParticles,fTracks))
  97      {
  98        Error("GetTrackEvent","Error in reading");
  99        return 0x0;
 100      }
 101    }
 102   return (fTracks)?fTracks->GetEvent(n):0x0;
 103  }
 104 /********************************************************************/
 105
 106 Int_t AliHBTReaderITSv2::GetNumberOfPartEvents()
 107  {
 108  //returns number of events of particles
 109   if (!fIsRead)
 110    {
 111     if (fParticles == 0x0) fParticles = new AliHBTRun();
 112     if (fTracks == 0x0) fTracks    = new AliHBTRun();
 113     if(Read(fParticles,fTracks))
 114      {
 115        Error("GetNumberOfPartEvents","Error in reading");
 116        return 0;
 117      }
 118    }
 119   return (fParticles)?fParticles->GetNumberOfEvents():0;
 120  }
 121
 122 /********************************************************************/
 123 Int_t AliHBTReaderITSv2::GetNumberOfTrackEvents()
 124  {
 125  //returns number of events of tracks
 126   if (!fIsRead)
 127    {
 128     if (fParticles == 0x0) fParticles = new AliHBTRun();
 129     if (fTracks == 0x0) fTracks    = new AliHBTRun();
 130     if(Read(fParticles,fTracks))
 131      {
 132        Error("GetNumberOfTrackEvents","Error in reading");
 133        return 0;
 134      }
 135    }
 136   return (fTracks)?fTracks->GetNumberOfEvents():0;
 137  }
 138 /********************************************************************/
 139 /********************************************************************/
 140
 141
 142 Int_t AliHBTReaderITSv2::Read(AliHBTRun* particles, AliHBTRun *tracks)
 143 {
 144 //reads data
 145  Int_t Nevents = 0; //number of events found in given directory
 146  Int_t Ndirs; //number of the directories to be read
 147  Int_t Ntracks; //number of tracks in current event
 148  Int_t currentdir = 0; //number of events in the current directory
 149  Int_t totalNevents = 0; //total number of events read from all directories up to now
 150  register Int_t i = 0; //iterator
 151
 152 // AliITStrackerV2 *tracker; // ITS tracker - used for cooking labels
 153  TTree *tracktree; // tree for tracks
 154
 155  Double_t xk;
 156  Double_t par[5]; //Kalman track parameters
 157  Float_t phi, lam, pt;//angles and transverse momentum
 158  Int_t label; //label of the current track
 159
 160  AliITStrackV2 *iotrack = 0x0; //buffer track for reading data from tree
 161
 162  if (!particles) //check if an object is instatiated
 163   {
 164     Error("Read"," particles object must instatiated before passing it to the reader");
 165   }
 166  if (!tracks)  //check if an object is instatiated
 167   {
 168     Error("Read"," tracks object must instatiated before passing it to the reader");
 169   }
 170  particles->Reset();//clear runs == delete all old events
 171  tracks->Reset();
 172
 173  if (fDirs) //if array with directories is supplied by user
 174   {
 175     Ndirs = fDirs->GetEntries(); //get the number if directories
 176   }
 177  else
 178   {
 179     Ndirs = 0; //if the array is not supplied read only from current directory
 180   }
 181
 182 // cout<<"Found "<<Ndirs<<" directory entries"<<endl;
 183
 184  do //do while is good even if Ndirs==0 (than read from current directory)
 185    {
 186     TString filename = GetDirName(currentdir);
 187     if (filename.IsNull())
 188      {
 189        Error("Read","Can not get directory name");
 190        currentdir++;
 191        continue;
 192      }
 193     filename = filename +"/"+ fFileName;
 194     AliRunLoader* rl = AliRunLoader::Open(filename);
 195     if( rl == 0x0)
 196      {
 197        Error("Read","Exiting due to problems with opening files.");
 198        currentdir++;
 199        continue;
 200      }
 201
 202     rl->LoadHeader();
 203     rl->LoadKinematics();
 204     rl->LoadgAlice();
 205     gAlice = rl->GetAliRun();
 206     AliITS* its = (AliITS*)gAlice->GetModule("ITS");
 207
 208     AliLoader* itsl = rl->GetLoader("ITSLoader");
 209
 210     if ((its == 0x0) || ( itsl== 0x0))
 211      {
 212        Error("Read","Can not found ITS in this run");
 213        delete rl;
 214        rl = 0x0;
 215        currentdir++;
 216        continue;
 217      }
 218     Nevents = rl->GetNumberOfEvents();
 219
 220     if (Nevents > 0)//check if tree E exists
 221      {
 222       Info("Read","________________________________________________________");
 223       Info("Read","Found %d event(s) in directory %s",Nevents,GetDirName(currentdir).Data());
 224       Float_t mf;
 225       if (fUseMagFFromRun)
 226        {
 227          mf = gAlice->Field()->SolenoidField();
 228          Info("Read","Setting Magnetic Field from run: B=%fT",mf/10.);
 229        }
 230       else
 231        {
 232          Info("Read","Setting Own Magnetic Field: B=%fT",fMagneticField);
 233          mf = fMagneticField*10.;
 234        }
 235       AliKalmanTrack::SetConvConst(1000/0.299792458/mf);
 236       if (iotrack == 0x0) iotrack = new AliITStrackV2();
 237      }
 238     else
 239      {//if not return an error
 240        Error("Read","No events in this run");
 241        delete rl;
 242        rl = 0x0;
 243        currentdir++;
 244        continue;
 245      }
 246
 247     AliITSgeom *geom= its->GetITSgeom();
 248     if (!geom)
 249      {
 250        Error("Read","Can't get the ITS geometry!");
 251        delete rl;
 252        rl = 0x0;
 253        currentdir++;
 254        continue;
 255      }
 256
 257     itsl->LoadTracks();
 258
 259     for(Int_t currentEvent =0; currentEvent<Nevents;currentEvent++)//loop over all events
 260      {
 261        cout<<"Reading Event "<<currentEvent<<endl;
 262        rl->GetEvent(currentEvent);
 263        tracktree=itsl->TreeT();
 264
 265        if (!tracktree)
 266          {
 267            Error("Read","Can't get a tree with ITS tracks");
 268            continue;
 269          }
 270        TBranch *tbranch=tracktree->GetBranch("tracks");
 271        Ntracks=(Int_t)tracktree->GetEntries();
 272
 273        Int_t accepted = 0;
 274        Int_t tpcfault = 0;
 275        Int_t itsfault = 0;
 276        for (i=0; i<Ntracks; i++) //loop over all tpc tracks
 277         {
 278           if(i%100 == 0)cout<<"all: "<<i<<"   accepted: "<<accepted<<"   tpc faults: "<<tpcfault<<"\r";
 279
 280           tbranch->SetAddress(&iotrack);
 281           tracktree->GetEvent(i);
 282
 283           label=iotrack->GetLabel();
 284           if (label < 0)
 285            {
 286              tpcfault++;
 287              continue;
 288            }
 289
 290           TParticle *p = (TParticle*)gAlice->Particle(label);
 291           if(p == 0x0) continue; //if returned pointer is NULL
 292           if(p->GetPDG() == 0x0) continue; //if particle has crezy PDG code (not known to our database)
 293
 294           if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type
 295                                               //if not take next partilce
 296
 297           AliHBTParticle* part = new AliHBTParticle(*p,i);
 298           if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
 299                                                   //if it does not delete it and take next good track
 300
 301           iotrack->PropagateTo(3.,0.0028,65.19);
 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(), i, 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           accepted++;
 329         }//end of loop over tracks in the event
 330
 331        totalNevents++;
 332        cout<<"all: "<<i<<"   accepted: "<<accepted<<"   tpc faults: "<<tpcfault<<"   its faults: "<<itsfault<<endl;
 333
 334      }//end of loop over events in current directory
 335     delete rl;
 336     currentdir++;
 337    }while(currentdir < Ndirs);//end of loop over directories specified in fDirs Obj Array
 338
 339  delete iotrack;
 340  fIsRead = kTRUE;
 341  return 0;
 342 }
 343
 344 /********************************************************************/
 345 /********************************************************************/
 346
 347