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