]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - FMD/AliFMDBaseDigitizer.cxx
change default max distance R to 0.1
[u/mrichter/AliRoot.git] / FMD / AliFMDBaseDigitizer.cxx
index a04cd4fa2f8cb4dd6e0c842b1e6f23dae5cd0e78..7dcf7e416decb52782dd8768b4ce09d451f01a34 100644 (file)
 //                    -1 + B + exp(-B)
 //
 
+#include <TMath.h>
 #include <TTree.h>             // ROOT_TTree
 //#include <TRandom.h>         // ROOT_TRandom
-#include <AliLog.h>            // ALILOG_H
+// #include <AliLog.h>         // ALILOG_H
+#include "AliFMDDebug.h" // Better debug macros
 #include "AliFMDBaseDigitizer.h" // ALIFMDDIGITIZER_H
 #include "AliFMD.h"            // ALIFMD_H
 #include "AliFMDGeometry.h"    // ALIFMDGEOMETRY_H
 #include "AliFMDHit.h"         // ALIFMDHIT_H
 // #include "AliFMDDigit.h"    // ALIFMDDIGIT_H
 #include "AliFMDParameters.h"   // ALIFMDPARAMETERS_H
-// #include <AliRunDigitizer.h>        // ALIRUNDIGITIZER_H
+// #include <AliDigitizationInput.h>   // ALIRUNDIGITIZER_H
 //#include <AliRun.h>          // ALIRUN_H
 #include <AliLoader.h>         // ALILOADER_H
+#include <AliRun.h>            // ALILOADER_H
 #include <AliRunLoader.h>      // ALIRUNLOADER_H
+#include <TRandom.h>
     
 //====================================================================
 ClassImp(AliFMDBaseDigitizer)
@@ -220,28 +224,32 @@ ClassImp(AliFMDBaseDigitizer)
 
 //____________________________________________________________________
 AliFMDBaseDigitizer::AliFMDBaseDigitizer()  
-  : fRunLoader(0),
+  : fFMD(0),
+    fRunLoader(0),
     fEdep(AliFMDMap::kMaxDetectors, 
          AliFMDMap::kMaxRings, 
          AliFMDMap::kMaxSectors, 
          AliFMDMap::kMaxStrips),
-    fShapingTime(0)
+    fShapingTime(6),
+    fStoreTrackRefs(kTRUE), 
+    fIgnoredLabels(0)
 {
+  AliFMDDebug(1, ("Constructed"));
   // Default ctor - don't use it
 }
 
 //____________________________________________________________________
-AliFMDBaseDigitizer::AliFMDBaseDigitizer(AliRunDigitizer* manager) 
-  : AliDigitizer(manager, "AliFMDBaseDigitizer", "FMD Digitizer base class"), 
+AliFMDBaseDigitizer::AliFMDBaseDigitizer(AliDigitizationInput* digInput) 
+  : AliDigitizer(digInput, "AliFMDBaseDigitizer", "FMD Digitizer base class"), 
+    fFMD(0),
     fRunLoader(0),
-    fEdep(AliFMDMap::kMaxDetectors, 
-         AliFMDMap::kMaxRings, 
-         AliFMDMap::kMaxSectors, 
-         AliFMDMap::kMaxStrips), 
-    fShapingTime(0)
+    fEdep(0),        // nDet==0 means 51200 slots
+    fShapingTime(6),
+    fStoreTrackRefs(kTRUE), 
+    fIgnoredLabels(0)
 {
   // Normal CTOR
-  AliDebug(1," processed");
+  AliFMDDebug(1, ("Constructed"));
   SetShapingTime();
 }
 
@@ -249,14 +257,15 @@ AliFMDBaseDigitizer::AliFMDBaseDigitizer(AliRunDigitizer* manager)
 AliFMDBaseDigitizer::AliFMDBaseDigitizer(const Char_t* name, 
                                         const Char_t* title) 
   : AliDigitizer(name, title),
+    fFMD(0),
     fRunLoader(0),
-    fEdep(AliFMDMap::kMaxDetectors, 
-         AliFMDMap::kMaxRings, 
-         AliFMDMap::kMaxSectors
-         AliFMDMap::kMaxStrips)
+    fEdep(0),        // nDet==0 means 51200 slots
+    fShapingTime(6),
+    fStoreTrackRefs(kTRUE)
+    fIgnoredLabels(0)
 {
   // Normal CTOR
-  AliDebug(1," processed");
+  AliFMDDebug(1, (" Constructed"));
   SetShapingTime();
 }
 
@@ -266,130 +275,140 @@ AliFMDBaseDigitizer::~AliFMDBaseDigitizer()
   // Destructor
 }
 
+//____________________________________________________________________
+AliFMDBaseDigitizer&
+AliFMDBaseDigitizer::operator=(const AliFMDBaseDigitizer& o) 
+{ 
+  // 
+  // Assignment operator
+  // 
+  // Return:
+  //    Reference to this object 
+  //
+  if (&o == this) return *this; 
+  AliDigitizer::operator=(o);
+  fRunLoader      = o.fRunLoader;
+  fEdep           = o.fEdep;
+  fShapingTime    = o.fShapingTime;
+  fStoreTrackRefs = o.fStoreTrackRefs;
+  fIgnoredLabels  = o.fIgnoredLabels;
+  return *this; 
+}
+
 //____________________________________________________________________
 Bool_t 
 AliFMDBaseDigitizer::Init()
 {
-  // Initialization
+  // Initialization.   Get a pointer to the parameter manager, and
+  // initialize it.  
   AliFMDParameters::Instance()->Init();
+  if (AliLog::GetDebugLevel("FMD","") >= 15) 
+    AliFMDParameters::Instance()->Print("");
   return kTRUE;
 }
 
 //____________________________________________________________________
 UShort_t
-AliFMDBaseDigitizer::MakePedestal(UShort_t, 
-                                 Char_t, 
-                                 UShort_t, 
-                                 UShort_t) const 
+AliFMDBaseDigitizer::MakePedestal(UShort_t detector
+                                 Char_t   ring
+                                 UShort_t sector
+                                 UShort_t strip) const 
 { 
-  // Make a pedestal
-  return 0; 
+  // Make a pedestal.  The pedestal value is drawn from a Gaussian
+  // distribution.  The mean of the distribution is the measured
+  // pedestal, and the width is the measured noise. 
+  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.));
 }
 
 //____________________________________________________________________
 void
-AliFMDBaseDigitizer::SumContributions(AliFMD* fmd) 
+AliFMDBaseDigitizer::AddContribution(UShort_t detector, 
+                                    Char_t   ring, 
+                                    UShort_t sector, 
+                                    UShort_t strip, 
+                                    Float_t  edep, 
+                                    Bool_t   isPrimary,
+                                    Int_t    nTrack,
+                                    Int_t*   tracknos)
 {
-  // Sum energy deposited contributions from each hit in a cache
-  // (fEdep).  
-  if (!fRunLoader) 
-    Fatal("SumContributions", "no run loader");
-  
-  // Clear array of deposited energies 
-  fEdep.Reset();
-  
-  // Get the FMD loader 
-  AliLoader* inFMD = fRunLoader->GetLoader("FMDLoader");
-  // And load the hits 
-  inFMD->LoadHits("READ");
-  
-  // Get the tree of hits 
-  TTree* hitsTree = inFMD->TreeH();
-  if (!hitsTree)  {
-    // Try again 
-    inFMD->LoadHits("READ");
-    hitsTree = inFMD->TreeH();
+  // Add edep contribution from (detector,ring,sector,strip) to cache
+  AliFMDParameters* param = AliFMDParameters::Instance();
+  AliFMDDebug(10, ("Adding contribution %7.5f for FMD%d%c[%2d,%3d] "
+                 " from %d tracks (%s)", 
+                 edep,
+                 detector, 
+                 ring,
+                 sector, 
+                 strip, 
+                 nTrack, 
+                 (isPrimary ? "primary" : "secondary")));
+  // Check if strip is `dead' 
+  if (param->IsDead(detector, ring, sector, strip)) { 
+    AliFMDDebug(5, ("FMD%d%c[%2d,%3d] is marked as dead", 
+                   detector, ring, sector, strip));
+    return;
   }
+  // Check if strip is out-side read-out range 
+  // if (strip < minstrip || strip > maxstrip) {
+  //   AliFMDDebug(5, ("FMD%d%c[%2d,%3d] is outside range [%3d,%3d]", 
+  //               detector,ring,sector,strip,minstrip,maxstrip));
+  //   continue;
+  // }
   
-  // Get the FMD branch 
-  TBranch* hitsBranch = hitsTree->GetBranch("FMD");
-  if (hitsBranch) fmd->SetHitsAddressBranch(hitsBranch);
-  else            AliFatal("Branch FMD hit not found");
-  
-  // Get a list of hits from the FMD manager 
-  TClonesArray *fmdHits = fmd->Hits();
-  
-  // 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' 
-    read += hitsBranch->GetEntry(track);
-    
-    // Get the number of hits 
-    Int_t nhits = fmdHits->GetEntries ();
-    for (Int_t hit = 0; hit < nhits; hit++) {
-      // Get the hit number `hit'
-      AliFMDHit* fmdHit = 
-       static_cast<AliFMDHit*>(fmdHits->UncheckedAt(hit));
-      
-      // Extract parameters 
-      UShort_t detector = fmdHit->Detector();
-      Char_t   ring     = fmdHit->Ring();
-      UShort_t sector   = fmdHit->Sector();
-      UShort_t strip    = fmdHit->Strip();
-      Float_t  edep     = fmdHit->Edep();
-      // UShort_t minstrip = param->GetMinStrip(detector, ring, sector, strip);
-      // UShort_t maxstrip = param->GetMaxStrip(detector, ring, sector, strip);
-      // Check if strip is `dead' 
-      if (param->IsDead(detector, ring, sector, strip)) { 
-       AliDebug(5, Form("FMD%d%c[%2d,%3d] is marked as dead", 
-                        detector, ring, sector, strip));
-       continue;
-      }
-      // Check if strip is out-side read-out range 
-      // if (strip < minstrip || strip > maxstrip) {
-      //   AliDebug(5, Form("FMD%d%c[%2d,%3d] is outside range [%3d,%3d]", 
-      //                   detector,ring,sector,strip,minstrip,maxstrip));
-      //   continue;
-      // }
-       
-      // Give warning in case of double hit 
-      if (fEdep(detector, ring, sector, strip).fEdep != 0)
-       AliDebug(5, Form("Double hit in %d%c(%d,%d)", 
-                        detector, ring, sector, strip));
+  AliFMDEdepHitPair& entry = fEdep(detector, ring, sector, strip);
+
+  // Give warning in case of double sdigit 
+  if (entry.fEdep != 0)
+    AliFMDDebug(5, ("Double digit in FMD%d%c[%2d,%3d]", 
+                   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));
+  // Sum energy deposition
+  Int_t oldN  =  entry.fN;
+  entry.fEdep += edep;
+  entry.fN    += nTrack;
+  if (isPrimary) entry.fNPrim += nTrack;
+  if (fStoreTrackRefs) { 
+    if (entry.fLabels.fN < entry.fN) {
+      AliFMDDebug(15, ("== New label array size %d, was %d, added %d", 
+                      entry.fN, entry.fLabels.fN, nTrack));
+      entry.fLabels.Set(entry.fN);
+    }
+    for (Int_t i = 0; i < nTrack; i++) {
+      AliFMDDebug(15, ("=> Setting track label # %d", oldN+i));
+      entry.fLabels[oldN + i] = tracknos[i];
+      AliFMDDebug(15, ("<= Setting track label # %d", oldN+i));
+    }
+  }
+  AliFMDDebug(15,("Adding contribution %f to FMD%d%c[%2d,%3d] (%f) track %d", 
+                 edep, detector, ring, sector, strip,
+                 entry.fEdep, (nTrack > 0 ? tracknos[0] : -1)));
+  
 }
 
 //____________________________________________________________________
 void
-AliFMDBaseDigitizer::DigitizeHits(AliFMD* fmd) const
+AliFMDBaseDigitizer::DigitizeHits() const
 {
   // For the stored energy contributions in the cache (fEdep), convert
   // the energy signal to ADC counts, and store the created digit in
   // the digits array (AliFMD::fDigits)
   //
+  AliFMDDebug(5, ("Will now digitize all the summed signals"));
+  fIgnoredLabels = 0;
   AliFMDGeometry* geometry = AliFMDGeometry::Instance();
   
-  TArrayI counts(3);
+  TArrayI counts(4);
   for (UShort_t detector=1; detector <= 3; detector++) {
+    AliFMDDebug(10, ("Processing hits in FMD%d", detector));
     // Get pointer to subdetector 
     AliFMDDetector* det = geometry->GetDetector(detector);
     if (!det) continue;
     for (UShort_t ringi = 0; ringi <= 1; ringi++) {
       Char_t ring = ringi == 0 ? 'I' : 'O';
+      AliFMDDebug(10, (" Processing hits in FMD%d%c", detector,ring));
       // Get pointer to Ring
       AliFMDRing* r = det->GetRing(ring);
       if (!r) continue;
@@ -398,6 +417,8 @@ AliFMDBaseDigitizer::DigitizeHits(AliFMD* fmd) const
       UShort_t nSectors = UShort_t(360. / r->GetTheta());
       // Loop over the number of sectors 
       for (UShort_t sector = 0; sector < nSectors; sector++) {
+       AliFMDDebug(10, ("  Processing hits in FMD%d%c[%2d]", 
+                       detector,ring,sector));
        // Get number of strips 
        UShort_t nStrips = r->GetNStrips();
        // Loop over the stips 
@@ -409,12 +430,41 @@ AliFMDBaseDigitizer::DigitizeHits(AliFMD* fmd) const
          // VA1_ALICE channel. 
          if (strip % 128 == 0) last = 0;
          
-         Float_t edep = fEdep(detector, ring, sector, strip).fEdep;
+         const AliFMDEdepHitPair& entry  = fEdep(detector,ring,sector,strip);
+         Float_t                  edep   = entry.fEdep;
+         UShort_t                 ntot   = entry.fN;
+         UShort_t                 nprim  = entry.fNPrim;
+         const TArrayI&           labels = entry.fLabels;
+         if (edep > 0)
+           AliFMDDebug(15, ("Edep = %f for FMD%d%c[%2d,%3d]", 
+                            edep, detector, ring, sector, strip));
          ConvertToCount(edep, last, detector, ring, sector, strip, counts);
          last = edep;
-         AddDigit(fmd, detector, ring, sector, strip, edep, 
+         
+
+         // The following line was introduced - wrongly - by Peter
+         // Hristov.  It _will_ break the digitisation and the
+         // following reconstruction.  The behviour of the
+         // digitisation models exactly the front-end as it should
+         // (no matter what memory concuption it may entail).  The
+         // check should be on zero suppression, since that's what
+         // models the front-end - if zero suppression is turned on
+         // in the front-end, then we can suppress empty digits -
+         // otherwise we shoud never do that.  Note, that the line
+         // affects _both_ normal digitisation and digitisation for
+         // summable digits, since the condition is on the energy
+         // deposition and not on the actual number of counts.  If
+         // this line should go anywhere, it should be in the
+         // possible overloaded AliFMDSDigitizer::AddDigit - not
+         // here. 
+         // 
+         //   if (edep<=0) continue;
+         AddDigit(detector, ring, sector, strip, edep, 
                   UShort_t(counts[0]), Short_t(counts[1]), 
-                  Short_t(counts[2]));
+                  Short_t(counts[2]), Short_t(counts[3]), 
+                  ntot, nprim, labels);
+         AliFMDDebug(15, ("   Adding digit in FMD%d%c[%2d,%3d]=%d", 
+                          detector,ring,sector,strip,counts[0]));
 #if 0
          // This checks if the digit created will give the `right'
          // number of particles when reconstructed, using a naiive
@@ -427,6 +477,10 @@ AliFMDBaseDigitizer::DigitizeHits(AliFMD* fmd) const
       } // Sector 
     } // Ring 
   } // Detector 
+  if (fIgnoredLabels > 0) 
+    AliWarning(Form("%d track labels could not be associated with digits "
+                   "due to limited storage facilities in AliDigit", 
+                   fIgnoredLabels));
 }
 
 //____________________________________________________________________
@@ -477,7 +531,8 @@ AliFMDBaseDigitizer::ConvertToCount(Float_t   edep,
   //                  = E + (l - E) * ext(-B * t)
   // 
   AliFMDParameters* param = AliFMDParameters::Instance();
-  Float_t  convF          = 1./param->GetPulseGain(detector,ring,sector,strip);
+  Float_t  convF          = (param->GetDACPerMIP() / param->GetEdepMip() *
+                            param->GetPulseGain(detector,ring,sector,strip));
   Int_t    ped            = MakePedestal(detector,ring,sector,strip);
   Int_t    maxAdc         = param->GetAltroChannelSize()-1;
   if (maxAdc < 0) {
@@ -485,31 +540,98 @@ AliFMDBaseDigitizer::ConvertToCount(Float_t   edep,
     maxAdc = 1023;
   }
   UShort_t rate           = param->GetSampleRate(detector,ring,sector,strip);
-  if (rate < 1 || rate > 3) rate = 1;
-  
+  AliFMDDebug(15, ("Sample rate for FMD%d%c[%2d,%3d] = %d", 
+                  detector, ring, sector, strip, rate));
+  if (rate < 1 || rate > 4) {
+    AliWarning(Form("Invalid sample rate for for FMD%d%c[%2d,%3d] = %d", 
+                   detector, ring, sector, strip, rate));
+    rate = 1;
+  }
+
   // In case we don't oversample, just return the end value. 
   if (rate == 1) {
     Float_t    a = edep * convF + ped;
     if (a < 0) a = 0;
     counts[0]    = UShort_t(TMath::Min(a, Float_t(maxAdc)));
-    AliDebug(2, Form("FMD%d%c[%2d,%3d]: converting ELoss %f to "
+    AliFMDDebug(15, ("FMD%d%c[%2d,%3d]: converting ELoss %f to "
                     "ADC %4d (%f,%d)",
                     detector,ring,sector,strip,edep,counts[0],convF,ped));
     return;
   }
+
   
   // Create a pedestal 
   Float_t b = fShapingTime;
   for (Ssiz_t i = 0; i < rate;  i++) {
-    Float_t t  = Float_t(i) / rate;
+    Float_t t  = Float_t(i) / rate + 1./rate;
     Float_t s  = edep + (last - edep) * TMath::Exp(-b * t);
     Float_t a  = Int_t(s * convF + ped);
     if (a < 0) a = 0;
     counts[i]  = UShort_t(TMath::Min(a, Float_t(maxAdc)));
   }
+  AliFMDDebug(15, ("Converted edep = %f to ADC (%x,%x,%x,%x) "
+                  "[gain: %f=(%f/%f*%f), pedestal: %d, rate: %d]", 
+                  edep, counts[0], counts[1], counts[2], counts[3], 
+                  convF, param->GetDACPerMIP(),param->GetEdepMip(),
+                  param->GetPulseGain(detector,ring,sector,strip), 
+                  ped, rate));
 }
 
+//____________________________________________________________________
+void
+AliFMDBaseDigitizer::AddDigit(UShort_t        detector, 
+                             Char_t          ring,
+                             UShort_t        sector, 
+                             UShort_t        strip, 
+                             Float_t         /* edep */, 
+                             UShort_t        count1, 
+                             Short_t         count2, 
+                             Short_t         count3,
+                             Short_t         count4,
+                             UShort_t        ntot, 
+                             UShort_t        /* nprim */,
+                             const TArrayI&  refs) const
+{
+  // Add a digit or summable digit
+  fFMD->AddDigitByFields(detector, ring, sector, strip, 
+                        count1, count2, count3, count4, 
+                        ntot, fStoreTrackRefs ? refs.fArray : 0);
+  if (fStoreTrackRefs && ntot > 3) fIgnoredLabels += ntot - 3;
+}
 
+//____________________________________________________________________
+TTree*
+AliFMDBaseDigitizer::MakeOutputTree(AliLoader* loader)
+{
+  // Create output tree using loader.   If the passed loader differs
+  // from the currently set loader in the FMD object, reset the FMD
+  // loader to be the passed loader.   This is for the cases wher the
+  // output is different from the output. 
+  AliFMDDebug(5, ("Making digits tree"));
+  loader->LoadDigits("UPDATE"); // "RECREATE");
+  TTree* out = loader->TreeD();
+  if (!out) loader->MakeTree("D");
+  out = loader->TreeD(); 
+  if (out) { 
+    out->Reset();
+    if (loader != fFMD->GetLoader()) 
+      fFMD->SetLoader(loader);
+    fFMD->MakeBranch("D");
+  }
+  return out;
+}
+
+//____________________________________________________________________
+void
+AliFMDBaseDigitizer::StoreDigits(const AliLoader* loader)
+{
+  // Write the digits to disk 
+  AliFMDDebug(5, ("Storing %d digits",   fFMD->Digits()->GetEntries()));
+  loader->WriteDigits("OVERWRITE");
+  loader->UnloadDigits();
+  // Reset the digits in the AliFMD object 
+  fFMD->ResetDigits();
+}
 
 //____________________________________________________________________
 //