]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - FMD/AliFMDInput.cxx
Fixes for reading zero-suppressed data. These should be propagated to
[u/mrichter/AliRoot.git] / FMD / AliFMDInput.cxx
index 1a77c393a4e22c3b6b258e27e1409e18749a3d2c..4496982720aa754bdee0b3cedfe832d37bcd167f 100644 (file)
 // Latest changes by Christian Holm Christensen
 //
 #include "AliFMDInput.h"       // ALIFMDHIT_H
-#include "AliFMDDebug.h"               // ALIFMDDEBUG_H ALILOG_H
+#include "AliFMDDebug.h"       // ALIFMDDEBUG_H ALILOG_H
 #include "AliLoader.h"          // ALILOADER_H
 #include "AliRunLoader.h"       // ALIRUNLOADER_H
 #include "AliRun.h"             // ALIRUN_H
 #include "AliStack.h"           // ALISTACK_H
 #include "AliRawReaderFile.h"   // ALIRAWREADERFILE_H
+#include "AliRawReaderRoot.h"   // ALIRAWREADERROOT_H
+#include "AliRawReaderDate.h"   // ALIRAWREADERDATE_H
 #include "AliFMD.h"             // ALIFMD_H
 #include "AliFMDHit.h"         // ALIFMDHIT_H
 #include "AliFMDDigit.h"       // ALIFMDDigit_H
@@ -58,6 +60,8 @@
 #include <Riostream.h>         // ROOT_Riostream
 #include <TFile.h>              // ROOT_TFile
 #include <TStreamerInfo.h>
+#include <TArrayF.h>
+
 //____________________________________________________________________
 ClassImp(AliFMDInput)
 #if 0
@@ -67,12 +71,14 @@ ClassImp(AliFMDInput)
 
 //____________________________________________________________________
 AliFMDInput::AliFMDInput()
-  : fGAliceFile(""), 
+  : TNamed("AliFMDInput", "Input handler for various FMD data"), 
+    fGAliceFile(""), 
     fLoader(0),
     fRun(0), 
     fStack(0),
     fFMDLoader(0), 
     fReader(0),
+    fFMDReader(0),
     fFMD(0),
     fESD(0),
     fESDEvent(0),
@@ -89,9 +95,12 @@ AliFMDInput::AliFMDInput()
     fArrayS(0), 
     fArrayR(0), 
     fArrayA(0), 
+    fHeader(0),
     fGeoManager(0),
     fTreeMask(0), 
-    fIsInit(kFALSE)
+    fRawFile(""),
+    fIsInit(kFALSE),
+    fEventCount(0)
 {
 
   // Constructor of an FMD input object.  Specify what data to read in
@@ -104,12 +113,14 @@ AliFMDInput::AliFMDInput()
 
 //____________________________________________________________________
 AliFMDInput::AliFMDInput(const char* gAliceFile)
-  : fGAliceFile(gAliceFile),
+  : TNamed("AliFMDInput", "Input handler for various FMD data"), 
+    fGAliceFile(gAliceFile),
     fLoader(0),
     fRun(0), 
     fStack(0),
     fFMDLoader(0), 
     fReader(0),
+    fFMDReader(0),
     fFMD(0),
     fESD(0),
     fESDEvent(0),
@@ -126,9 +137,12 @@ AliFMDInput::AliFMDInput(const char* gAliceFile)
     fArrayS(0), 
     fArrayR(0), 
     fArrayA(0), 
+    fHeader(0),
     fGeoManager(0),
     fTreeMask(0), 
-    fIsInit(kFALSE)
+    fRawFile(""),
+    fIsInit(kFALSE),
+    fEventCount(0)
 {
   
   // Constructor of an FMD input object.  Specify what data to read in
@@ -142,6 +156,7 @@ Int_t
 AliFMDInput::NEvents() const 
 {
   // Get number of events
+  if (TESTBIT(fTreeMask, kRaw)) return -1;
   if (fTreeE) return fTreeE->GetEntries();
   return -1;
 }
@@ -157,36 +172,38 @@ AliFMDInput::Init()
     AliWarning("Already initialized");
     return fIsInit;
   }
-  if (fGAliceFile.IsNull()) fGAliceFile = "galice.root";
-  // Get the loader
-  fLoader = AliRunLoader::Open(fGAliceFile.Data(), "Alice", "read");
-  if (!fLoader) {
-    AliError(Form("Coulnd't read the file %s", fGAliceFile.Data()));
-    return kFALSE;
-  }
-  
   // Get the run 
-  if  (fLoader->LoadgAlice()) return kFALSE;
-  fRun = fLoader->GetAliRun();
+  if (!TESTBIT(fTreeMask, kRaw)) { 
+    if (fGAliceFile.IsNull()) fGAliceFile = "galice.root";
+    // Get the loader
+    fLoader = AliRunLoader::Open(fGAliceFile.Data(), "Alice", "read");
+    if (!fLoader) {
+      AliError(Form("Coulnd't read the file %s", fGAliceFile.Data()));
+      return kFALSE;
+    }
   
-  // Get the FMD 
-  fFMD = static_cast<AliFMD*>(fRun->GetDetector("FMD"));
-  if (!fFMD) {
-    AliError("Failed to get detector FMD from loader");
-    return kFALSE;
-  }
+    if  (fLoader->LoadgAlice()) return kFALSE;
+    fRun = fLoader->GetAliRun();
   
-  // Get the FMD loader
-  fFMDLoader = fLoader->GetLoader("FMDLoader");
-  if (!fFMDLoader) {
-    AliError("Failed to get detector FMD loader from loader");
-    return kFALSE;
-  }
-  if (fLoader->LoadHeader()) { 
-    AliError("Failed to get event header information from loader");
-    return kFALSE;
+    // Get the FMD 
+    fFMD = static_cast<AliFMD*>(fRun->GetDetector("FMD"));
+    if (!fFMD) {
+      AliError("Failed to get detector FMD from loader");
+      return kFALSE;
+    }
+  
+    // Get the FMD loader
+    fFMDLoader = fLoader->GetLoader("FMDLoader");
+    if (!fFMDLoader) {
+      AliError("Failed to get detector FMD loader from loader");
+      return kFALSE;
+    }
+    if (fLoader->LoadHeader()) { 
+      AliError("Failed to get event header information from loader");
+      return kFALSE;
+    }
+    fTreeE = fLoader->TreeE();
   }
-  fTreeE = fLoader->TreeE();
 
   // Optionally, get the ESD files
   if (TESTBIT(fTreeMask, kESD)) {
@@ -213,21 +230,39 @@ AliFMDInput::Init()
   if (TESTBIT(fTreeMask, kRaw)) {
     AliInfo("Getting FMD raw data digits");
     fArrayA = new TClonesArray("AliFMDDigit");
-    fReader = new AliRawReaderFile(-1);
+#if 0
+    if (!fRawFile.IsNull() && fRawFile.EndsWith(".root"))
+      fReader = new AliRawReaderRoot(fRawFile.Data());
+    else if (!fRawFile.IsNull() && fRawFile.EndsWith(".raw"))
+      fReader = new AliRawReaderDate(fRawFile.Data());
+    else
+      fReader = new AliRawReaderFile(-1);
+#else
+    if(!fRawFile.IsNull()) 
+      fReader = AliRawReader::Create(fRawFile.Data());
+    else 
+      fReader = new AliRawReaderFile(-1);
+#endif
+    fFMDReader = new AliFMDRawReader(fReader, 0);
   }
   
   // Optionally, get the geometry 
   if (TESTBIT(fTreeMask, kGeometry)) {
-    TString fname(fRun->GetGeometryFileName());
-    if (fname.IsNull()) {
-      Warning("Init", "No file name for the geometry from AliRun");
-      fname = gSystem->DirName(fGAliceFile);
-      fname.Append("/geometry.root");
+    if (fRun) {
+      TString fname(fRun->GetGeometryFileName());
+      if (fname.IsNull()) {
+       Warning("Init", "No file name for the geometry from AliRun");
+       fname = gSystem->DirName(fGAliceFile);
+       fname.Append("/geometry.root");
+      }
+      fGeoManager = TGeoManager::Import(fname.Data());
+      if (!fGeoManager) {
+       Fatal("Init", "No geometry manager found");
+       return kFALSE;
+      }
     }
-    fGeoManager = TGeoManager::Import(fname.Data());
-    if (!fGeoManager) {
-      Fatal("Init", "No geometry manager found");
-      return kFALSE;
+    else { 
+      AliGeomManager::LoadGeometry();
     }
     AliCDBManager* cdb   = AliCDBManager::Instance();
     AliCDBEntry*   align = cdb->Get("FMD/Align/Data");
@@ -250,7 +285,7 @@ AliFMDInput::Init()
     }
   }
 
-  
+  fEventCount = 0;
   fIsInit = kTRUE;
   return fIsInit;
 }
@@ -269,27 +304,37 @@ AliFMDInput::Begin(Int_t event)
     AliError("Not initialized");
     return fIsInit;
   }
+
   // Get the event 
-  if (fLoader->GetEvent(event)) return kFALSE;
-  AliInfo(Form("Now in event %d/%d", event, NEvents()));
+  if (fLoader && fLoader->GetEvent(event)) return kFALSE;
+  AliInfo(Form("Now in event %8d/%8d", event, NEvents()));
 
   // Possibly load global kinematics information 
   if (TESTBIT(fTreeMask, kKinematics) || TESTBIT(fTreeMask, kTracks)) {
-    AliInfo("Getting kinematics");
-    if (fLoader->LoadKinematics()) return kFALSE;
+    // AliInfo("Getting kinematics");
+    if (fLoader->LoadKinematics("READ")) return kFALSE;
     fStack = fLoader->Stack();
   }
+
   // Possibly load FMD Hit information 
   if (TESTBIT(fTreeMask, kHits) || TESTBIT(fTreeMask, kTracks)) {
-    AliInfo("Getting FMD hits");
-    if (fFMDLoader->LoadHits()) return kFALSE;
+    // AliInfo("Getting FMD hits");
+    if (!fFMDLoader || fFMDLoader->LoadHits("READ")) return kFALSE;
     fTreeH = fFMDLoader->TreeH();
     if (!fArrayH) fArrayH = fFMD->Hits(); 
   }
+
+  // Possibly load heaedr information 
+  if (TESTBIT(fTreeMask, kHeader)) {
+    // AliInfo("Getting FMD hits");
+    if (!fLoader /* || fLoader->LoadHeader()*/) return kFALSE;
+    fHeader = fLoader->GetHeader();
+  }
+
   // Possibly load FMD Digit information 
   if (TESTBIT(fTreeMask, kDigits)) {
-    AliInfo("Getting FMD digits");
-    if (fFMDLoader->LoadDigits()) return kFALSE;
+    // AliInfo("Getting FMD digits");
+    if (!fFMDLoader || fFMDLoader->LoadDigits("READ")) return kFALSE;
     fTreeD = fFMDLoader->TreeD();
     if (fTreeD) {
       if (!fArrayD) fArrayD = fFMD->Digits();
@@ -299,27 +344,32 @@ AliFMDInput::Begin(Int_t event)
       AliWarning(Form("Failed to load FMD Digits"));
     } 
   }
+
   // Possibly load FMD Sdigit information 
   if (TESTBIT(fTreeMask, kSDigits)) {
-    AliInfo("Getting FMD summable digits");
-    if (fFMDLoader->LoadSDigits()) return kFALSE;
+    // AliInfo("Getting FMD summable digits");
+    if (!fFMDLoader || fFMDLoader->LoadSDigits("READ")) return kFALSE;
     fTreeS = fFMDLoader->TreeS();
     if (!fArrayS) fArrayS = fFMD->SDigits();
   }
+
   // Possibly load FMD RecPoints information 
   if (TESTBIT(fTreeMask, kRecPoints)) {
-    AliInfo("Getting FMD reconstructed points");
-    if (fFMDLoader->LoadRecPoints()) return kFALSE;
+    // AliInfo("Getting FMD reconstructed points");
+    if (!fFMDLoader || fFMDLoader->LoadRecPoints("READ")) return kFALSE;
     fTreeR = fFMDLoader->TreeR();
     if (!fArrayR) fArrayR = new TClonesArray("AliFMDRecPoint");
     fTreeR->SetBranchAddress("FMD",  &fArrayR);
-  }  // Possibly load FMD ESD information 
+  }  
+
+  // Possibly load FMD ESD information 
   if (TESTBIT(fTreeMask, kESD)) {
-    AliInfo("Getting FMD event summary data");
+    // AliInfo("Getting FMD event summary data");
     Int_t read = fChainE->GetEntry(event);
     if (read <= 0) return kFALSE;
     fESD = fESDEvent->GetFMDData();
     if (!fESD) return kFALSE;
+#if 0
     TFile* f = fChainE->GetFile();
     if (f) {
       TObject* o = f->GetStreamerInfoList()->FindObject("AliFMDMap");
@@ -330,16 +380,19 @@ AliFMDInput::Begin(Int_t event)
       }
     }
     // fESD->CheckNeedUShort(fChainE->GetFile());
+#endif
   }
+
   // Possibly load FMD Digit information 
   if (TESTBIT(fTreeMask, kRaw)) {
-    AliInfo("Getting FMD raw data digits");
+    // AliInfo("Getting FMD raw data digits");
     if (!fReader->NextEvent()) return kFALSE;
-    AliFMDRawReader r(fReader, 0);
+    // AliFMDRawReader r(fReader, 0);
     fArrayA->Clear();
-    r.ReadAdcs(fArrayA);
+    fFMDReader->ReadAdcs(fArrayA);
+    AliFMDDebug(1, ("Got a total of %d digits", fArrayA->GetEntriesFast()));
   }
-  
+  fEventCount++;
   return kTRUE;
 }
 
@@ -384,19 +437,23 @@ AliFMDInput::ProcessHits()
     AliError("No hit tree defined");
     return kFALSE;
   }
+  if (!fArrayH) {
+    AliError("No hit array defined");
+    return kFALSE;
+  }
+
   Int_t nTracks = fTreeH->GetEntries();
   for (Int_t i = 0; i < nTracks; i++) {
     Int_t hitRead  = fTreeH->GetEntry(i);
     if (hitRead <= 0) continue;
-    if (!fArrayH) {
-      AliError("No hit array defined");
-      return kFALSE;
-    }
+
     Int_t nHit = fArrayH->GetEntries();
     if (nHit <= 0) continue;
+  
     for (Int_t j = 0; j < nHit; j++) {
       AliFMDHit* hit = static_cast<AliFMDHit*>(fArrayH->At(j));
       if (!hit) continue;
+
       TParticle* track = 0;
       if (TESTBIT(fTreeMask, kKinematics) && fStack) {
        Int_t trackno = hit->Track();
@@ -413,7 +470,7 @@ Bool_t
 AliFMDInput::ProcessTracks()
 {
   // Read the hit tree, and pass each hit to the member function
-  // ProcessHit.
+  // ProcessTrack.
   if (!fStack) {
     AliError("No track tree defined");
     return kFALSE;
@@ -422,25 +479,32 @@ AliFMDInput::ProcessTracks()
     AliError("No hit tree defined");
     return kFALSE;
   }
+  if (!fArrayH) {
+    AliError("No hit array defined");
+    return kFALSE;
+  }
+
+  // Int_t nTracks = fStack->GetNtrack();
   Int_t nTracks = fTreeH->GetEntries();
   for (Int_t i = 0; i < nTracks; i++) {
-    TParticle* track = fStack->Particle(i);
+    Int_t      trackno = nTracks - i - 1;
+    TParticle* track   = fStack->Particle(trackno);
     if (!track) continue;
+
+    // Get the hits for this track. 
     Int_t hitRead  = fTreeH->GetEntry(i);
-    if (hitRead <= 0) continue;
-    if (!fArrayH) {
-      AliError("No hit array defined");
-      return kFALSE;
+    Int_t nHit     = fArrayH->GetEntries();
+    if (nHit == 0 || hitRead <= 0) { 
+      // Let user code see the track, even if there's no hits. 
+      if (!ProcessTrack(trackno, track, 0)) return kFALSE;
+      continue;
     }
-    Int_t nHit = fArrayH->GetEntries();
-    if (nHit <= 0) continue;
 
+    // Loop over the hits corresponding to this track.
     for (Int_t j = 0; j < nHit; j++) {
       AliFMDHit* hit = static_cast<AliFMDHit*>(fArrayH->At(j));
-      if (!hit) continue;
-      if (!ProcessTrack(i, track, hit)) return kFALSE;
-    }    
-    // if (!ProcessTrack(i, track, fArrayH)) return kFALSE;
+      if (!ProcessTrack(trackno, track, hit)) return kFALSE;
+    }   
   }
   return kTRUE;
 }
@@ -451,11 +515,21 @@ AliFMDInput::ProcessDigits()
 {
   // Read the digit tree, and pass each digit to the member function
   // ProcessDigit.
+  if (!fTreeD) {
+    AliError("No digit tree defined");
+    return kFALSE;
+  }
+  if (!fArrayD) {
+    AliError("No digit array defined");
+    return kFALSE;
+  }
+
   Int_t nEv = fTreeD->GetEntries();
   for (Int_t i = 0; i < nEv; i++) {
     Int_t digitRead  = fTreeD->GetEntry(i);
     if (digitRead <= 0) continue;
     Int_t nDigit = fArrayD->GetEntries();
+    AliFMDDebug(0, ("Got %5d digits for this event", nDigit));
     if (nDigit <= 0) continue;
     for (Int_t j = 0; j < nDigit; j++) {
       AliFMDDigit* digit = static_cast<AliFMDDigit*>(fArrayD->At(j));
@@ -472,11 +546,21 @@ AliFMDInput::ProcessSDigits()
 {
   // Read the summable digit tree, and pass each sumable digit to the
   // member function ProcessSdigit.
-  Int_t nEv = fTreeD->GetEntries();
+  if (!fTreeS) {
+    AliWarning("No sdigit tree defined");
+    return kTRUE; // Empty SDigits is fine
+  }
+  if (!fArrayS) {
+    AliWarning("No sdigit array defined");
+    return kTRUE; // Empty SDigits is fine
+  }
+
+  Int_t nEv = fTreeS->GetEntries();
   for (Int_t i = 0; i < nEv; i++) {
     Int_t sdigitRead  = fTreeS->GetEntry(i);
     if (sdigitRead <= 0) continue;
     Int_t nSdigit = fArrayS->GetEntries();
+    AliFMDDebug(0, ("Got %5d digits for this event", nSdigit));
     if (nSdigit <= 0) continue;
     for (Int_t j = 0; j < nSdigit; j++) {
       AliFMDSDigit* sdigit = static_cast<AliFMDSDigit*>(fArrayS->At(j));
@@ -493,11 +577,18 @@ AliFMDInput::ProcessRawDigits()
 {
   // Read the digit tree, and pass each digit to the member function
   // ProcessDigit.
+  if (!fArrayA) {
+    AliError("No raw digit array defined");
+    return kFALSE;
+  }
+
   Int_t nDigit = fArrayA->GetEntries();
   if (nDigit <= 0) return kTRUE;
   for (Int_t j = 0; j < nDigit; j++) {
     AliFMDDigit* digit = static_cast<AliFMDDigit*>(fArrayA->At(j));
     if (!digit) continue;
+    if (AliLog::GetDebugLevel("FMD","") >= 40 && j < 30) 
+      digit->Print();
     if (!ProcessRawDigit(digit)) return kFALSE;
   }    
   return kTRUE;
@@ -509,6 +600,15 @@ AliFMDInput::ProcessRecPoints()
 {
   // Read the reconstrcted points tree, and pass each reconstruction
   // object (AliFMDRecPoint) to either ProcessRecPoint.
+  if (!fTreeR) {
+    AliError("No recpoint tree defined");
+    return kFALSE;
+  }
+  if (!fArrayR) {
+    AliError("No recpoints array defined");
+    return kFALSE;
+  }
+
   Int_t nEv = fTreeR->GetEntries();
   for (Int_t i = 0; i < nEv; i++) {
     Int_t recRead  = fTreeR->GetEntry(i);
@@ -588,7 +688,7 @@ AliFMDInput::End()
     fFMDLoader->UnloadRecPoints();
     fTreeR = 0;
   }
-  AliInfo("Now out event");
+  // AliInfo("Now out event");
   return kTRUE;
 }
 
@@ -602,7 +702,7 @@ AliFMDInput::Run()
   if (!(retval = Init())) return retval;
 
   Int_t nEvents = NEvents();
-  for (Int_t event = 0; event < nEvents; event++) {
+  for (Int_t event = 0; nEvents < 0 || event < nEvents; event++) {
     if (!(retval = Begin(event))) break;
     if (!(retval = Event())) break;
     if (!(retval = End())) break;
@@ -612,6 +712,31 @@ AliFMDInput::Run()
   return retval;
 }
 
+//__________________________________________________________________
+TArrayF 
+AliFMDInput::MakeLogScale(Int_t n, Double_t min, Double_t max) 
+{
+  // Service function to define a logarithmic axis. 
+  // Parameters: 
+  //   n    Number of bins 
+  //   min  Minimum of axis 
+  //   max  Maximum of axis 
+  TArrayF bins(n+1);
+  bins[0]      = min;
+  if (n <= 20) {
+    for (Int_t i = 1; i < n+1; i++) bins[i] = bins[i-1] + (max-min)/n;
+    return bins;
+  }
+  Float_t dp   = n / TMath::Log10(max / min);
+  Float_t pmin = TMath::Log10(min);
+  for (Int_t i = 1; i < n+1; i++) {
+    Float_t p = pmin + i / dp;
+    bins[i]   = TMath::Power(10, p);
+  }
+  return bins;
+}
+
+
 
 //____________________________________________________________________
 //