]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EMCAL/AliEMCALGetter.cxx
Logging message with AliLog
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALGetter.cxx
index b785bcd3a85acb50d97e993e923ae461304c25ba..7543cda1bfdad762fcb9504f5ed878c1e63e78ba 100644 (file)
 
 // --- ROOT system ---
 
-#include "TSystem.h"
-#include "TROOT.h"
-
+#include <TFile.h>
+#include <TROOT.h>
+#include <TSystem.h>
+#include <TF1.h>
+#include <TGraph.h>
+//#include <TCanvas.h>
+//#include <TFrame.h>
 
 // --- Standard library ---
 
 // --- AliRoot header files ---
 
-#include "AliEMCALGetter.h"
 #include "AliEMCAL.h"
-#include "AliRunLoader.h"
-#include "AliStack.h"  
-#include "AliHeader.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)
   
@@ -65,10 +71,7 @@ Int_t AliEMCALGetter::fgDebug = 0;
 //  TFile * AliEMCALGetter::fgFile = 0 ; 
 
 //____________________________________________________________________________ 
-AliEMCALGetter::AliEMCALGetter(const char* headerFile, const char* version, Option_t * openingOption)
-{
-  // ctor only called by Instance()
-
+AliEMCALGetter::AliEMCALGetter(const char* headerFile, const char* version, Option_t * openingOption) {
   AliRunLoader* rl = AliRunLoader::GetRunLoader(version) ; 
   if (!rl) {
     rl = AliRunLoader::Open(headerFile, version, openingOption);
@@ -80,46 +83,38 @@ 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
+
+   // initialize data members
   SetDebug(0) ; 
-  //fBTE = 0 ; 
-  fPrimaries = 0 ; 
   fLoadingStatus = "" ; 
+
 }
 
 //____________________________________________________________________________ 
 AliEMCALGetter::~AliEMCALGetter()
 {
-  // dtor
-  delete fgEmcalLoader ;
-  fgEmcalLoader = 0 ;
-  //delete fBTE ; 
-  // fBTE = 0 ; 
-  fPrimaries->Delete() ; 
-  delete fPrimaries ; 
-  fgObjGetter = 0; 
+  if (fgEmcalLoader) {
+    delete fgEmcalLoader;
+    fgEmcalLoader = 0 ;
+  }
+  fgObjGetter = 0;
 }
 
+
 //____________________________________________________________________________ 
 void AliEMCALGetter::Reset()
 {
   // resets things in case the getter is called consecutively with different files
   // the EMCAL Loader is already deleted by the Run Loader
+  fgEmcalLoader = 0;
+  fgObjGetter = 0;
 
-  if (fPrimaries) { 
-    fPrimaries->Delete() ; 
-    delete fPrimaries ;
-  } 
-  fgEmcalLoader = 0; 
-  fgObjGetter = 0; 
 }
 
 //____________________________________________________________________________ 
@@ -226,10 +221,10 @@ void AliEMCALGetter::Event(Int_t event, const char* opt)
 {
   // Reads the content of all Tree's S, D and R
 
-  if ( event >= MaxEvent() ) {
-    Error("Event", "%d not found in TreeE !", event) ; 
-    return ; 
-  }
+//   if ( event >= MaxEvent() ) {
+//     Error("Event", "%d not found in TreeE !", event) ; 
+//     return ; 
+//   }
 
   AliRunLoader * rl = AliRunLoader::GetRunLoader(EmcalLoader()->GetTitle());
   // checks if we are dealing with test-beam data
@@ -271,7 +266,10 @@ void AliEMCALGetter::Event(Int_t event, const char* opt)
 
   if( strstr(opt,"P") )
     ReadTreeP() ;
+
+  if( strstr(opt,"W") )
+    ReadRaw(event) ;
+  
 }
 
 
@@ -303,7 +301,7 @@ 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) ;
   }
@@ -314,7 +312,7 @@ AliEMCALGetter * AliEMCALGetter::Instance(const char* alirunFileName, const char
       TFile * galiceFile = dynamic_cast<TFile *>(gROOT->FindObject(rl->GetFileName()) ) ; 
       
       if ( !galiceFile ) 
-       fgObjGetter = new AliEMCALGetter(alirunFileName, version, openingOption) ;
+       fgObjGetter = new AliEMCALGetter(alirunFileName, version, openingOption);
       
       else {  // the file is already open check the version name
        TString currentVersionName = rl->GetEventFolder()->GetName() ; 
@@ -322,14 +320,14 @@ 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 {    
+         fgObjGetter = new AliEMCALGetter(alirunFileName, version, openingOption) ; 
        }
       }
     }
     else { 
       AliRunLoader * rl = AliRunLoader::GetRunLoader(fgEmcalLoader->GetTitle());
-      if ( strstr(version, AliConfig::fgkDefaultEventFolderName) ) // false in case of merging
+      if ( strstr(version, AliConfig::GetDefaultEventFolderName()) ) // false in case of merging
        delete rl ; 
       fgObjGetter = new AliEMCALGetter(alirunFileName, version, openingOption) ;      
     }
@@ -339,7 +337,7 @@ AliEMCALGetter * AliEMCALGetter::Instance(const char* alirunFileName, const char
   else 
     if (fgDebug)
       Print() ;
-  
+
   return fgObjGetter ;
 }
 
@@ -371,11 +369,27 @@ 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  
 {
   // returns the EMCAL object 
-  AliEMCAL * emcal = dynamic_cast<AliEMCAL*>(EmcalLoader()->GetModulesFolder()->FindObject("EMCAL")) ;  
+  AliEMCALLoader *    loader = 0;
+  static AliEMCALLoader * oldloader = 0;
+  static AliEMCAL * emcal = 0;
+
+  loader = EmcalLoader();
+
+  if (loader != oldloader ) {
+    emcal = dynamic_cast<AliEMCAL*>(loader->GetModulesFolder()->FindObject("EMCAL")) ;  
+    oldloader = loader;
+  }
   if (!emcal) 
     if (fgDebug)
       Warning("EMCAL", "EMCAL module not found in module folders: %s", EmcalLoader()->GetModulesFolder()->GetName() ) ; 
@@ -408,18 +422,6 @@ AliEMCALGeometry * AliEMCALGetter::EMCALGeometry() const
   return rv ; 
 } 
 
-//____________________________________________________________________________ 
-TClonesArray * AliEMCALGetter::Primaries()  
-{
-  // creates the Primaries container if needed
-  if ( !fPrimaries ) {
-    if (fgDebug) 
-      Info("Primaries", "Creating a new TClonesArray for primaries") ; 
-    fPrimaries = new TClonesArray("TParticle", 1000) ;
-  } 
-  return fPrimaries ; 
-}
-
 //____________________________________________________________________________ 
 void  AliEMCALGetter::Print() 
 {
@@ -440,45 +442,202 @@ void AliEMCALGetter::ReadPrimaries()
   if ( ! rl->TreeK() )  // load treeK the first time
     rl->LoadKinematics() ;
   
-  fNPrimaries = (rl->GetHeader())->GetNtrack(); 
   if (fgDebug) 
-    Info("ReadTreeK", "Found %d particles in event # %d", fNPrimaries, EventNumber() ) ; 
+    Info("ReadTreeK", "Found %d particles in event # %d", NPrimaries(), EventNumber() ) ; 
+}
+
+//____________________________________________________________________________ 
+void AliEMCALGetter::FitRaw(Bool_t lowGainFlag, TGraph * gLowGain, TGraph * gHighGain, TF1* signalF, Int_t & amp, Double_t & time)
+{
+  // Fits the raw signal time distribution 
+
+  const Int_t kNoiseThreshold = 0 ;
+  Double_t timezero1 = 0., timezero2 = 0., timemax = 0. ;
+  Double_t signal = 0., signalmax = 0. ;       
+  Double_t energy = time = 0. ; 
 
-  // first time creates the container
-  if ( Primaries() ) 
-    fPrimaries->Clear() ; 
+  if (lowGainFlag) {
+    timezero1 = timezero2 = signalmax = timemax = 0. ;
+    signalF->FixParameter(0, EMCAL()->GetRawFormatLowCharge()) ; 
+    signalF->FixParameter(1, EMCAL()->GetRawFormatLowGain()) ; 
+    Int_t index ; 
+    for (index = 0; index < EMCAL()->GetRawFormatTimeBins(); index++) {
+      gLowGain->GetPoint(index, time, signal) ; 
+      if (signal > kNoiseThreshold && timezero1 == 0.) 
+       timezero1 = time ;
+      if (signal <= kNoiseThreshold && timezero1 > 0. && timezero2 == 0.)
+       timezero2 = time ; 
+      if (signal > signalmax) {
+       signalmax = signal ; 
+       timemax   = time ; 
+      }
+    }
+    signalmax /= EMCAL()->RawResponseFunctionMax(EMCAL()->GetRawFormatLowCharge(), 
+                                               EMCAL()->GetRawFormatLowGain()) ;
+    if ( timezero1 + EMCAL()->GetRawFormatTimePeak() < EMCAL()->GetRawFormatTimeMax() * 0.4 ) { // else its noise 
+      signalF->SetParameter(2, signalmax) ; 
+      signalF->SetParameter(3, timezero1) ;                
+      gLowGain->Fit(signalF, "QRON", "", 0., timezero2); //, "QRON") ; 
+      energy = signalF->GetParameter(2) ; 
+      time   = signalF->GetMaximumX() - EMCAL()->GetRawFormatTimePeak() - EMCAL()->GetRawFormatTimeTrigger() ;
+    }
+  } else {
+    timezero1 = timezero2 = signalmax = timemax = 0. ;
+    signalF->FixParameter(0, EMCAL()->GetRawFormatHighCharge()) ; 
+    signalF->FixParameter(1, EMCAL()->GetRawFormatHighGain()) ; 
+    Int_t index ; 
+    for (index = 0; index < EMCAL()->GetRawFormatTimeBins(); index++) {
+      gHighGain->GetPoint(index, time, signal) ;               
+      if (signal > kNoiseThreshold && timezero1 == 0.) 
+       timezero1 = time ;
+      if (signal <= kNoiseThreshold && timezero1 > 0. && timezero2 == 0.)
+       timezero2 = time ; 
+      if (signal > signalmax) {
+       signalmax = signal ;   
+       timemax   = time ; 
+      }
+    }
+    signalmax /= EMCAL()->RawResponseFunctionMax(EMCAL()->GetRawFormatHighCharge(), 
+                                               EMCAL()->GetRawFormatHighGain()) ;;
+    if ( timezero1 + EMCAL()->GetRawFormatTimePeak() < EMCAL()->GetRawFormatTimeMax() * 0.4 ) { // else its noise  
+      signalF->SetParameter(2, signalmax) ; 
+      signalF->SetParameter(3, timezero1) ;               
+      gHighGain->Fit(signalF, "QRON", "", 0., timezero2) ; 
+      energy = signalF->GetParameter(2) ; 
+      time   = signalF->GetMaximumX() - EMCAL()->GetRawFormatTimePeak() - EMCAL()->GetRawFormatTimeTrigger() ;
+    }
+  }
   
-  Int_t index = 0 ; 
-  for (index = 0 ; index < fNPrimaries; index++) { 
-    new ((*fPrimaries)[index]) TParticle(*(Primary(index)));
+  if (time == 0. && energy == 0.) 
+    amp = 0 ; 
+  else {
+  AliEMCALDigitizer * digitizer = Digitizer() ; 
+  amp = static_cast<Int_t>( (energy - digitizer->GetECApedestal()) / digitizer->GetECAchannel() + 0.5 ) ; 
   }
-}
+  // dessin
+//   TCanvas * c1 = new TCanvas("c1","A Simple Graph Example",200,10,700,500);
+//   c1->SetFillColor(42);
+//   c1->SetGrid();
+//   gLowGain->SetLineColor(2);
+//   gLowGain->SetLineWidth(4);
+//   gLowGain->SetMarkerColor(4);
+//   gLowGain->SetMarkerStyle(21);
+//   gLowGain->SetTitle("Lowgain");
+//   gLowGain->GetXaxis()->SetTitle("X title");
+//   gLowGain->GetYaxis()->SetTitle("Y title");
+//   gLowGain->Draw("ACP");
+  
+//   c1->Update();
+//   c1->GetFrame()->SetFillColor(21);
+//   c1->GetFrame()->SetBorderSize(12);
+//   c1->Modified();
+  
+//   TCanvas * c2 = new TCanvas("c2","A Simple Graph Example",200,10,700,500);
+//   c2->SetFillColor(42);
+//   c2->SetGrid();
+//   gHighGain->SetLineColor(2);
+//   gHighGain->SetLineWidth(4);
+//   gHighGain->SetMarkerColor(4);
+//   gHighGain->SetMarkerStyle(21);
+//   gHighGain->SetTitle("Highgain");
+//   gHighGain->GetXaxis()->SetTitle("X title");
+//   gHighGain->GetYaxis()->SetTitle("Y title");
+//     gHighGain->Draw("ACP");
+  
+//   c2->Update();
+//   c2->GetFrame()->SetFillColor(21);
+//   c2->GetFrame()->SetBorderSize(12);
+//   c2->Modified();
+}
+
+//____________________________________________________________________________ 
+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);
+  
+  Bool_t first = kTRUE ;
+  TF1 * signalF = new TF1("signal", AliEMCAL::RawResponseFunction, 0, EMCAL()->GetRawFormatTimeMax(), 4);
+  signalF->SetParNames("Charge", "Gain", "Amplitude", "TimeZero") ; 
+  
+  Int_t id = -1;
+  Bool_t lowGainFlag = kFALSE ; 
+
+  TClonesArray * digits = Digits() ;
+  digits->Clear() ; 
+  Int_t idigit = 0 ; 
+  Int_t amp = 0 ; 
+  Double_t time = 0. ; 
+
+  TGraph * gLowGain = new TGraph(EMCAL()->GetRawFormatTimeBins()) ; 
+  TGraph * gHighGain= new TGraph(EMCAL()->GetRawFormatTimeBins()) ;  
+
+  while ( in.Next() ) { // EMCAL entries loop 
+    if ( in.IsNewId() ) {
+      if (!first) {
+       FitRaw(lowGainFlag, gLowGain, gHighGain, signalF, amp, time) ; 
+       if (amp > 0) {
+         new((*digits)[idigit]) AliEMCALDigit( -1, -1, id, amp, time) ;        
+         idigit++ ; 
+       }
+       Int_t index ; 
+       for (index = 0; index < EMCAL()->GetRawFormatTimeBins(); index++) {
+         gLowGain->SetPoint(index, index * EMCAL()->GetRawFormatTimeMax() / EMCAL()->GetRawFormatTimeBins(), 0) ;  
+         gHighGain->SetPoint(index, index * EMCAL()->GetRawFormatTimeMax() / EMCAL()->GetRawFormatTimeBins(), 0) ; 
+       } 
+      }  
+      first = kFALSE ; 
+      id = in.GetId() ; 
+      if (in.GetModule() == EMCAL()->GetRawFormatLowGainOffset() ) 
+       lowGainFlag = kTRUE ; 
+      else 
+       lowGainFlag = kFALSE ; 
+    }
+    if (lowGainFlag)
+      gLowGain->SetPoint(in.GetTime(), 
+                        in.GetTime()* EMCAL()->GetRawFormatTimeMax() / EMCAL()->GetRawFormatTimeBins(), 
+                        in.GetSignal()) ;
+    else 
+      gHighGain->SetPoint(in.GetTime(), 
+                         in.GetTime() * EMCAL()->GetRawFormatTimeMax() / EMCAL()->GetRawFormatTimeBins(), 
+                         in.GetSignal() ) ;
+    
+  } // EMCAL entries loop
+  digits->Sort() ; 
 
-//____________________________________________________________________________ 
+  delete signalF ; 
+  delete gLowGain;
+  delete gHighGain ; 
+    
+  return Digits()->GetEntriesFast() ; 
+}
+  
+  //____________________________________________________________________________ 
 Int_t AliEMCALGetter::ReadTreeD()
 {
   // Read the Digits
   
-  
-  // gets TreeD from the root file (EMCAL.SDigits.root)
-  if ( !IsLoaded("D") ) {
-    EmcalLoader()->LoadDigits("UPDATE") ;
-    EmcalLoader()->LoadDigitizer("UPDATE") ;
-    SetLoaded("D") ; 
-  } 
-  return Digits()->GetEntries() ; 
+   EmcalLoader()->CleanDigits() ; 
+   EmcalLoader()->LoadDigits("UPDATE") ;
+   EmcalLoader()->LoadDigitizer("UPDATE") ;
+   return Digits()->GetEntries() ; 
 }
 
 //____________________________________________________________________________ 
 Int_t AliEMCALGetter::ReadTreeH()
 {
   // Read the Hits
-    
+  EmcalLoader()->CleanHits() ; 
   // gets TreeH from the root file (EMCAL.Hit.root)
-  if ( !IsLoaded("H") ) {
-    EmcalLoader()->LoadHits("READ") ;
-    SetLoaded("H") ; 
-  }  
+  //if ( !IsLoaded("H") ) {
+    EmcalLoader()->LoadHits("UPDATE") ;
+    //SetLoaded("H") ; 
+    //}  
   return Hits()->GetEntries() ; 
 }
 
@@ -486,14 +645,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() ; 
 }
@@ -503,28 +662,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() ; 
 }
@@ -533,13 +692,13 @@ Int_t AliEMCALGetter::ReadTreeS()
 {
   // Read the SDigits
   
-  
+  //  EmcalLoader()->CleanSDigits() ; 
   // gets TreeS from the root file (EMCAL.SDigits.root)
-   if ( !IsLoaded("S") ) {
+  // if ( !IsLoaded("S") ) {
     EmcalLoader()->LoadSDigits("READ") ;
     EmcalLoader()->LoadSDigitizer("READ") ;
-    SetLoaded("S") ; 
-  }
+    //  SetLoaded("S") ; 
+    //}
 
   return SDigits()->GetEntries() ; 
 }
@@ -710,7 +869,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()) ;
@@ -722,7 +881,7 @@ 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()) ;