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