Fix range and binning on histograms
[u/mrichter/AliRoot.git] / FMD / AliFMDBaseDigitizer.cxx
index 4bc8eca..9ae092c 100644 (file)
 // #include <AliRunDigitizer.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)
@@ -222,28 +224,32 @@ ClassImp(AliFMDBaseDigitizer)
 
 //____________________________________________________________________
 AliFMDBaseDigitizer::AliFMDBaseDigitizer()  
-  : fRunLoader(0),
+  : fFMD(0),
+    fRunLoader(0),
     fEdep(AliFMDMap::kMaxDetectors, 
          AliFMDMap::kMaxRings, 
          AliFMDMap::kMaxSectors, 
          AliFMDMap::kMaxStrips),
-    fShapingTime(6)
+    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"), 
+    fFMD(0),
     fRunLoader(0),
-    fEdep(AliFMDMap::kMaxDetectors, 
-         AliFMDMap::kMaxRings, 
-         AliFMDMap::kMaxSectors, 
-         AliFMDMap::kMaxStrips), 
-    fShapingTime(6)
+    fEdep(0),        // nDet==0 means 51200 slots
+    fShapingTime(6),
+    fStoreTrackRefs(kTRUE), 
+    fIgnoredLabels(0)
 {
   // Normal CTOR
-  AliFMDDebug(1, (" processed"));
+  AliFMDDebug(1, ("Constructed"));
   SetShapingTime();
 }
 
@@ -251,15 +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),
-    fShapingTime(6)
+    fEdep(0),        // nDet==0 means 51200 slots
+    fShapingTime(6),
+    fStoreTrackRefs(kTRUE), 
+    fIgnoredLabels(0)
 {
   // Normal CTOR
-  AliFMDDebug(1, (" processed"));
+  AliFMDDebug(1, (" Constructed"));
   SetShapingTime();
 }
 
@@ -270,135 +276,138 @@ AliFMDBaseDigitizer::~AliFMDBaseDigitizer()
 }
 
 //____________________________________________________________________
+AliFMDBaseDigitizer&
+AliFMDBaseDigitizer::operator=(const AliFMDBaseDigitizer& o) 
+{ 
+  // 
+  // Assignment operator
+  // 
+  // Return:
+  //    Reference to this object 
+  //
+  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","") >= 10) 
-    AliFMDParameters::Instance()->Print("ALL");
+  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' 
-      AliFMDDebug(2, ("Hit in FMD%d%c[%2d,%3d]=%f",
-                     detector, ring, sector, strip, edep));
-      if (param->IsDead(detector, ring, sector, strip)) { 
-       AliFMDDebug(1, ("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) {
-      //   AliFMDDebug(5, ("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)
-       AliFMDDebug(5, ("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
-  AliFMDDebug(1, ("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(4);
   for (UShort_t detector=1; detector <= 3; detector++) {
-    AliFMDDebug(5, ("Processing hits in FMD%d", 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(5, (" Processing hits in FMD%d%c", detector,ring));
+      AliFMDDebug(10, (" Processing hits in FMD%d%c", detector,ring));
       // Get pointer to Ring
       AliFMDRing* r = det->GetRing(ring);
       if (!r) continue;
@@ -407,7 +416,7 @@ 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(5, ("  Processing hits in FMD%d%c[%2d]", 
+       AliFMDDebug(10, ("  Processing hits in FMD%d%c[%2d]", 
                        detector,ring,sector));
        // Get number of strips 
        UShort_t nStrips = r->GetNStrips();
@@ -420,10 +429,18 @@ 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;
          
+
          // The following line was introduced - wrongly - by Peter
          // Hristov.  It _will_ break the digitisation and the
          // following reconstruction.  The behviour of the
@@ -441,11 +458,12 @@ AliFMDBaseDigitizer::DigitizeHits(AliFMD* fmd) const
          // here. 
          // 
          //   if (edep<=0) continue;
-         AddDigit(fmd, detector, ring, sector, strip, edep, 
+         AddDigit(detector, ring, sector, strip, edep, 
                   UShort_t(counts[0]), Short_t(counts[1]), 
-                  Short_t(counts[2]), Short_t(counts[3]));
-         AliFMDDebug(10, ("   Adding digit in FMD%d%c[%2d,%3d]=%d", 
-                         detector,ring,sector,strip,counts[0]));
+                  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
@@ -458,6 +476,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));
 }
 
 //____________________________________________________________________
@@ -508,7 +530,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) {
@@ -516,18 +539,25 @@ AliFMDBaseDigitizer::ConvertToCount(Float_t   edep,
     maxAdc = 1023;
   }
   UShort_t rate           = param->GetSampleRate(detector,ring,sector,strip);
-  if (rate < 1 || rate > 4) 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)));
-    AliFMDDebug(2, ("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;
@@ -538,9 +568,69 @@ AliFMDBaseDigitizer::ConvertToCount(Float_t   edep,
     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();
+}
 
 //____________________________________________________________________
 //