HBTAN/AliHBTReaderITSv1.cxx

   1
   2 #include "AliHBTReaderITSv1.h"
   3 #include "AliHBTEvent.h"
   4 #include "AliHBTRun.h"
   5 #include "AliHBTParticle.h"
   6 #include "AliHBTParticleCut.h"
   7
   8 #include <iostream.h>
   9
  10 #include <TROOT.h>
  11 #include <TFile.h>
  12 #include <TTree.h>
  13 #include <TBranch.h>
  14 #include <TObjArray.h>
  15 #include <TParticle.h>
  16 #include <TString.h>
  17 #include <TObjString.h>
  18
  19 #include <AliRun.h>
  20 #include <AliMagF.h>
  21 #include <AliKalmanTrack.h>
  22 #include <AliITSIOTrack.h>
  23
  24 ClassImp(AliHBTReaderITSv1)
  25 /********************************************************************/
  26
  27 AliHBTReaderITSv1::
  28 AliHBTReaderITSv1(const Char_t* tracksfilename,const Char_t* galicefilename):
  29                  fITSTracksFileName(tracksfilename),fGAliceFileName(galicefilename)
  30  {
  31      fParticles = new AliHBTRun();
  32      fTracks    = new AliHBTRun();
  33      fIsRead = kFALSE;
  34  }
  35 /********************************************************************/
  36
  37 AliHBTReaderITSv1::
  38 AliHBTReaderITSv1(TObjArray* dirs, const Char_t* tracksfilename,const Char_t* galicefilename):
  39                  AliHBTReader(dirs),
  40                  fITSTracksFileName(tracksfilename),fGAliceFileName(galicefilename)
  41  {
  42    fParticles = new AliHBTRun();
  43    fTracks    = new AliHBTRun();
  44    fIsRead    = kFALSE;
  45  }
  46 /********************************************************************/
  47
  48 AliHBTReaderITSv1::~AliHBTReaderITSv1()
  49 {
  50    delete fParticles;
  51    delete fTracks;
  52 }
  53 /********************************************************************/
  54
  55 AliHBTEvent* AliHBTReaderITSv1::GetParticleEvent(Int_t n)
  56  {
  57  //returns Nth event with simulated particles
  58    if (!fIsRead)
  59     if(Read(fParticles,fTracks))
  60      {
  61        Error("GetParticleEvent","Error in reading");
  62        return 0x0;
  63      }
  64
  65    return fParticles->GetEvent(n);
  66  }
  67 /********************************************************************/
  68
  69 AliHBTEvent* AliHBTReaderITSv1::GetTrackEvent(Int_t n)
  70  {
  71  //returns Nth event with reconstructed tracks
  72    if (!fIsRead)
  73     if(Read(fParticles,fTracks))
  74      {
  75        Error("GetTrackEvent","Error in reading");
  76        return 0x0;
  77      }
  78    return fTracks->GetEvent(n);
  79  }
  80 /********************************************************************/
  81
  82 Int_t AliHBTReaderITSv1::GetNumberOfPartEvents()
  83  {
  84  //returns number of events of particles
  85    if (!fIsRead)
  86     if(Read(fParticles,fTracks))
  87      {
  88        Error("GetNumberOfPartEvents","Error in reading");
  89        return 0;
  90      }
  91    return fParticles->GetNumberOfEvents();
  92  }
  93
  94 /********************************************************************/
  95 Int_t AliHBTReaderITSv1::GetNumberOfTrackEvents()
  96  {
  97  //returns number of events of tracks
  98   if (!fIsRead)
  99     if(Read(fParticles,fTracks))
 100      {
 101        Error("GetNumberOfTrackEvents","Error in reading");
 102        return 0;
 103      }
 104   return fTracks->GetNumberOfEvents();
 105  }
 106 /********************************************************************/
 107
 108 Int_t AliHBTReaderITSv1::Read(AliHBTRun* particles, AliHBTRun *tracks)
 109 {
 110  cout<<"AliHBTReaderITSv1::Read()"<<endl;
 111  Int_t Nevents = 0;
 112  AliITSIOTrack *iotrack=new AliITSIOTrack;
 113  Int_t currentdir = 0;
 114  Int_t Ndirs;
 115
 116  if (fDirs)
 117   {
 118     Ndirs = fDirs->GetEntries();
 119   }
 120  else
 121   {
 122     Ndirs = 0;
 123   }
 124
 125  do //do while is good even if
 126   {
 127    TFile* gAliceFile = OpenGAliceFile(currentdir);
 128    if(gAliceFile == 0x0)
 129     {
 130        Error("Read","Can not open the file with gAlice");
 131        delete iotrack;
 132        return 1;
 133     }
 134    TFile *file = OpenTrackFile(currentdir);
 135    if(file == 0x0)
 136     {
 137        Error("Read","Can not open the file with ITS tracks V1");
 138        delete iotrack;
 139        return 2;
 140     }
 141
 142    if (gAlice->TreeE())//check if tree E exists
 143      {
 144       Nevents = (Int_t)gAlice->TreeE()->GetEntries();//if yes get number of events in gAlice
 145       cout<<"________________________________________________________\n";
 146       cout<<"Found "<<Nevents<<" event(s) in directory "<<GetDirName(currentdir)<<endl;
 147       cout<<"Setting Magnetic Field. Factor is "<<gAlice->Field()->Factor()<<endl;
 148       AliKalmanTrack::SetConvConst(100/0.299792458/0.2/gAlice->Field()->Factor());
 149      }
 150     else
 151      {//if not return an error
 152        Error("Read","Can not find Header tree (TreeE) in gAlice");
 153        delete iotrack;
 154        return 4;
 155      }
 156
 157    char tname[30];
 158    Int_t totalNevents = 0;
 159    for (Int_t currentEvent = 0; currentEvent < Nevents; currentEvent++)
 160     {
 161       cout<<"Reading Event "<<currentEvent<<endl;
 162
 163       sprintf(tname,"TreeT%d",currentEvent);
 164       file->cd();
 165       TTree *tracktree=(TTree*)file->Get(tname);
 166       TBranch *tbranch=tracktree->GetBranch("ITStracks");
 167       tbranch->SetAddress(&iotrack);
 168
 169       gAliceFile->cd();
 170       gAlice->GetEvent(currentEvent);
 171       gAlice->Particles();
 172
 173       Int_t nentr=(Int_t)tracktree->GetEntries();
 174
 175       for (Int_t i=0; i<nentr; i++)
 176        {
 177
 178         tracktree->GetEvent(i);
 179         if(!iotrack) continue;
 180         Int_t label = iotrack->GetLabel();
 181         if (label < 0)
 182          {
 183            continue;
 184            delete iotrack;
 185          }
 186         TParticle *p = (TParticle*)gAlice->Particle(label);
 187         if(!p)
 188          {
 189            Warning("Read","Can not get particle with label &d",label);
 190            continue;
 191          }
 192         if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type
 193                                            //if not take next partilce
 194
 195         AliHBTParticle* part = new AliHBTParticle(*p);
 196         if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
 197                                                 //if it does not delete it and take next good track
 198
 199         Double_t px=iotrack->GetPx();
 200         Double_t py=iotrack->GetPy();
 201         Double_t pz=iotrack->GetPz();
 202         Double_t mass = p->GetMass();
 203         Double_t tEtot = TMath::Sqrt(px*px + py*py + pz*pz + mass*mass);//total energy of the track
 204
 205         Double_t x= iotrack->GetX();
 206         Double_t y= iotrack->GetY();
 207         Double_t z= iotrack->GetZ();
 208
 209         AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(), px, py , pz, tEtot, x, y, z, 0.);
 210         if(Pass(track)) { delete  track;continue;}//check if meets all criteria of any of our cuts
 211                                                   //if it does not delete it and take next good track
 212
 213         particles->AddParticle(totalNevents,part);//put track and particle on the run
 214         tracks->AddParticle(totalNevents,track);
 215        }//end loop over tracks in the event
 216
 217        totalNevents++;
 218
 219      }//end of loop over events in current directory
 220
 221     gAliceFile->Close();
 222     delete gAliceFile;
 223     gAliceFile = 0;
 224
 225     file->Close();
 226     delete file;
 227     file = 0;
 228     currentdir++;
 229    }while(currentdir < Ndirs);//end of loop over directories specified in fDirs Obj Array
 230
 231
 232   delete iotrack;
 233   fIsRead = kTRUE;
 234   return 0;
 235
 236  }
 237 /********************************************************************/
 238
 239 TFile* AliHBTReaderITSv1::OpenTrackFile(Int_t ndir)
 240 {
 241 //opens files to be read for given directoru nomber in fDirs Array
 242    const TString& dirname = GetDirName(ndir);
 243    if (dirname == "")
 244     {
 245       Error("OpenGAliceFile","Can not get directory name");
 246       return 0x0;
 247     }
 248    TString filename = dirname + "/" + fITSTracksFileName;
 249    TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(filename.Data());
 250    if (!file) file = new TFile(filename.Data());
 251
 252    if (!file)
 253     {
 254       Error("Read","Can not open file %s",filename.Data());
 255       return 0x0;
 256     }
 257    if (!file->IsOpen())
 258     {
 259       Error("Read","Can not open file %s",filename.Data());
 260       return 0x0;
 261     }
 262
 263    return file;
 264 }
 265
 266
 267 /********************************************************************/
 268 TFile* AliHBTReaderITSv1::OpenGAliceFile(Int_t ndir)
 269 {
 270   const TString& dirname = GetDirName(ndir);
 271    if (dirname == "")
 272     {
 273       Error("OpenGAliceFile","Can not get directory name");
 274       return 0x0;
 275     }
 276
 277   TString filename = dirname + "/" + fGAliceFileName;
 278
 279   TFile* gAliceFile = TFile::Open(filename.Data());
 280   if ( gAliceFile== 0x0)
 281    {
 282      Error("OpenFiles","Can't open file named %s",filename.Data());
 283      return 0x0;
 284    }
 285   if (!gAliceFile->IsOpen())
 286    {
 287      Error("OpenFiles","Can't open file named %s",filename.Data());
 288      return 0x0;
 289    }
 290
 291   if (!(gAlice=(AliRun*)gAliceFile->Get("gAlice")))
 292    {
 293      Error("OpenFiles","gAlice have not been found on %s !\n",filename.Data());
 294      gAliceFile->Close();
 295      delete gAliceFile;
 296      return 0x0;
 297    }
 298
 299   return gAliceFile;
 300 }
 301
 302 /********************************************************************/
 303 /********************************************************************/
 304
 305