]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EMCAL/AliEMCALGetter.cxx
Corrected posting and cleaning of digitizers (T.Kuhr)
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALGetter.cxx
index 18cc5775d975c4ce6ede428862e945d6b4d50158..fb48153b7c758f1a17f8e7fe20fab9a02af71389 100644 (file)
 
 // --- ROOT system ---
 
-#include "TSystem.h"
-#include "TFile.h"
-#include "TROOT.h"
-
+#include <TFile.h>
+#include <TROOT.h>
+#include <TSystem.h>
+#include <TH1D.h>
+#include <TF1.h>
 
 // --- Standard library ---
 
 // --- AliRoot header files ---
 
-#include "AliEMCALGetter.h"
 #include "AliEMCAL.h"
-#include "AliRunLoader.h"
-#include "AliStack.h"  
+#include "AliEMCALGetter.h"
 #include "AliEMCALLoader.h"
+#include "AliHeader.h"  
 #include "AliMC.h"
+#include "AliRunLoader.h"
+#include "AliStack.h"  
+#include "AliEMCALRawStream.h"
+#include "AliRawReaderFile.h"
 
 ClassImp(AliEMCALGetter)
   
 AliEMCALGetter * AliEMCALGetter::fgObjGetter = 0 ; 
 AliEMCALLoader * AliEMCALGetter::fgEmcalLoader = 0;
 Int_t AliEMCALGetter::fgDebug = 0;
+TString AliEMCALGetter::fVersion = "";
 
 //  TFile * AliEMCALGetter::fgFile = 0 ; 
 
 //____________________________________________________________________________ 
-AliEMCALGetter::AliEMCALGetter(const char* headerFile, const char* version, Option_t * openingOption)
+AliEMCALGetter::AliEMCALGetter(const char* headerFile, const char* version, Option_t * openingOption) 
 {
   // ctor only called by Instance()
 
+  // initialize data members
+  SetDebug(0) ; 
+  //fBTE = 0 ; 
+  
+  fLoadingStatus = "" ; 
+
+  fgObjGetter=this;
+  
+  OpenFile(headerFile,version,openingOption);
+}
+
+
+//____________________________________________________________________________ 
+AliEMCALGetter::~AliEMCALGetter()
+{
+  //PH  AliRunLoader * rl = AliRunLoader::GetRunLoader(fgEmcalLoader->GetTitle());
+  //PH  delete rl;
+  fgEmcalLoader = 0 ;
+  fgObjGetter = 0; 
+  fVersion = "";
+}
+
+//____________________________________________________________________________ 
+void AliEMCALGetter::OpenFile(const char* headerFile, const char* version, Option_t * openingOption) {
+  fVersion = version;
   AliRunLoader* rl = AliRunLoader::GetRunLoader(version) ; 
   if (!rl) {
     rl = AliRunLoader::Open(headerFile, version, openingOption);
@@ -80,31 +110,20 @@ AliEMCALGetter::AliEMCALGetter(const char* headerFile, const char* version, Opti
       rl->LoadgAlice();
       gAlice = rl->GetAliRun(); // should be removed
     }
-  }
+  } 
   fgEmcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetLoader("EMCALLoader"));
   if ( !fgEmcalLoader ) 
     Error("AliEMCALGetter", "Could not find EMCALLoader") ; 
   else 
     fgEmcalLoader->SetTitle(version);
-  
-  
-  // initialize data members
-  SetDebug(0) ; 
-  //fBTE = 0 ; 
-  fPrimaries = 0 ; 
-  fLoadingStatus = "" ; 
 }
 
 //____________________________________________________________________________ 
-AliEMCALGetter::~AliEMCALGetter()
+void AliEMCALGetter::Reset()
 {
-  // dtor
-  delete fgEmcalLoader ;
-  fgEmcalLoader = 0 ;
-  //delete fBTE ; 
-  // fBTE = 0 ; 
-  fPrimaries->Delete() ; 
-  delete fPrimaries ; 
+  // resets things in case the getter is called consecutively with different files
+  // the EMCAL Loader is already deleted by the Run Loader
+
 }
 
 //____________________________________________________________________________ 
@@ -207,7 +226,7 @@ TClonesArray * AliEMCALGetter::RecParticles() const
   return rv ; 
 }
 //____________________________________________________________________________ 
-void AliEMCALGetter::Event(const Int_t event, const char* opt) 
+void AliEMCALGetter::Event(Int_t event, const char* opt) 
 {
   // Reads the content of all Tree's S, D and R
 
@@ -217,7 +236,6 @@ void AliEMCALGetter::Event(const Int_t event, const char* opt)
   }
 
   AliRunLoader * rl = AliRunLoader::GetRunLoader(EmcalLoader()->GetTitle());
-
   // checks if we are dealing with test-beam data
 //   TBranch * btb = rl->TreeE()->GetBranch("AliEMCALBeamTestEvent") ;
 //   if(btb){
@@ -257,7 +275,10 @@ void AliEMCALGetter::Event(const Int_t event, const char* opt)
 
   if( strstr(opt,"P") )
     ReadTreeP() ;
+
+  if( strstr(opt,"W") )
+    ReadRaw(event) ;
+  
 }
 
 
@@ -289,18 +310,19 @@ AliEMCALGetter * AliEMCALGetter::Instance(const char* alirunFileName, const char
 {
   // Creates and returns the pointer of the unique instance
   // Must be called only when the environment has changed
-  
+
   if(!fgObjGetter){ // first time the getter is called 
     fgObjGetter = new AliEMCALGetter(alirunFileName, version, openingOption) ;
   }
   else { // the getter has been called previously
-    AliRunLoader * rl = AliRunLoader::GetRunLoader(fgEmcalLoader->GetTitle());
-    if ( rl->GetFileName() == alirunFileName ) {// the alirunFile has the same name
+    AliRunLoader * rl = AliRunLoader::GetRunLoader(fVersion);
+    if (rl == 0)  fgObjGetter->OpenFile(alirunFileName, version, openingOption) ; 
+    else if ( rl->GetFileName() == alirunFileName ) {// the alirunFile has the same name
       // check if the file is already open
       TFile * galiceFile = dynamic_cast<TFile *>(gROOT->FindObject(rl->GetFileName()) ) ; 
       
       if ( !galiceFile ) 
-       fgObjGetter = new AliEMCALGetter(alirunFileName, version, openingOption) ;
+       fgObjGetter->OpenFile(alirunFileName, version, openingOption);
       
       else {  // the file is already open check the version name
        TString currentVersionName = rl->GetEventFolder()->GetName() ; 
@@ -308,20 +330,24 @@ AliEMCALGetter * AliEMCALGetter::Instance(const char* alirunFileName, const char
        if (currentVersionName == newVersionName) 
          if(fgDebug)
            ::Warning( "Instance", "Files with version %s already open", currentVersionName.Data() ) ;  
-       else {
-         fgObjGetter = new AliEMCALGetter(alirunFileName, version, openingOption) ;      
+       else {    
+         fgEmcalLoader->SetTitle(version); fVersion = version;
        }
       }
     }
-    else 
-      fgObjGetter = new AliEMCALGetter(alirunFileName, version, openingOption) ;      
+    else { 
+      AliRunLoader * rl = AliRunLoader::GetRunLoader(fVersion);
+      if ( strstr(version, AliConfig::GetDefaultEventFolderName()) ) // false in case of merging
+       delete rl ; 
+      fgObjGetter->OpenFile(alirunFileName, version, openingOption) ;      
+    }
   }
   if (!fgObjGetter) 
-    ::Error("Instance", "Failed to create the EMCAL Getter object") ;
+    ::Error("AliEMCALGetter::Instance", "Failed to create the EMCAL Getter object") ;
   else 
     if (fgDebug)
       Print() ;
-  
+
   return fgObjGetter ;
 }
 
@@ -330,8 +356,8 @@ AliEMCALGetter *  AliEMCALGetter::Instance()
 {
   // Returns the pointer of the unique instance already defined
   
-  if(!fgObjGetter)
-     ::Error("Instance", "Getter not initialized") ;
+  if(!fgObjGetter && fgDebug)
+     ::Warning("AliEMCALGetter::Instance", "Getter not initialized") ;
 
    return fgObjGetter ;
            
@@ -353,6 +379,13 @@ TParticle * AliEMCALGetter::Primary(Int_t index) const
   return rl->Stack()->Particle(index) ; 
 } 
 
+//____________________________________________________________________________ 
+Int_t AliEMCALGetter::NPrimaries() const
+{
+  AliRunLoader * rl = AliRunLoader::GetRunLoader(EmcalLoader()->GetTitle());
+  return (rl->GetHeader())->GetNtrack(); 
+} 
+
 //____________________________________________________________________________ 
 AliEMCAL * AliEMCALGetter:: EMCAL() const  
 {
@@ -390,25 +423,13 @@ AliEMCALGeometry * AliEMCALGetter::EMCALGeometry() const
   return rv ; 
 } 
 
-//____________________________________________________________________________ 
-TClonesArray * AliEMCALGetter::Primaries()  
-{
-  // creates the Primaries container if needed
-  if ( !fPrimaries ) {
-    if (fgDebug) 
-      printf("Primaries: Creating a new TClonesArray for primaries") ; 
-    fPrimaries = new TClonesArray("TParticle", 1000) ;
-  } 
-  return fPrimaries ; 
-}
-
 //____________________________________________________________________________ 
 void  AliEMCALGetter::Print() 
 {
   // Print usefull information about the getter
     
   AliRunLoader * rl = AliRunLoader::GetRunLoader(fgEmcalLoader->GetTitle());
-  printf("Print: gAlice file is %s -- version name is %s", (rl->GetFileName()).Data(), rl->GetEventFolder()->GetName() ) ; 
+  ::Info("Print", "gAlice file is %s -- version name is %s", (rl->GetFileName()).Data(), rl->GetEventFolder()->GetName() ) ; 
 }
 
 //____________________________________________________________________________ 
@@ -422,34 +443,101 @@ void AliEMCALGetter::ReadPrimaries()
   if ( ! rl->TreeK() )  // load treeK the first time
     rl->LoadKinematics() ;
   
-  fNPrimaries = rl->Stack()->GetNtrack() ; 
-
   if (fgDebug) 
-    printf("ReadTreeK: Found %d particles in event # %d", fNPrimaries, EventNumber() ) ; 
-
-
-  // first time creates the container
-  if ( Primaries() ) 
-    fPrimaries->Clear() ; 
-  
-  Int_t index = 0 ; 
-  for (index = 0 ; index < fNPrimaries; index++) { 
-    new ((*fPrimaries)[index]) TParticle(*(Primary(index)));
-  }
+    Info("ReadTreeK", "Found %d particles in event # %d", NPrimaries(), EventNumber() ) ; 
 }
 
 //____________________________________________________________________________ 
+Int_t AliEMCALGetter::ReadRaw(Int_t event)
+{
+  // reads the raw format data, converts it into digits format and store digits in Digits()
+  // container.
+  
+  AliRawReaderFile rawReader(event) ; 
+  AliEMCALRawStream in(&rawReader);
+  
+  const Int_t kHighGainIdOffset = EMCALGeometry()->GetNTowers()
+    * EMCALGeometry()->GetNPhi()
+    * EMCALGeometry()->GetNZ() 
+    * 2 ;
+  
+  Bool_t first = kTRUE ;
+  
+  TH1D hLowGain("hLowGain", "Low Gain", 1000, 0., EMCAL()->GetRawFormatTimeMax()) ; 
+  TH1D hHighGain("hHighGain", "High Gain", 1000, 0., EMCAL()->GetRawFormatTimeMax()) ; 
+  
+  // fit half the gaussian decay rather than AliEMCAL::RawResponseFunction because thiswould give a floating point
+  // exception during error calculation ... To solve... 
+  TF1 * gauss = new TF1("gauss", "gaus", 
+                       EMCAL()->GetRawFormatTimePeak(), 
+                       EMCAL()->GetRawFormatTimeMax() ) ;   
+  
+  Int_t id = -1;
+  Bool_t hgflag = kFALSE ; 
+  
+  TClonesArray * digits = Digits() ;
+  digits->Clear() ; 
+  Int_t idigit = 0 ; 
+  
+  while ( in.Next() ) { // EMCAL entries loop 
+    
+    if ( in.IsNewId() ) {
+      if (!first) {
+       hLowGain.Fit(gauss, "QRON") ; 
+       Int_t ampL =  static_cast<Int_t>(gauss->Eval(gauss->GetParameter(2)) + 0.5) ; 
+       Double_t timeL  = EMCAL()->GetRawFormatTimePeak() - gauss->GetParameter(2) ;
+       if (timeL < 0 ) // happens with noise 
+         timeL =  EMCAL()->GetRawFormatTimePeak() ;  
+       if (hgflag) {
+         hHighGain.Fit(gauss, "QRON") ; 
+         Int_t ampH =  static_cast<Int_t>(gauss->Eval(gauss->GetParameter(2)) + 0.5) ; 
+         Double_t timeH  = EMCAL()->GetRawFormatTimePeak() - gauss->GetParameter(2) ; 
+         if (timeH < 0 ) // happens with noise 
+           timeH =  EMCAL()->GetRawFormatTimePeak() ;  
+         new((*digits)[idigit]) AliEMCALDigit( -1, -1, id+kHighGainIdOffset, ampH, timeH) ;
+         idigit++ ; 
+       }
+       else {
+         new((*digits)[idigit]) AliEMCALDigit( -1, -1, id, ampL, timeL) ;
+         idigit++ ; 
+       }
+      }
+      first = kFALSE ; 
+      hLowGain.Reset() ; 
+      hHighGain.Reset() ; 
+      id = in.GetId() ; 
+      if (id > 9999 ) { // fixme 
+       hgflag = kTRUE ; 
+      } else 
+       hgflag = kFALSE ; 
+    }
+    if (hgflag)
+      hHighGain.Fill(
+                    in.GetTime() * EMCAL()->GetRawFormatTimeMax() / EMCAL()->GetRawFormatTimeBins(), 
+                    in. GetSignal() * EMCAL()->GetRawFormatHighGainFactor() ) ;
+    else 
+      hLowGain.Fill(
+                   in.GetTime() * EMCAL()->GetRawFormatTimeMax() / EMCAL()->GetRawFormatTimeBins(), 
+                   in. GetSignal() ) ;
+  } // EMCAL entries loop
+  
+  delete gauss ; 
+  
+  return Digits()->GetEntriesFast() ; 
+}
+  
+  //____________________________________________________________________________ 
 Int_t AliEMCALGetter::ReadTreeD()
 {
   // Read the Digits
   
-  
-  // gets TreeD from the root file (EMCAL.SDigits.root)
-  if ( !IsLoaded("D") ) {
+   EmcalLoader()->CleanDigits() ; 
+  // gets TreeD from the root file (EMCAL.Digits.root)
+  //if ( !IsLoaded("D") ) {
     EmcalLoader()->LoadDigits("UPDATE") ;
     EmcalLoader()->LoadDigitizer("UPDATE") ;
-    SetLoaded("D") ; 
-  } 
+    //  SetLoaded("D") ; 
+    //
   return Digits()->GetEntries() ; 
 }
 
@@ -457,12 +545,12 @@ Int_t AliEMCALGetter::ReadTreeD()
 Int_t AliEMCALGetter::ReadTreeH()
 {
   // Read the Hits
-    
+  EmcalLoader()->CleanHits() ; 
   // gets TreeH from the root file (EMCAL.Hit.root)
-  if ( !IsLoaded("H") ) {
-    EmcalLoader()->LoadHits("UPDATE") ;
+  //if ( !IsLoaded("H") ) {
+    EmcalLoader()->LoadHits("READ") ;
     SetLoaded("H") ; 
-  }  
+    //}  
   return Hits()->GetEntries() ; 
 }
 
@@ -470,14 +558,14 @@ Int_t AliEMCALGetter::ReadTreeH()
 Int_t AliEMCALGetter::ReadTreeR()
 {
   // Read the RecPoints
-  
-  
+
+   EmcalLoader()->CleanRecPoints() ; 
   // gets TreeR from the root file (EMCAL.RecPoints.root)
-  if ( !IsLoaded("R") ) {
+  //if ( !IsLoaded("R") ) {
     EmcalLoader()->LoadRecPoints("UPDATE") ;
     EmcalLoader()->LoadClusterizer("UPDATE") ;
-    SetLoaded("R") ; 
-  }
+    //  SetLoaded("R") ; 
+    //}
 
   return ECARecPoints()->GetEntries() ; 
 }
@@ -487,28 +575,28 @@ Int_t AliEMCALGetter::ReadTreeT()
 {
   // Read the TrackSegments
   
-  
+  EmcalLoader()->CleanTracks() ; 
   // gets TreeT from the root file (EMCAL.TrackSegments.root)
-  if ( !IsLoaded("T") ) {
+  //if ( !IsLoaded("T") ) {
     EmcalLoader()->LoadTracks("UPDATE") ;
     EmcalLoader()->LoadTrackSegmentMaker("UPDATE") ;
-    SetLoaded("T") ; 
-  }
+    //  SetLoaded("T") ; 
+    //}
 
   return TrackSegments()->GetEntries() ; 
 }
 //____________________________________________________________________________ 
 Int_t AliEMCALGetter::ReadTreeP()
 {
-  // Read the TrackSegments
-  
+  // Read the RecParticles
   
-  // gets TreeT from the root file (EMCAL.TrackSegments.root)
-  if ( !IsLoaded("P") ) {
+  EmcalLoader()->CleanRecParticles() ; 
+  // gets TreeP from the root file (EMCAL.RecParticles.root)
+  //  if ( !IsLoaded("P") ) {
     EmcalLoader()->LoadRecParticles("UPDATE") ;
     EmcalLoader()->LoadPID("UPDATE") ;
-    SetLoaded("P") ; 
-  }
+    //  SetLoaded("P") ; 
+    // }
 
   return RecParticles()->GetEntries() ; 
 }
@@ -517,13 +605,13 @@ Int_t AliEMCALGetter::ReadTreeS()
 {
   // Read the SDigits
   
-  
+  //  EmcalLoader()->CleanSDigits() ; 
   // gets TreeS from the root file (EMCAL.SDigits.root)
-  if ( !IsLoaded("S") ) {
-    EmcalLoader()->LoadSDigits("UPDATE") ;
-    EmcalLoader()->LoadSDigitizer("UPDATE") ;
-    SetLoaded("S") ; 
-  }
+  // if ( !IsLoaded("S") ) {
+    EmcalLoader()->LoadSDigits("READ") ;
+    EmcalLoader()->LoadSDigitizer("READ") ;
+    //  SetLoaded("S") ; 
+    //}
 
   return SDigits()->GetEntries() ; 
 }
@@ -546,7 +634,7 @@ TClonesArray * AliEMCALGetter::SDigits() const
 //____________________________________________________________________________ 
 AliEMCALSDigitizer * AliEMCALGetter::SDigitizer() 
 { 
-  // return pointer to SDigitizer Tree
+  // Return pointer to SDigitizer task
   AliEMCALSDigitizer * rv ; 
   rv =  dynamic_cast<AliEMCALSDigitizer *>(EmcalLoader()->SDigitizer()) ;
   if (!rv) {
@@ -557,7 +645,7 @@ AliEMCALSDigitizer * AliEMCALGetter::SDigitizer()
 }
 
 //____________________________________________________________________________ 
-TParticle * AliEMCALGetter::Secondary(const TParticle* p, const Int_t index) const
+TParticle * AliEMCALGetter::Secondary(const TParticle* p, Int_t index) const
 {
   // Return first (index=1) or second (index=2) secondary particle of primary particle p 
 
@@ -576,7 +664,7 @@ TParticle * AliEMCALGetter::Secondary(const TParticle* p, const Int_t index) con
 }
 
 //____________________________________________________________________________ 
-void AliEMCALGetter::Track(const Int_t itrack) 
+void AliEMCALGetter::Track(Int_t itrack) 
 {
   // Read the first entry of EMCAL branch in hit tree gAlice->TreeH()
  
@@ -694,7 +782,7 @@ Bool_t AliEMCALGetter::VersionExists(TString & opt) const
   if ( opt == "sdigits") {
     // add the version name to the root file name
     TString fileName( EmcalLoader()->GetSDigitsFileName() ) ; 
-    if (version != AliConfig::fgkDefaultEventFolderName) // only if not the default folder name 
+    if (version != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name 
       fileName = fileName.ReplaceAll(".root", "") + "_" + version + ".root" ;
     if ( !(gSystem->AccessPathName(fileName)) ) { 
       Warning("VersionExists", "The file %s already exists", fileName.Data()) ;
@@ -706,13 +794,12 @@ Bool_t AliEMCALGetter::VersionExists(TString & opt) const
   if ( opt == "digits") {
     // add the version name to the root file name
     TString fileName( EmcalLoader()->GetDigitsFileName() ) ; 
-    if (version != AliConfig::fgkDefaultEventFolderName) // only if not the default folder name 
+    if (version != AliConfig::GetDefaultEventFolderName()) // only if not the default folder name 
       fileName = fileName.ReplaceAll(".root", "") + "_" + version + ".root" ;
     if ( !(gSystem->AccessPathName(fileName)) ) {
       Warning("VersionExists", "The file %s already exists", fileName.Data()) ;  
       rv = kTRUE ; 
     }
-    EmcalLoader()->SetDigitsFileName(fileName) ;
   }
 
   return rv ;