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  Int_t totalNevents = 0;
 116
 117  if (fDirs)
 118   {
 119     Ndirs = fDirs->GetEntries();
 120   }
 121  else
 122   {
 123     Ndirs = 0;
 124   }
 125
 126  do //do while is good even if
 127   {
 128    TFile* gAliceFile = OpenGAliceFile(currentdir);
 129    if(gAliceFile == 0x0)
 130     {
 131        Error("Read","Can not open the file with gAlice");
 132        delete iotrack;
 133        return 1;
 134     }
 135    if (gAlice->TreeE())//check if tree E exists
 136      {
 137       Nevents = (Int_t)gAlice->TreeE()->GetEntries();//if yes get number of events in gAlice
 138       cout<<"________________________________________________________\n";
 139       cout<<"Found "<<Nevents<<" event(s) in directory "<<GetDirName(currentdir)<<endl;
 140       cout<<"Setting Magnetic Field. Factor is "<<gAlice->Field()->Factor()<<endl;
 141       AliKalmanTrack::SetConvConst(100/0.299792458/0.2/gAlice->Field()->Factor());
 142      }
 143     else
 144      {//if not return an error
 145        Error("Read","Can not find Header tree (TreeE) in gAlice");
 146        delete iotrack;
 147        return 4;
 148      }
 149
 150    TFile *file = OpenTrackFile(currentdir);
 151    if(file == 0x0)
 152     {
 153        Error("Read","Can not open the file with ITS tracks V1");
 154        delete iotrack;
 155        return 2;
 156     }
 157
 158    Int_t naccepted = 0;
 159    char tname[30];
 160
 161    for (Int_t currentEvent = 0; currentEvent < Nevents; currentEvent++)
 162     {
 163       cout<<"Reading Event "<<currentEvent;
 164
 165       sprintf(tname,"TreeT%d",currentEvent);
 166       file->cd();
 167       TTree *tracktree=(TTree*)file->Get(tname);
 168       TBranch *tbranch=tracktree->GetBranch("ITStracks");
 169       tbranch->SetAddress(&iotrack);
 170
 171       gAliceFile->cd();
 172       gAlice->GetEvent(currentEvent);
 173       gAlice->Particles();
 174
 175       Int_t nentr=(Int_t)tracktree->GetEntries();
 176
 177       cout<<".  Found "<<nentr<<" tracks.";
 178       fflush(0);
 179
 180       for (Int_t i=0; i<nentr; i++)
 181        {
 182
 183         tracktree->GetEvent(i);
 184         if(!iotrack) continue;
 185         Int_t label = iotrack->GetLabel();
 186         if (label < 0)
 187          {
 188            continue;
 189          }
 190
 191         TParticle *p = (TParticle*)gAlice->Particle(label);
 192         if(!p)
 193          {
 194            Warning("Read","Can not get particle with label &d",label);
 195            continue;
 196          }
 197         if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type
 198                                            //if not take next partilce
 199
 200         AliHBTParticle* part = new AliHBTParticle(*p);
 201         if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
 202                                                 //if it does not delete it and take next good track
 203
 204         Double_t px=iotrack->GetPx();
 205         Double_t py=iotrack->GetPy();
 206         Double_t pz=iotrack->GetPz();
 207         Double_t mass = p->GetMass();
 208         Double_t tEtot = TMath::Sqrt(px*px + py*py + pz*pz + mass*mass);//total energy of the track
 209
 210         Double_t x= iotrack->GetX();
 211         Double_t y= iotrack->GetY();
 212         Double_t z= iotrack->GetZ();
 213
 214         AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(), px, py , pz, tEtot, x, y, z, 0.);
 215         if(Pass(track)) { delete  track;continue;}//check if meets all criteria of any of our cuts
 216                                                   //if it does not delete it and take next good track
 217
 218         particles->AddParticle(totalNevents,part);//put track and particle on the run
 219         tracks->AddParticle(totalNevents,track);
 220         naccepted++;
 221        }//end loop over tracks in the event
 222
 223        totalNevents++;
 224        cout<<"  Accepted "<<naccepted<<" tracks"<<endl;
 225      }//end of loop over events in current directory
 226
 227     gAliceFile->Close();
 228     delete gAliceFile;
 229     gAliceFile = 0;
 230
 231     file->Close();
 232     delete file;
 233     file = 0;
 234     currentdir++;
 235    }while(currentdir < Ndirs);//end of loop over directories specified in fDirs Obj Array
 236
 237
 238   delete iotrack;
 239   fIsRead = kTRUE;
 240   return 0;
 241
 242  }
 243 /********************************************************************/
 244
 245 TFile* AliHBTReaderITSv1::OpenTrackFile(Int_t ndir)
 246 {
 247 //opens files to be read for given directoru nomber in fDirs Array
 248    const TString& dirname = GetDirName(ndir);
 249    if (dirname == "")
 250     {
 251       Error("OpenGAliceFile","Can not get directory name");
 252       return 0x0;
 253     }
 254    TString filename = dirname + "/" + fITSTracksFileName;
 255
 256    TFile *file = TFile::Open(filename.Data());
 257    if (!file)
 258     {
 259       Error("Read","Can not open file %s",filename.Data());
 260       return 0x0;
 261     }
 262    if (!file->IsOpen())
 263     {
 264       Error("Read","Can not open file %s",filename.Data());
 265       return 0x0;
 266     }
 267
 268    return file;
 269 }
 270
 271
 272 /********************************************************************/
 273 TFile* AliHBTReaderITSv1::OpenGAliceFile(Int_t ndir)
 274 {
 275   const TString& dirname = GetDirName(ndir);
 276    if (dirname == "")
 277     {
 278       Error("OpenGAliceFile","Can not get directory name");
 279       return 0x0;
 280     }
 281
 282   TString filename = dirname + "/" + fGAliceFileName;
 283
 284   TFile* gAliceFile = TFile::Open(filename.Data());
 285   if ( gAliceFile== 0x0)
 286    {
 287      Error("OpenFiles","Can't open file named %s",filename.Data());
 288      return 0x0;
 289    }
 290   if (!gAliceFile->IsOpen())
 291    {
 292      Error("OpenFiles","Can't open file named %s",filename.Data());
 293      return 0x0;
 294    }
 295
 296   if (!(gAlice=(AliRun*)gAliceFile->Get("gAlice")))
 297    {
 298      Error("OpenFiles","gAlice have not been found on %s !\n",filename.Data());
 299      gAliceFile->Close();
 300      delete gAliceFile;
 301      return 0x0;
 302    }
 303
 304   return gAliceFile;
 305 }
 306
 307 /********************************************************************/
 308 /********************************************************************/
 309
 310