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