HBTAN/AliHBTReaderITSv2.cxx

   1
   2 #include "AliHBTReaderITSv2.h"
   3
   4 #include <iostream.h>
   5 #include <fstream.h>
   6 #include <TString.h>
   7 #include <TObjString.h>
   8 #include <TTree.h>
   9 #include <TFile.h>
  10 #include <TParticle.h>
  11
  12 #include <AliRun.h>
  13 #include <AliMagF.h>
  14 #include <AliITStrackV2.h>
  15 //#include <AliITSParam.h>
  16 #include <AliITStrackerV2.h>
  17 #include <AliITSgeom.h>
  18
  19 #include "AliHBTRun.h"
  20 #include "AliHBTEvent.h"
  21 #include "AliHBTParticle.h"
  22 #include "AliHBTParticleCut.h"
  23
  24
  25 ClassImp(AliHBTReaderITSv2)
  26
  27 AliHBTReaderITSv2::
  28  AliHBTReaderITSv2(const Char_t* trackfilename, const Char_t* clusterfilename,
  29         const Char_t* galicefilename)
  30                   :fTrackFileName(trackfilename),fClusterFileName(clusterfilename),
  31                     fGAliceFileName(galicefilename)
  32 {
  33   //constructor,
  34   //Defaults:
  35   //  trackfilename = "AliITStracksV2.root"
  36   //  clusterfilename = "AliITSclustersV2.root"
  37   //  galicefilename = "galice.root"
  38   fParticles = new AliHBTRun();
  39   fTracks    = new AliHBTRun();
  40   fIsRead = kFALSE;
  41 }
  42 /********************************************************************/
  43
  44 AliHBTReaderITSv2::
  45  AliHBTReaderITSv2(TObjArray* dirs, const Char_t* trackfilename,
  46                    const Char_t* clusterfilename, const Char_t* galicefilename)
  47                   : AliHBTReader(dirs),
  48                     fTrackFileName(trackfilename),fClusterFileName(clusterfilename),
  49                     fGAliceFileName(galicefilename)
  50 {
  51   //constructor,
  52   //Defaults:
  53   //  trackfilename = "AliITStracksV2.root"
  54   //  clusterfilename = "AliITSclustersV2.root"
  55   //  galicefilename = "galice.root"
  56
  57   fParticles = new AliHBTRun();
  58   fTracks    = new AliHBTRun();
  59   fIsRead = kFALSE;
  60 }
  61 /********************************************************************/
  62
  63 AliHBTReaderITSv2::~AliHBTReaderITSv2()
  64  {
  65    if (fParticles) delete fParticles;
  66    if (fTracks) delete fTracks;
  67  }
  68 /********************************************************************/
  69 /********************************************************************/
  70
  71 AliHBTEvent* AliHBTReaderITSv2::GetParticleEvent(Int_t n)
  72  {
  73  //returns Nth event with simulated particles
  74    if (!fIsRead)
  75     if(Read(fParticles,fTracks))
  76      {
  77        Error("GetParticleEvent","Error in reading");
  78        return 0x0;
  79      }
  80
  81    return fParticles->GetEvent(n);
  82  }
  83 /********************************************************************/
  84
  85 AliHBTEvent* AliHBTReaderITSv2::GetTrackEvent(Int_t n)
  86  {
  87  //returns Nth event with reconstructed tracks
  88    if (!fIsRead)
  89     if(Read(fParticles,fTracks))
  90      {
  91        Error("GetTrackEvent","Error in reading");
  92        return 0x0;
  93      }
  94    return fTracks->GetEvent(n);
  95  }
  96 /********************************************************************/
  97
  98 Int_t AliHBTReaderITSv2::GetNumberOfPartEvents()
  99  {
 100  //returns number of events of particles
 101    if (!fIsRead)
 102     if(Read(fParticles,fTracks))
 103      {
 104        Error("GetNumberOfPartEvents","Error in reading");
 105        return 0;
 106      }
 107    return fParticles->GetNumberOfEvents();
 108  }
 109
 110 /********************************************************************/
 111 Int_t AliHBTReaderITSv2::GetNumberOfTrackEvents()
 112  {
 113  //returns number of events of tracks
 114   if (!fIsRead)
 115     if(Read(fParticles,fTracks))
 116      {
 117        Error("GetNumberOfTrackEvents","Error in reading");
 118        return 0;
 119      }
 120   return fTracks->GetNumberOfEvents();
 121  }
 122
 123
 124 /********************************************************************/
 125 /********************************************************************/
 126 Int_t AliHBTReaderITSv2::Read(AliHBTRun* particles, AliHBTRun *tracks)
 127 {
 128  Int_t Nevents = 0; //number of events found in given directory
 129  Int_t Ndirs; //number of the directories to be read
 130  Int_t Ntracks; //number of tracks in current event
 131  Int_t currentdir = 0; //number of events in the current directory
 132  Int_t totalNevents = 0; //total number of events read from all directories up to now
 133  register Int_t i; //iterator
 134
 135  TFile *aTracksFile;//file with tracks
 136  TFile *aClustersFile;//file with clusters
 137  TFile *aGAliceFile;//file name with galice
 138
 139 // AliITStrackerV2 *tracker; // ITS tracker - used for cooking labels
 140  TTree *tracktree; // tree for tracks
 141
 142  Double_t xk;
 143  Double_t par[5]; //Kalman track parameters
 144  Float_t phi, lam, pt;//angles and transverse momentum
 145  Int_t label; //label of the current track
 146
 147  char tname[100]; //buffer for tree name
 148  AliITStrackV2 *iotrack= new AliITStrackV2(); //buffer track for reading data from tree
 149
 150  if (!particles) //check if an object is instatiated
 151   {
 152     Error("Read"," particles object must instatiated before passing it to the reader");
 153   }
 154  if (!tracks)  //check if an object is instatiated
 155   {
 156     Error("Read"," tracks object must instatiated before passing it to the reader");
 157   }
 158  particles->Reset();//clear runs == delete all old events
 159  tracks->Reset();
 160
 161  if (fDirs) //if array with directories is supplied by user
 162   {
 163     Ndirs = fDirs->GetEntries(); //get the number if directories
 164   }
 165  else
 166   {
 167     Ndirs = 0; //if the array is not supplied read only from current directory
 168   }
 169
 170 // cout<<"Found "<<Ndirs<<" directory entries"<<endl;
 171
 172  do //do while is good even if Ndirs==0 (than read from current directory)
 173    {
 174     if( (i=OpenFiles(aTracksFile,aClustersFile,aGAliceFile,currentdir)) )
 175      {
 176        Error("Read","Exiting due to problems with opening files. Errorcode %d",i);
 177        delete iotrack;
 178        return i;
 179      }
 180
 181     if (gAlice->TreeE())//check if tree E exists
 182      {
 183       Nevents = (Int_t)gAlice->TreeE()->GetEntries();//if yes get number of events in gAlice
 184       cout<<"________________________________________________________\n";
 185       cout<<"Found "<<Nevents<<" event(s) in directory "<<GetDirName(currentdir)<<endl;
 186       cout<<"Setting Magnetic Field. Factor is "<<gAlice->Field()->Factor()<<endl;
 187       AliKalmanTrack::SetConvConst(100/0.299792458/0.2/gAlice->Field()->Factor());
 188      }
 189     else
 190      {//if not return an error
 191        Error("Read","Can not find Header tree (TreeE) in gAlice");
 192        delete iotrack;
 193        return 1;
 194      }
 195
 196     AliITSgeom *geom=(AliITSgeom*)aClustersFile->Get("AliITSgeom");
 197     if (!geom)
 198      {
 199        Error("Read","Can't get the ITS geometry!");
 200        delete iotrack;
 201        return 2;
 202      }
 203
 204     for(Int_t currentEvent =0; currentEvent<Nevents;currentEvent++)//loop over all events
 205      {
 206        cout<<"Reading Event "<<currentEvent<<endl;
 207
 208        aGAliceFile->cd();
 209        gAlice->GetEvent(currentEvent);
 210        TParticle * part = gAlice->Particle(0);
 211        Double_t orz=part->Vz();
 212
 213        aClustersFile->cd();
 214 //       tracker = new AliITStrackerV2(geom,currentEvent,orz); //<---- this is for Massimo version
 215 //       tracker = new AliITStrackerV2(geom,currentEvent);
 216        sprintf(tname,"TreeT_ITS_%d",currentEvent);
 217
 218        tracktree=(TTree*)aTracksFile->Get(tname);
 219        if (!tracktree)
 220          {
 221            Error("Read","Can't get a tree with ITS tracks");
 222            delete iotrack;
 223  //          delete tracker;
 224            return 4;
 225          }
 226        TBranch *tbranch=tracktree->GetBranch("tracks");
 227
 228        Ntracks=(Int_t)tracktree->GetEntries();
 229
 230        Int_t accepted = 0;
 231        Int_t tpcfault = 0;
 232        Int_t itsfault = 0;
 233        for (i=0; i<Ntracks; i++) //loop over all tpc tracks
 234         {
 235           if(i%100 == 0)cout<<"all: "<<i<<"   accepted: "<<accepted<<"   tpc faults: "<<tpcfault<<"   its faults: "<<itsfault<<"\r";
 236
 237           tbranch->SetAddress(&iotrack);
 238           tracktree->GetEvent(i);
 239
 240           label=iotrack->GetLabel();
 241           if (label < 0)
 242            {
 243              tpcfault++;
 244              continue;
 245            }
 246 //          tracker->CookLabel(iotrack,0.);
 247 //          Int_t itsLabel=iotrack->GetLabel();
 248  //         if (itsLabel != label)
 249  //          {
 250  //           itsfault++;
 251  //           continue;
 252  //          }
 253
 254           TParticle *p = (TParticle*)gAlice->Particle(label);
 255           if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type
 256                                               //if not take next partilce
 257
 258           AliHBTParticle* part = new AliHBTParticle(*p);
 259           if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
 260                                                   //if it does not delete it and take next good track
 261
 262
 263           iotrack->Propagate(iotrack->GetAlpha(),3.,0.1/65.19*1.848,0.1*1.848);
 264           iotrack->PropagateToVertex();
 265
 266           iotrack->GetExternalParameters(xk,par);     //get properties of the track
 267           phi=TMath::ASin(par[2]) + iotrack->GetAlpha();
 268           if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
 269           if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
 270           lam=par[3];
 271           pt=1.0/TMath::Abs(par[4]);
 272
 273           Double_t tpx = pt * TMath::Cos(phi); //track x coordinate of momentum
 274           Double_t tpy = pt * TMath::Sin(phi); //track y coordinate of momentum
 275           Double_t tpz = pt * lam; //track z coordinate of momentum
 276
 277           Double_t mass = p->GetMass();
 278           Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);//total energy of the track
 279
 280           AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(), tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
 281           if(Pass(track))//check if meets all criteria of any of our cuts
 282                          //if it does not delete it and take next good track
 283            {
 284             delete track;
 285             delete part;
 286             continue;
 287            }
 288           particles->AddParticle(totalNevents,part);//put track and particle on the run
 289           tracks->AddParticle(totalNevents,track);
 290           accepted++;
 291         }//end of loop over tracks in the event
 292
 293        aTracksFile->Delete(tname);
 294        aTracksFile->Delete("tracks");
 295 //       delete tracker;
 296
 297        totalNevents++;
 298        CloseFiles(aTracksFile,aClustersFile,aGAliceFile);
 299        cout<<"all: "<<i<<"   accepted: "<<accepted<<"   tpc faults: "<<tpcfault<<"   its faults: "<<itsfault<<endl;
 300
 301      }//end of loop over events in current directory
 302      currentdir++;
 303    }while(currentdir < Ndirs);//end of loop over directories specified in fDirs Obj Array
 304
 305  delete iotrack;
 306  fIsRead = kTRUE;
 307  return 0;
 308 }
 309
 310 /********************************************************************/
 311 Int_t AliHBTReaderITSv2::OpenFiles
 312 (TFile*& aTracksFile, TFile*& aClustersFile, TFile*& agAliceFile,Int_t event)
 313 {
 314  //opens all the files
 315
 316
 317    const TString& dirname = GetDirName(event);
 318    if (dirname == "")
 319     {
 320       Error("OpenFiles","Can not get directory name");
 321       return 4;
 322     }
 323
 324    TString filename = dirname +"/"+ fTrackFileName;
 325    aTracksFile = TFile::Open(filename.Data());
 326    if ( aTracksFile  == 0x0 )
 327      {
 328        Error("OpenFiles","Can't open file with tacks named %s",filename.Data());
 329        return 1;
 330      }
 331    if (!aTracksFile->IsOpen())
 332      {
 333        Error("OpenFiles","Can't open file with tacks named %s",filename.Data());
 334        return 1;
 335      }
 336
 337    filename = dirname +"/"+ fClusterFileName;
 338    aClustersFile = TFile::Open(filename.Data());
 339    if ( aClustersFile == 0x0 )
 340     {
 341       Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
 342       return 2;
 343     }
 344    if (!aClustersFile->IsOpen())
 345     {
 346       Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
 347       return 2;
 348     }
 349
 350    filename = dirname +"/"+ fGAliceFileName;
 351    agAliceFile = TFile::Open(filename.Data());
 352    if ( agAliceFile== 0x0)
 353     {
 354       Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
 355       return 3;
 356     }
 357    if (!agAliceFile->IsOpen())
 358     {
 359       Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
 360       return 3;
 361     }
 362
 363    if (!(gAlice=(AliRun*)agAliceFile->Get("gAlice")))
 364     {
 365       Error("OpenFiles","gAlice have not been found on %s !\n",filename.Data());
 366       return 5;
 367     }
 368
 369    return 0;
 370 }
 371 /********************************************************************/
 372
 373 /********************************************************************/
 374
 375 void AliHBTReaderITSv2::CloseFiles(TFile*& tracksFile, TFile*& clustersFile, TFile*& gAliceFile)
 376 {
 377   //closes the files
 378   tracksFile->Close();
 379   delete tracksFile;
 380   tracksFile = 0x0;
 381
 382   clustersFile->Close();
 383   delete clustersFile;
 384   clustersFile = 0x0;
 385
 386   delete gAlice;
 387   gAlice = 0;
 388
 389   if (gAliceFile)
 390    {
 391      gAliceFile->Close();
 392      delete gAliceFile;
 393      gAliceFile = 0x0;
 394    }
 395 }
 396
 397 /********************************************************************/
 398
 399