]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - HBTAN/AliHBTReaderITSv2.cxx
Transition to NewIO
[u/mrichter/AliRoot.git] / HBTAN / AliHBTReaderITSv2.cxx
index 7867873a5ae5acaa60063436cf5e866d0f82abf2..dd56bda19c255b13845a40411766e33c512d8099 100644 (file)
@@ -1,8 +1,8 @@
 
 #include "AliHBTReaderITSv2.h"
 
-#include <iostream.h>
-#include <fstream.h>
+#include <Riostream.h>
+#include <Riostream.h>
 #include <TString.h>
 #include <TObjString.h>
 #include <TTree.h>
 #include <TParticle.h>
 
 #include <AliRun.h>
+#include <AliRunLoader.h>
+#include <AliLoader.h>
 #include <AliMagF.h>
+#include <AliITS.h>
 #include <AliITStrackV2.h>
-//#include <AliITSParam.h>
 #include <AliITStrackerV2.h>
 #include <AliITSgeom.h>
 
 
 ClassImp(AliHBTReaderITSv2)
 
-AliHBTReaderITSv2::
- AliHBTReaderITSv2(const Char_t* trackfilename, const Char_t* clusterfilename,
-       const Char_t* galicefilename)
-                  :fTrackFileName(trackfilename),fClusterFileName(clusterfilename),
-                    fGAliceFileName(galicefilename)
+AliHBTReaderITSv2::AliHBTReaderITSv2():fFileName("galice.root")
+{
+  //constructor, 
+  //Defaults:
+  //  galicefilename = "galice.root"
+  fParticles = 0x0;
+  fTracks    = 0x0;
+  fIsRead = kFALSE;
+}
+/********************************************************************/
+
+AliHBTReaderITSv2::AliHBTReaderITSv2(const Char_t* galicefilename):fFileName(galicefilename)
 {
   //constructor, 
   //Defaults:
-  //  trackfilename = "AliITStracksV2.root"
-  //  clusterfilename = "AliITSclustersV2.root"
   //  galicefilename = "galice.root"
   fParticles = new AliHBTRun();
   fTracks    = new AliHBTRun();
@@ -41,17 +48,11 @@ AliHBTReaderITSv2::
 }
 /********************************************************************/
 
-AliHBTReaderITSv2::
- AliHBTReaderITSv2(TObjArray* dirs, const Char_t* trackfilename, 
-                   const Char_t* clusterfilename, const Char_t* galicefilename)
-                  : AliHBTReader(dirs),
-                    fTrackFileName(trackfilename),fClusterFileName(clusterfilename),
-                    fGAliceFileName(galicefilename)
+AliHBTReaderITSv2::AliHBTReaderITSv2(TObjArray* dirs, const Char_t* galicefilename): 
+       AliHBTReader(dirs), fFileName(galicefilename)
 {
   //constructor, 
   //Defaults:
-  //  trackfilename = "AliITStracksV2.root"
-  //  clusterfilename = "AliITSclustersV2.root"
   //  galicefilename = "galice.root"
   
   fParticles = new AliHBTRun();
@@ -71,40 +72,51 @@ AliHBTReaderITSv2::~AliHBTReaderITSv2()
 AliHBTEvent* AliHBTReaderITSv2::GetParticleEvent(Int_t n)
  {
  //returns Nth event with simulated particles
-   if (!fIsRead) 
+ if (!fIsRead) 
+  { 
+    if (fParticles == 0x0) fParticles = new AliHBTRun();
+    if (fTracks == 0x0) fTracks    = new AliHBTRun();
     if(Read(fParticles,fTracks))
      {
        Error("GetParticleEvent","Error in reading");
        return 0x0;
      }
-
  return fParticles->GetEvent(n);
- }
+  }
return (fParticles)?fParticles->GetEvent(n):0x0;
+}
 /********************************************************************/
 
 AliHBTEvent* AliHBTReaderITSv2::GetTrackEvent(Int_t n)
  {
  //returns Nth event with reconstructed tracks
-   if (!fIsRead) 
+ if (!fIsRead) 
+  { 
+    if (fParticles == 0x0) fParticles = new AliHBTRun();
+    if (fTracks == 0x0) fTracks    = new AliHBTRun();
     if(Read(fParticles,fTracks))
      {
        Error("GetTrackEvent","Error in reading");
        return 0x0;
      }
-   return fTracks->GetEvent(n);
+   }
+  return (fTracks)?fTracks->GetEvent(n):0x0;
  }
 /********************************************************************/
 
 Int_t AliHBTReaderITSv2::GetNumberOfPartEvents()
  {
  //returns number of events of particles
-   if (!fIsRead)
+  if (!fIsRead)
+   {
+    if (fParticles == 0x0) fParticles = new AliHBTRun();
+    if (fTracks == 0x0) fTracks    = new AliHBTRun();
     if(Read(fParticles,fTracks))
      {
        Error("GetNumberOfPartEvents","Error in reading");
        return 0;
      }
-   return fParticles->GetNumberOfEvents();
+   }
+  return (fParticles)?fParticles->GetNumberOfEvents():0;
  }
 
 /********************************************************************/
@@ -112,30 +124,31 @@ Int_t AliHBTReaderITSv2::GetNumberOfTrackEvents()
  {
  //returns number of events of tracks
   if (!fIsRead) 
+   {
+    if (fParticles == 0x0) fParticles = new AliHBTRun();
+    if (fTracks == 0x0) fTracks    = new AliHBTRun();
     if(Read(fParticles,fTracks))
      {
        Error("GetNumberOfTrackEvents","Error in reading");
        return 0;
      }
-  return fTracks->GetNumberOfEvents();
+   }
+  return (fTracks)?fTracks->GetNumberOfEvents():0;
  }
-
-
 /********************************************************************/
 /********************************************************************/
+
 Int_t AliHBTReaderITSv2::Read(AliHBTRun* particles, AliHBTRun *tracks)
 {
+//reads data
  Int_t Nevents = 0; //number of events found in given directory
  Int_t Ndirs; //number of the directories to be read
  Int_t Ntracks; //number of tracks in current event
  Int_t currentdir = 0; //number of events in the current directory 
  Int_t totalNevents = 0; //total number of events read from all directories up to now
- register Int_t i; //iterator
+ register Int_t i = 0; //iterator
  
- TFile *aTracksFile;//file with tracks
- TFile *aClustersFile;//file with clusters
- TFile *aGAliceFile;//file name with galice
-
 // AliITStrackerV2 *tracker; // ITS tracker - used for cooking labels
  TTree *tracktree; // tree for tracks
  
@@ -144,8 +157,7 @@ Int_t AliHBTReaderITSv2::Read(AliHBTRun* particles, AliHBTRun *tracks)
  Float_t phi, lam, pt;//angles and transverse momentum
  Int_t label; //label of the current track
 
- char tname[100]; //buffer for tree name
- AliITStrackV2 *iotrack= new AliITStrackV2(); //buffer track for reading data from tree
+ AliITStrackV2 *iotrack = 0x0; //buffer track for reading data from tree
 
  if (!particles) //check if an object is instatiated
   {
@@ -171,60 +183,91 @@ Int_t AliHBTReaderITSv2::Read(AliHBTRun* particles, AliHBTRun *tracks)
  
  do //do while is good even if Ndirs==0 (than read from current directory)
    {
-    if( (i=OpenFiles(aTracksFile,aClustersFile,aGAliceFile,currentdir)) )
+    TString filename = GetDirName(currentdir);
+    if (filename.IsNull())
+     {
+       Error("Read","Can not get directory name");
+       currentdir++;
+       continue;
+     }
+    filename = filename +"/"+ fFileName;
+    AliRunLoader* rl = AliRunLoader::Open(filename);
+    if( rl == 0x0)
      {
-       Error("Read","Exiting due to problems with opening files. Errorcode %d",i);
-       delete iotrack;
-       return i;
+       Error("Read","Exiting due to problems with opening files.");
+       currentdir++;
+       continue;
      }
     
-    if (gAlice->TreeE())//check if tree E exists
+    rl->LoadHeader();
+    rl->LoadKinematics();
+    rl->LoadgAlice();
+    gAlice = rl->GetAliRun();
+    AliITS* its = (AliITS*)gAlice->GetModule("ITS");
+    
+    AliLoader* itsl = rl->GetLoader("ITSLoader");
+    
+    if ((its == 0x0) || ( itsl== 0x0))
+     {
+       Error("Read","Can not found ITS in this run");
+       delete rl;
+       rl = 0x0;
+       currentdir++;
+       continue;
+     }
+    Nevents = rl->GetNumberOfEvents();
+    if (Nevents > 0)//check if tree E exists
      {
-      Nevents = (Int_t)gAlice->TreeE()->GetEntries();//if yes get number of events in gAlice
-      cout<<"________________________________________________________\n";
-      cout<<"Found "<<Nevents<<" event(s) in directory "<<GetDirName(currentdir)<<endl;
-      cout<<"Setting Magnetic Field. Factor is "<<gAlice->Field()->Factor()<<endl;
-      AliKalmanTrack::SetConvConst(100/0.299792458/0.2/gAlice->Field()->Factor());  
+      Info("Read","________________________________________________________");
+      Info("Read","Found %d event(s) in directory %s",Nevents,GetDirName(currentdir).Data());
+      Float_t mf;
+      if (fUseMagFFromRun)
+       {
+         mf = gAlice->Field()->SolenoidField();
+         Info("Read","Setting Magnetic Field from run: B=%fT",mf/10.);
+       }
+      else
+       {
+         Info("Read","Setting Own Magnetic Field: B=%fT",fMagneticField);
+         mf = fMagneticField*10.;
+       }
+      AliKalmanTrack::SetConvConst(1000/0.299792458/mf);
+      if (iotrack == 0x0) iotrack = new AliITStrackV2();
      }
     else
      {//if not return an error
-       Error("Read","Can not find Header tree (TreeE) in gAlice");
-       delete iotrack;
-       return 1;
+       Error("Read","No events in this run");
+       delete rl;
+       rl = 0x0;
+       currentdir++;
+       continue;
      }
     
-    AliITSgeom *geom=(AliITSgeom*)aClustersFile->Get("AliITSgeom");
+    AliITSgeom *geom= its->GetITSgeom();
     if (!geom) 
      { 
        Error("Read","Can't get the ITS geometry!"); 
-       delete iotrack;
-       return 2; 
+       delete rl;
+       rl = 0x0;
+       currentdir++;
+       continue;
      }
 
+    itsl->LoadTracks();
+
     for(Int_t currentEvent =0; currentEvent<Nevents;currentEvent++)//loop over all events
      {
        cout<<"Reading Event "<<currentEvent<<endl;
+       rl->GetEvent(currentEvent);
+       tracktree=itsl->TreeT();
        
-       aGAliceFile->cd();
-       gAlice->GetEvent(currentEvent);
-       TParticle * part = gAlice->Particle(0);
-       Double_t orz=part->Vz();
-      
-       aClustersFile->cd();
-//       tracker = new AliITStrackerV2(geom,currentEvent,orz); //<---- this is for Massimo version
-//       tracker = new AliITStrackerV2(geom,currentEvent);
-       sprintf(tname,"TreeT_ITS_%d",currentEvent);
-       
-       tracktree=(TTree*)aTracksFile->Get(tname);
        if (!tracktree) 
          {
            Error("Read","Can't get a tree with ITS tracks"); 
-           delete iotrack;
- //          delete tracker;
-           return 4;
+           continue;
          }
        TBranch *tbranch=tracktree->GetBranch("tracks");
-      
        Ntracks=(Int_t)tracktree->GetEntries();
 
        Int_t accepted = 0;
@@ -232,7 +275,7 @@ Int_t AliHBTReaderITSv2::Read(AliHBTRun* particles, AliHBTRun *tracks)
        Int_t itsfault = 0;
        for (i=0; i<Ntracks; i++) //loop over all tpc tracks
         { 
-          if(i%100 == 0)cout<<"all: "<<i<<"   accepted: "<<accepted<<"   tpc faults: "<<tpcfault<<"   its faults: "<<itsfault<<"\r";
+          if(i%100 == 0)cout<<"all: "<<i<<"   accepted: "<<accepted<<"   tpc faults: "<<tpcfault<<"\r";
           
           tbranch->SetAddress(&iotrack);
           tracktree->GetEvent(i);
@@ -243,24 +286,19 @@ Int_t AliHBTReaderITSv2::Read(AliHBTRun* particles, AliHBTRun *tracks)
              tpcfault++;
              continue;
            }
-//          tracker->CookLabel(iotrack,0.);
-//          Int_t itsLabel=iotrack->GetLabel();
- //         if (itsLabel != label)
- //          {
- //           itsfault++;
- //           continue; 
- //          }
 
           TParticle *p = (TParticle*)gAlice->Particle(label);
+          if(p == 0x0) continue; //if returned pointer is NULL
+          if(p->GetPDG() == 0x0) continue; //if particle has crezy PDG code (not known to our database)
+
           if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type 
                                               //if not take next partilce
             
-          AliHBTParticle* part = new AliHBTParticle(*p);
+          AliHBTParticle* part = new AliHBTParticle(*p,i);
           if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
                                                   //if it does not delete it and take next good track
 
-          
-          iotrack->Propagate(iotrack->GetAlpha(),3.,0.1/65.19*1.848,0.1*1.848);
+          iotrack->PropagateTo(3.,0.0028,65.19);
           iotrack->PropagateToVertex();
  
           iotrack->GetExternalParameters(xk,par);     //get properties of the track
@@ -277,7 +315,7 @@ Int_t AliHBTReaderITSv2::Read(AliHBTRun* particles, AliHBTRun *tracks)
           Double_t mass = p->GetMass();
           Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);//total energy of the track
             
-          AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(), tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
+          AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(), i, tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
           if(Pass(track))//check if meets all criteria of any of our cuts
                          //if it does not delete it and take next good track
            { 
@@ -289,17 +327,13 @@ Int_t AliHBTReaderITSv2::Read(AliHBTRun* particles, AliHBTRun *tracks)
           tracks->AddParticle(totalNevents,track);
           accepted++;
         }//end of loop over tracks in the event
-        
-       aTracksFile->Delete(tname);
-       aTracksFile->Delete("tracks");
-//       delete tracker;
        
        totalNevents++;
-       CloseFiles(aTracksFile,aClustersFile,aGAliceFile);     
        cout<<"all: "<<i<<"   accepted: "<<accepted<<"   tpc faults: "<<tpcfault<<"   its faults: "<<itsfault<<endl;
      
      }//end of loop over events in current directory
-     currentdir++;
+    delete rl;
+    currentdir++;
    }while(currentdir < Ndirs);//end of loop over directories specified in fDirs Obj Array
 
  delete iotrack;
@@ -308,92 +342,6 @@ Int_t AliHBTReaderITSv2::Read(AliHBTRun* particles, AliHBTRun *tracks)
 }
 
 /********************************************************************/
-Int_t AliHBTReaderITSv2::OpenFiles
-(TFile*& aTracksFile, TFile*& aClustersFile, TFile*& agAliceFile,Int_t event)
-{
- //opens all the files
-   
-   
-   const TString& dirname = GetDirName(event); 
-   if (dirname == "")
-    {
-      Error("OpenFiles","Can not get directory name");
-      return 4;
-    }
-   
-   TString filename = dirname +"/"+ fTrackFileName;
-   aTracksFile = TFile::Open(filename.Data());
-   if ( aTracksFile  == 0x0 ) 
-     {
-       Error("OpenFiles","Can't open file with tacks named %s",filename.Data());
-       return 1;
-     }
-   if (!aTracksFile->IsOpen())
-     {
-       Error("OpenFiles","Can't open file with tacks named %s",filename.Data());
-       return 1;
-     }
-  
-   filename = dirname +"/"+ fClusterFileName;
-   aClustersFile = TFile::Open(filename.Data());
-   if ( aClustersFile == 0x0 )
-    {
-      Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
-      return 2;
-    }
-   if (!aClustersFile->IsOpen())
-    {
-      Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
-      return 2;
-    }
-
-   filename = dirname +"/"+ fGAliceFileName;
-   agAliceFile = TFile::Open(filename.Data());
-   if ( agAliceFile== 0x0)
-    {
-      Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
-      return 3;
-    }
-   if (!agAliceFile->IsOpen())
-    {
-      Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
-      return 3;
-    } 
-   
-   if (!(gAlice=(AliRun*)agAliceFile->Get("gAlice"))) 
-    {
-      Error("OpenFiles","gAlice have not been found on %s !\n",filename.Data());
-      return 5;
-    }
-
-   return 0; 
-}
-/********************************************************************/
-
-/********************************************************************/
-  
-void AliHBTReaderITSv2::CloseFiles(TFile*& tracksFile, TFile*& clustersFile, TFile*& gAliceFile)
-{
-  //closes the files
-  tracksFile->Close();
-  delete tracksFile;
-  tracksFile = 0x0;
-  
-  clustersFile->Close();
-  delete clustersFile;
-  clustersFile = 0x0;
-  
-  delete gAlice;
-  gAlice = 0;
-
-  if (gAliceFile) 
-   {
-     gAliceFile->Close();
-     delete gAliceFile;
-     gAliceFile = 0x0;
-   }
-}
-
 /********************************************************************/