]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - FMD/AliFMDDigitizer.cxx
New RAW I/O. I rolled my own, because I wasn't happy with the old
[u/mrichter/AliRoot.git] / FMD / AliFMDDigitizer.cxx
index b55bbb4c330e94a4393ca5ce1c8906d12f3fe3bb..64f491855ea1b67d44297acb9f49186b4d8493b0 100644 (file)
 #include <AliLog.h>            // ALILOG_H
 #include "AliFMDDigitizer.h"   // ALIFMDDIGITIZER_H
 #include "AliFMD.h"            // ALIFMD_H
+#include "AliFMDGeometry.h"    // ALIFMDGEOMETRY_H
+#include "AliFMDDetector.h"    // ALIFMDDETECTOR_H
+#include "AliFMDRing.h"                // ALIFMDRING_H
 #include "AliFMDHit.h"         // ALIFMDHIT_H
 #include "AliFMDDigit.h"       // ALIFMDDIGIT_H
+#include "AliFMDParameters.h"   // ALIFMDPARAMETERS_H
 #include <AliRunDigitizer.h>   // ALIRUNDIGITIZER_H
 #include <AliRun.h>            // ALIRUN_H
 #include <AliLoader.h>         // ALILOADER_H
 #include <AliRunLoader.h>      // ALIRUNLOADER_H
     
-//____________________________________________________________________
-ClassImp(AliFMDEdepMap);
-
 //====================================================================
-ClassImp(AliFMDBaseDigitizer);
+ClassImp(AliFMDBaseDigitizer)
+#if 0
+  ; // This is here to keep Emacs for indenting the next line
+#endif
 
 //____________________________________________________________________
 AliFMDBaseDigitizer::AliFMDBaseDigitizer()  
@@ -228,9 +232,6 @@ AliFMDBaseDigitizer::AliFMDBaseDigitizer(AliRunDigitizer* manager)
 {
   // Normal CTOR
   AliDebug(1," processed");
-  SetVA1MipRange();
-  SetAltroChannelSize();
-  SetSampleRate();
   SetShapingTime();
 }
 
@@ -246,9 +247,6 @@ AliFMDBaseDigitizer::AliFMDBaseDigitizer(const Char_t* name,
 {
   // Normal CTOR
   AliDebug(1," processed");
-  SetVA1MipRange();
-  SetAltroChannelSize();
-  SetSampleRate();
   SetShapingTime();
 }
 
@@ -263,10 +261,21 @@ Bool_t
 AliFMDBaseDigitizer::Init()
 {
   // Initialization
+  AliFMDParameters::Instance()->Init();
   return kTRUE;
 }
  
 
+//____________________________________________________________________
+UShort_t
+AliFMDBaseDigitizer::MakePedestal(UShort_t, 
+                                 Char_t, 
+                                 UShort_t, 
+                                 UShort_t) const 
+{ 
+  return 0; 
+}
+
 //____________________________________________________________________
 void
 AliFMDBaseDigitizer::SumContributions(AliFMD* fmd) 
@@ -277,7 +286,7 @@ AliFMDBaseDigitizer::SumContributions(AliFMD* fmd)
     Fatal("SumContributions", "no run loader");
   
   // Clear array of deposited energies 
-  fEdep.Clear();
+  fEdep.Reset();
   
   // Get the FMD loader 
   AliLoader* inFMD = fRunLoader->GetLoader("FMDLoader");
@@ -303,10 +312,12 @@ AliFMDBaseDigitizer::SumContributions(AliFMD* fmd)
   // Get number of entries in the tree 
   Int_t ntracks  = Int_t(hitsTree->GetEntries());
   
+  AliFMDParameters* param = AliFMDParameters::Instance();
+  Int_t read = 0;
   // Loop over the tracks in the 
   for (Int_t track = 0; track < ntracks; track++)  {
     // Read in entry number `track' 
-    hitsBranch->GetEntry(track);
+    read += hitsBranch->GetEntry(track);
     
     // Get the number of hits 
     Int_t nhits = fmdHits->GetEntries ();
@@ -321,15 +332,23 @@ AliFMDBaseDigitizer::SumContributions(AliFMD* fmd)
       UShort_t sector   = fmdHit->Sector();
       UShort_t strip    = fmdHit->Strip();
       Float_t  edep     = fmdHit->Edep();
+
+      // Check if strip is `dead' 
+      if (param->IsDead(detector, ring, sector, strip)) continue;
+      
+      // Give warning in case of double hit 
       if (fEdep(detector, ring, sector, strip).fEdep != 0)
-       AliDebug(1, Form("Double hit in %d%c(%d,%d)", 
+       AliDebug(5, Form("Double hit in %d%c(%d,%d)", 
                         detector, ring, sector, strip));
       
+      // Sum energy deposition
       fEdep(detector, ring, sector, strip).fEdep  += edep;
       fEdep(detector, ring, sector, strip).fN     += 1;
       // Add this to the energy deposited for this strip
     }  // hit loop
   } // track loop
+  AliDebug(1, Form("Size of cache: %d bytes, read %d bytes", 
+                  sizeof(fEdep), read));
 }
 
 //____________________________________________________________________
@@ -340,23 +359,17 @@ AliFMDBaseDigitizer::DigitizeHits(AliFMD* fmd) const
   // the energy signal to ADC counts, and store the created digit in
   // the digits array (AliFMD::fDigits)
   //
+  AliFMDGeometry* geometry = AliFMDGeometry::Instance();
+  
   TArrayI counts(3);
   for (UShort_t detector=1; detector <= 3; detector++) {
     // Get pointer to subdetector 
-    AliFMDSubDetector* det = 0;
-    switch (detector) {
-    case 1: det = fmd->GetFMD1(); break;
-    case 2: det = fmd->GetFMD2(); break;
-    case 3: det = fmd->GetFMD3(); break;
-    }
+    AliFMDDetector* det = geometry->GetDetector(detector);
     if (!det) continue;
     for (UShort_t ringi = 0; ringi <= 1; ringi++) {
+      Char_t ring = ringi == 0 ? 'I' : 'O';
       // Get pointer to Ring
-      AliFMDRing* r = 0;
-      switch (ringi) {
-      case 0: if (det->GetInner()) r = det->GetInner(); break;
-      case 1: if (det->GetOuter()) r = det->GetOuter(); break;
-      }
+      AliFMDRing* r = det->GetRing(ring);
       if (!r) continue;
       
       // Get number of sectors 
@@ -374,20 +387,18 @@ AliFMDBaseDigitizer::DigitizeHits(AliFMD* fmd) const
          // VA1_ALICE channel. 
          if (strip % 128 == 0) last = 0;
          
-         Float_t edep = fEdep(detector, r->GetId(), sector, strip).fEdep;
-         ConvertToCount(edep, last, r->GetSiThickness(), 
-                        fmd->GetSiDensity(), counts);
+         Float_t edep = fEdep(detector, ring, sector, strip).fEdep;
+         ConvertToCount(edep, last, detector, ring, sector, strip, counts);
          last = edep;
-         AddDigit(fmd, detector, r->GetId(), sector, strip, 
-                  edep, UShort_t(counts[0]), 
-                  Short_t(counts[1]), Short_t(counts[2]));
+         AddDigit(fmd, detector, ring, sector, strip, edep, 
+                  UShort_t(counts[0]), Short_t(counts[1]), 
+                  Short_t(counts[2]));
 #if 0
          // This checks if the digit created will give the `right'
          // number of particles when reconstructed, using a naiive
          // approach.  It's here only as a quality check - nothing
          // else. 
-         CheckDigit(fEdep(detector, r->GetId(), sector, strip).fEdep,
-                    fEdep(detector, r->GetId(), sector, strip).fN,
+         CheckDigit(digit, fEdep(detector, ring, sector, strip).fN,
                     counts);
 #endif
        } // Strip
@@ -400,8 +411,10 @@ AliFMDBaseDigitizer::DigitizeHits(AliFMD* fmd) const
 void
 AliFMDBaseDigitizer::ConvertToCount(Float_t   edep, 
                                    Float_t   last,
-                                   Float_t   siThickness, 
-                                   Float_t   siDensity, 
+                                   UShort_t  detector, 
+                                   Char_t    ring, 
+                                   UShort_t  sector, 
+                                   UShort_t  strip,
                                    TArrayI&  counts) const
 {
   // Convert the total energy deposited to a (set of) ADC count(s). 
@@ -441,32 +454,31 @@ AliFMDBaseDigitizer::ConvertToCount(Float_t   edep,
   // 
   //                  = E + (l - E) * ext(-B * t)
   // 
-  const Float_t mipEnergy = 1.664 * siThickness * siDensity;
-  const Float_t convf     = (1 / mipEnergy * Float_t(fAltroChannelSize) 
-                            / fVA1MipRange);
-  UShort_t ped            = MakePedestal();
-
+  AliFMDParameters* param = AliFMDParameters::Instance();
+  Float_t  convF          = 1/param->GetPulseGain(detector,ring,sector,strip);
+  UShort_t ped            = MakePedestal(detector,ring,sector,strip);
+  UInt_t   maxAdc         = param->GetAltroChannelSize();
+  UShort_t rate           = param->GetSampleRate(AliFMDParameters::kBaseDDL);
+  UShort_t size           = param->GetAltroChannelSize();
+  
   // In case we don't oversample, just return the end value. 
-  if (fSampleRate == 1) {
-    counts[0] = UShort_t(TMath::Min(edep * convf + ped, 
-                                   Float_t(fAltroChannelSize)));
+  if (rate == 1) {
+    counts[0] = UShort_t(TMath::Min(edep * convF + ped, Float_t(size)));
     return;
   }
-
+  
   // Create a pedestal 
-  Int_t   n = fSampleRate;
-  Float_t B = fShapingTime;
-  for (Ssiz_t i = 0; i < n;  i++) {
-    Float_t t = Float_t(i) / n;
-    Float_t s = edep + (last - edep) * TMath::Exp(-B * t);
-    counts[i] = UShort_t(TMath::Min(s * convf + ped, 
-                                   Float_t(fAltroChannelSize))); 
+  Float_t b = fShapingTime;
+  for (Ssiz_t i = 0; i < rate;  i++) {
+    Float_t t = Float_t(i) / rate;
+    Float_t s = edep + (last - edep) * TMath::Exp(-b * t);
+    counts[i] = UShort_t(TMath::Min(s * convF + ped, Float_t(maxAdc)));
   }
 }
 
 
 //====================================================================
-ClassImp(AliFMDDigitizer);
+ClassImp(AliFMDDigitizer)
 
 //____________________________________________________________________
 AliFMDDigitizer::AliFMDDigitizer()  
@@ -481,7 +493,6 @@ AliFMDDigitizer::AliFMDDigitizer(AliRunDigitizer* manager)
 {
   // Normal CTOR
   AliDebug(1," processed");
-  SetPedestal();
 }
 
 //____________________________________________________________________
@@ -498,7 +509,7 @@ AliFMDDigitizer::Exec(Option_t*)
   // Get the input loader 
   TString inFolder(fManager->GetInputFolderName(0));
   fRunLoader = 
-    AliRunLoader::GetRunLoader(fManager->GetInputFolderName(0));
+    AliRunLoader::GetRunLoader(inFolder.Data());
   if (!fRunLoader) {
     AliError("Can not find Run Loader for input stream 0");
     return;
@@ -507,7 +518,8 @@ AliFMDDigitizer::Exec(Option_t*)
   if (!fRunLoader->GetAliRun()) fRunLoader->LoadgAlice();
 
   // Get the AliFMD object 
-  AliFMD* fmd = static_cast<AliFMD*>(fRunLoader->GetAliRun()->GetDetector("FMD"));
+  AliFMD* fmd = 
+    static_cast<AliFMD*>(fRunLoader->GetAliRun()->GetDetector("FMD"));
   if (!fmd) {
     AliError("Can not get FMD from gAlice");
     return;
@@ -541,7 +553,9 @@ AliFMDDigitizer::Exec(Option_t*)
   fmd->MakeBranchInTree(digitTree, fmd->GetName(), &(digits), 4000, 0);
   // TBranch* digitBranch = digitTree->GetBranch(fmd->GetName());
   // Fill the tree 
-  digitTree->Fill();
+  Int_t write = 0;
+  write = digitTree->Fill();
+  AliDebug(1, Form("Wrote %d bytes to digit tree", write));
   
   // Write the digits to disk 
   outFMD->WriteDigits("OVERWRITE");
@@ -555,9 +569,16 @@ AliFMDDigitizer::Exec(Option_t*)
 
 //____________________________________________________________________
 UShort_t
-AliFMDDigitizer::MakePedestal() const 
+AliFMDDigitizer::MakePedestal(UShort_t  detector, 
+                             Char_t    ring, 
+                             UShort_t  sector, 
+                             UShort_t  strip) const 
 {
-  return UShort_t(TMath::Max(gRandom->Gaus(fPedestal, fPedestalWidth), 0.));
+  // Make a pedestal 
+  AliFMDParameters* param =AliFMDParameters::Instance();
+  Float_t           mean  =param->GetPedestal(detector,ring,sector,strip);
+  Float_t           width =param->GetPedestalWidth(detector,ring,sector,strip);
+  return UShort_t(TMath::Max(gRandom->Gaus(mean, width), 0.));
 }
 
 //____________________________________________________________________
@@ -572,30 +593,41 @@ AliFMDDigitizer::AddDigit(AliFMD*  fmd,
                          Short_t  count2, 
                          Short_t  count3) const
 {
-  fmd->AddDigit(detector, ring, sector, strip, count1, count2, count3);
+  // Add a digit
+  fmd->AddDigitByFields(detector, ring, sector, strip, count1, count2, count3);
 }
 
 //____________________________________________________________________
 void
-AliFMDDigitizer::CheckDigit(Float_t         /* edep */, 
+AliFMDDigitizer::CheckDigit(AliFMDDigit*    digit,
                            UShort_t        nhits,
                            const TArrayI&  counts) 
 {
-  Int_t integral = counts[0];
+  // Check that digit is consistent
+  AliFMDParameters* param = AliFMDParameters::Instance();
+  UShort_t          det   = digit->Detector();
+  Char_t            ring  = digit->Ring();
+  UShort_t          sec   = digit->Sector();
+  UShort_t          str   = digit->Strip();
+  Float_t           mean  = param->GetPedestal(det,ring,sec,str);
+  Float_t           width = param->GetPedestalWidth(det,ring,sec,str);
+  UShort_t          range = param->GetVA1MipRange();
+  UShort_t          size  = param->GetAltroChannelSize();
+  Int_t             integral = counts[0];
   if (counts[1] >= 0) integral += counts[1];
   if (counts[2] >= 0) integral += counts[2];
-  integral -= Int_t(fPedestal + 2 * fPedestalWidth);
+  integral -= Int_t(mean + 2 * width);
   if (integral < 0) integral = 0;
   
-  Float_t convf = Float_t(fVA1MipRange) / fAltroChannelSize;
-  Float_t mips  = integral * convf;
+  Float_t convF = Float_t(range) / size;
+  Float_t mips  = integral * convF;
   if (mips > Float_t(nhits) + .5 || mips < Float_t(nhits) - .5) 
     Warning("CheckDigit", "Digit -> %4.2f MIPS != %d +/- .5 hits", 
            mips, nhits);
 }
 
 //====================================================================
-ClassImp(AliFMDSDigitizer);
+ClassImp(AliFMDSDigitizer)
 
 //____________________________________________________________________
 AliFMDSDigitizer::AliFMDSDigitizer()  
@@ -627,6 +659,7 @@ AliFMDSDigitizer::AliFMDSDigitizer(const Char_t* headerFile,
 //____________________________________________________________________
 AliFMDSDigitizer::~AliFMDSDigitizer() 
 {
+  // Destructor
   AliLoader* loader = fRunLoader->GetLoader("FMDLoader");
   loader->CleanSDigitizer();
 }
@@ -688,7 +721,9 @@ AliFMDSDigitizer::AddDigit(AliFMD*  fmd,
                           Short_t  count2, 
                           Short_t  count3) const
 {
-  fmd->AddSDigit(detector, ring, sector, strip, edep, count1, count2, count3);
+  // Add a summable digit
+  fmd->AddSDigitByFields(detector, ring, sector, strip, edep, 
+                        count1, count2, count3); 
 }