HBTAN/AliHBTReaderITSv1.cxx

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