]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - HBTAN/AliHBTReaderTPC.cxx
Transition to NewIO
[u/mrichter/AliRoot.git] / HBTAN / AliHBTReaderTPC.cxx
index aee2bd74e92b98ba3ce16cced1518924f33f846a..92dbd82ae3d227dce5aa33c56fa5acdd5be2ce5c 100644 (file)
@@ -1,13 +1,16 @@
 #include "AliHBTReaderTPC.h"
+
 #include <TTree.h>
 #include <TFile.h>
 #include <TParticle.h>
 
 #include <AliRun.h>
+#include <AliLoader.h>
+#include <AliStack.h>
 #include <AliMagF.h>
 #include <AliTPCtrack.h>
 #include <AliTPCParam.h>
+#include <AliTPCtracker.h>
 
 #include "AliHBTRun.h"
 #include "AliHBTEvent.h"
@@ -21,56 +24,46 @@ ClassImp(AliHBTReaderTPC)
 // class AliHBTReaderTPC
 //
 //reader for TPC tracking
-//needs galice.root, AliTPCtracks.root
+//needs galice.root, AliTPCtracks.root, AliTPCclusters.root
 //
 //more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html
 //Piotr.Skowronski@cern.ch
-//
-
+AliHBTReaderTPC::AliHBTReaderTPC():fFileName("galice.root")
+{
+  //constructor, 
+  //Defaults:
+  //  galicefilename = ""  - this means: Do not open gAlice file - 
+  //                         just leave the global pointer untouched
+  
+  fParticles = 0x0;
+  fTracks    = 0x0;
+  fIsRead = kFALSE;
+}
 
-AliHBTReaderTPC:: AliHBTReaderTPC(const Char_t* trackfilename,
-                                  const Char_t* clusterfilename,
-                                  const Char_t* galicefilename):
- fParticles(new AliHBTRun()),
- fTracks(new AliHBTRun()),
- fTrackFileName(trackfilename),
- fClusterFileName(clusterfilename),
- fGAliceFileName(galicefilename),
- fIsRead(kFALSE),
- fMagneticField(0.4),
- fUseMagFFromRun(kTRUE)
+AliHBTReaderTPC::AliHBTReaderTPC(const Char_t* galicefilename):fFileName(galicefilename)
 {
   //constructor, 
   //Defaults:
-  //  trackfilename = "AliTPCtracks.root"
-  //  clusterfilename = "AliTPCclusters.root"
   //  galicefilename = ""  - this means: Do not open gAlice file - 
-  //                         just leave the global pointer untached
-  //this class is not supposed to be written, that is why I let myself and user
-  //for comfort of having default constructor which allocates dynamic memmory
-  //in case it is going to be written, it should be changed
+  //                         just leave the global pointer untouched
   
+  fParticles = new AliHBTRun();
+  fTracks    = new AliHBTRun();
+  fIsRead = kFALSE;
 }
 /********************************************************************/
-AliHBTReaderTPC::AliHBTReaderTPC(TObjArray* dirs,
-                  const Char_t* trackfilename, const Char_t* clusterfilename,
-                  const Char_t* galicefilename):
- AliHBTReader(dirs), 
- fParticles(new AliHBTRun()),
- fTracks(new AliHBTRun()),
- fTrackFileName(trackfilename),
- fClusterFileName(clusterfilename),
- fGAliceFileName(galicefilename),
- fIsRead(kFALSE),
- fMagneticField(0.4),
- fUseMagFFromRun(kTRUE)
+AliHBTReaderTPC::AliHBTReaderTPC(TObjArray* dirs, const Char_t* galicefilename):
+                  AliHBTReader(dirs), fFileName(galicefilename)
+
 {
   //constructor, 
   //Defaults:
-  //  trackfilename = "AliTPCtracks.root"
-  //  clusterfilename = "AliTPCclusters.root"
   //  galicefilename = ""  - this means: Do not open gAlice file - 
   //                         just leave the global pointer untached
+  
+  fParticles = new AliHBTRun();
+  fTracks    = new AliHBTRun();
+  fIsRead = kFALSE;
 }
 /********************************************************************/
 
@@ -85,13 +78,18 @@ AliHBTReaderTPC::~AliHBTReaderTPC()
 AliHBTEvent* AliHBTReaderTPC::GetParticleEvent(Int_t n)
  {
  //returns Nth event with simulated particles
-   if (!fIsRead) 
-    if(Read(fParticles,fTracks))
-     {
-       Error("GetParticleEvent","Error in reading");
-       return 0x0;
-     }
-   return fParticles->GetEvent(n);
+   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)?fParticles->GetEvent(n):0x0;
  }
 /********************************************************************/
 
@@ -99,12 +97,16 @@ AliHBTEvent* AliHBTReaderTPC::GetTrackEvent(Int_t n)
  {
  //returns Nth event with reconstructed tracks
    if (!fIsRead) 
-    if(Read(fParticles,fTracks))
-     {
-       Error("GetTrackEvent","Error in reading");
-       return 0x0;
-     }
-   return fTracks->GetEvent(n);
+    {
+      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)?fTracks->GetEvent(n):0x0;
  }
 /********************************************************************/
 
@@ -112,12 +114,16 @@ Int_t AliHBTReaderTPC::GetNumberOfPartEvents()
  {
  //returns number of events of particles
    if (!fIsRead) 
-    if ( Read(fParticles,fTracks))
-     {
-       Error("GetNumberOfPartEvents","Error in reading");
-       return 0;
-     }
-   return fParticles->GetNumberOfEvents();
+    {
+      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)?fParticles->GetNumberOfEvents():0;
  }
 
 /********************************************************************/
@@ -125,12 +131,16 @@ Int_t AliHBTReaderTPC::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;
  }
 /********************************************************************/
 
@@ -141,12 +151,8 @@ Int_t AliHBTReaderTPC::Read(AliHBTRun* particles, AliHBTRun *tracks)
  //reurns 0 if everything is OK
  //
   Info("Read","");
-  Int_t i; //iterator and some temprary values
   Int_t Nevents = 0;
   Int_t totalNevents = 0;
-  TFile *aTracksFile;//file with tracks
-  TFile *aClustersFile;//file with clusters
-  TFile *aGAliceFile;//!ile name with galice
 
   if (!particles) //check if an object is instatiated
    {
@@ -177,32 +183,41 @@ Int_t AliHBTReaderTPC::Read(AliHBTRun* particles, AliHBTRun *tracks)
 
   do  //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
    {
-    
-    if( (i=OpenFiles(aTracksFile,aClustersFile,aGAliceFile,currentdir)) )
+    TString filename = GetDirName(currentdir);
+    if (filename.IsNull())
      {
-       Error("Read","Exiting due to problems with opening files. Errorcode %d",i);
+       Error("Read","Can not get directory name");
+       return 4;
+     }
+    filename = filename +"/"+ fFileName;
+    AliRunLoader* rl = AliRunLoader::Open(filename,AliConfig::fgkDefaultEventFolderName);
+    if( rl == 0x0)
+     {
+       Error("Read","Can not open session.");
        currentdir++;
        continue;
      }
-  
     
-    if (gAlice->TreeE())//check if tree E exists
+    rl->LoadHeader();
+    rl->LoadKinematics();
+    AliLoader* tpcl = rl->GetLoader("TPCLoader");
+    
+    if ( tpcl== 0x0)
+     {
+       Error("Read","Exiting due to problems with opening files.");
+       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
       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);
+      rl->LoadgAlice();
+      Info("Read","Setting Magnetic Field: B=%fT",rl->GetAliRun()->Field()->SolenoidField());
+      AliKalmanTrack::SetConvConst(1000/0.299792458/rl->GetAliRun()->Field()->SolenoidField());
+      rl->UnloadgAlice();
      }
     else
      {//if not return an error
@@ -210,20 +225,24 @@ Int_t AliHBTReaderTPC::Read(AliHBTRun* particles, AliHBTRun *tracks)
        currentdir++;
        continue;
      }
-  
-    aGAliceFile->cd();//set cluster file active 
-    AliTPCParam *TPCParam= (AliTPCParam*)aGAliceFile->Get("75x40_100x60_150x60");
-    if (!TPCParam) 
+    
+   rl->CdGAFile();
+   AliTPCParam *TPCParam= (AliTPCParam*)gDirectory->Get("75x40_100x60");
+   
+   if (!TPCParam) 
+    {
+     TPCParam=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
+     if (!TPCParam) 
       { 
-       TPCParam= (AliTPCParam*)aGAliceFile->Get("75x40_100x60");
-       if (!TPCParam) 
-        { 
-          Error("Read","TPC parameters have not been found !\n");
-          currentdir++;
-          continue;
-        }
+        Error("Read","TPC parameters have not been found !\n");
+        delete rl;
+        rl = 0x0;
+        currentdir++;
+        continue;
       }
+    }
 
+    tpcl->LoadTracks();
   
     for(Int_t currentEvent =0; currentEvent<Nevents;currentEvent++)//loop over all events
      {
@@ -231,15 +250,8 @@ Int_t AliHBTReaderTPC::Read(AliHBTRun* particles, AliHBTRun *tracks)
        /**************************************/
         /**************************************/
          /**************************************/ 
-         
-         aTracksFile->cd();//set track file active
-          
-         Char_t  treename[100];
-         sprintf(treename,"TreeT_TPC_%d",currentEvent);//prepare name of the tree
-   
-         TTree *tracktree=0;
-         
-         tracktree=(TTree*)aTracksFile->Get(treename);//get the tree 
+         rl->GetEvent(currentEvent);
+         TTree *tracktree = tpcl->TreeT();//get the tree 
          if (!tracktree) //check if we got the tree
           {//if not return with error
             Error("Read","Can't get a tree with TPC tracks !\n"); 
@@ -258,19 +270,28 @@ Int_t AliHBTReaderTPC::Read(AliHBTRun* particles, AliHBTRun *tracks)
          
          AliTPCtrack *iotrack=0;
          
-//         aClustersFile->cd();//set cluster file active 
+         AliTPCtracker *tracker = new AliTPCtracker(TPCParam,currentEvent,AliConfig::fgkDefaultEventFolderName);//create the tacker for this event
+         if (!tracker) //check if it has created succeffuly
+          {//if not return with error
+            Error("Read","Can't get a tracker !\n"); 
+            continue;
+          }
+         tracker->LoadInnerSectors();
+         tracker->LoadOuterSectors();
    
-         for (i=0; i<NTPCtracks; i++) //loop over all tpc tracks
+         for (Int_t i=0; i<NTPCtracks; i++) //loop over all tpc tracks
           {
             iotrack=new AliTPCtrack;   //create new tracks
             trackbranch->SetAddress(&iotrack); //tell the branch ehere to put track data from tree(file)
             tracktree->GetEvent(i); //stream track i to the iotrack
+            tracker->CookLabel(iotrack,0.1); //calculate (cook) the label of the tpc track
+                                             //which is the label of corresponding simulated particle 
             tarray->AddLast(iotrack); //put the track in the array
           }
          
-         aTracksFile->Delete(treename);//delete tree from memmory (and leave untached on disk)- we do not need it any more
-         aTracksFile->Delete("tracks");//delete branch from memmory
+         delete tracker; //delete tracker
          
+         tracker = 0x0;
          trackbranch = 0x0;
          tracktree = 0x0;
    
@@ -278,27 +299,24 @@ Int_t AliHBTReaderTPC::Read(AliHBTRun* particles, AliHBTRun *tracks)
          Double_t par[5];
          Float_t phi, lam, pt;//angles and transverse momentum
          Int_t label; //label of the current track
-         
-         aGAliceFile->cd();
-         gAlice->GetEvent(currentEvent); 
 
-         gAlice->Particles();
+         rl->Stack()->Particles();
          
-         for (i=0; i<NTPCtracks; i++) //loop over all good tracks
+         for (Int_t i=0; i<NTPCtracks; i++) //loop over all good tracks
           { 
             iotrack = (AliTPCtrack*)tarray->At(i);
             label = iotrack->GetLabel();
-            
+
             if (label < 0) continue;
             
-            
-            TParticle *p = (TParticle*)gAlice->Particle(label);
-            
+            TParticle *p = (TParticle*)rl->Stack()->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,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
@@ -319,7 +337,7 @@ Int_t AliHBTReaderTPC::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(), i, 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
              { 
@@ -329,6 +347,7 @@ Int_t AliHBTReaderTPC::Read(AliHBTRun* particles, AliHBTRun *tracks)
              }
             particles->AddParticle(totalNevents,part);//put track and particle on the run
             tracks->AddParticle(totalNevents,track);
+
           }
          tarray->Clear(); //clear the array
          
@@ -337,107 +356,19 @@ Int_t AliHBTReaderTPC::Read(AliHBTRun* particles, AliHBTRun *tracks)
       /**************************************/  
      totalNevents++;
     }
+  
     //save environment (resouces) --
     //clean your place after the work
-    CloseFiles(aTracksFile,aClustersFile,aGAliceFile); 
+    delete rl;
     currentdir++;
    }while(currentdir < Ndirs);
 
-
   delete tarray;
   fIsRead = kTRUE;
+  
   return 0;
  }
 
-/********************************************************************/
-Int_t AliHBTReaderTPC::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 AliHBTReaderTPC::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;
-   }
-}
-
-/********************************************************************/
 
 /********************************************************************/
 /********************************************************************/