Adding PiKP-only histograms and eliminating a number of switches where histograms...
[u/mrichter/AliRoot.git] / MUON / AliMUONDigitCalibrator.cxx
index b214d4b..86a4d62 100644 (file)
@@ -1,24 +1,27 @@
 /**************************************************************************
-* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
-*                                                                        *
-* Author: The ALICE Off-line Project.                                    *
-* Contributors are mentioned in the code where appropriate.              *
-*                                                                        *
-* Permission to use, copy, modify and distribute this software and its   *
-* documentation strictly for non-commercial purposes is hereby granted   *
-* without fee, provided that the above copyright notice appears in all   *
-* copies and that both the copyright notice and this permission notice   *
-* appear in the supporting documentation. The authors make no claims     *
-* about the suitability of this software for any purpose. It is          *
-* provided "as is" without express or implied warranty.                  *
-**************************************************************************/
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
 
 // $Id$
 
 #include "AliMUONDigitCalibrator.h"
 
+#include "AliCDBEntry.h"
+#include "AliCDBManager.h"
 #include "AliLog.h"
 #include "AliMUONCalibrationData.h"
+#include "AliMUONConstants.h"
 #include "AliMUONLogger.h"
 #include "AliMUONPadStatusMaker.h"
 #include "AliMUONPadStatusMapMaker.h"
 #include "AliMUONVStore.h"
 #include "AliMpBusPatch.h"
 #include "AliMpConstants.h"
+#include "AliMpCDB.h"
 #include "AliMpDDLStore.h"
 #include "AliMpDEIterator.h"
 #include "AliMpDetElement.h"
+#include "AliMpManuStore.h"
 
 //-----------------------------------------------------------------------------
 /// \class AliMUONDigitCalibrator
@@ -60,11 +65,10 @@ ClassImp(AliMUONDigitCalibrator)
 const Int_t AliMUONDigitCalibrator::fgkNoGain(0);
 const Int_t AliMUONDigitCalibrator::fgkGainConstantCapa(1);
 const Int_t AliMUONDigitCalibrator::fgkGain(2);
+const Int_t AliMUONDigitCalibrator::fgkInjectionGain(3);
 
 //_____________________________________________________________________________
-AliMUONDigitCalibrator::AliMUONDigitCalibrator(const AliMUONCalibrationData& calib,
-                                                                                                                                                                                        const AliMUONRecoParam* recoParams,
-                                               const char* calibMode)
+AliMUONDigitCalibrator::AliMUONDigitCalibrator(Int_t runNumber)
 : TObject(),
 fLogger(new AliMUONLogger(20000)),
 fStatusMaker(0x0),
@@ -75,16 +79,49 @@ fApplyGains(0),
 fCapacitances(0x0),
 fNumberOfBadPads(0),
 fNumberOfPads(0),
-fChargeSigmaCut(0)
+fChargeSigmaCut(0),
+fMask(0)
 {
   /// ctor
   
-  Ctor(calibMode,calib,recoParams);
+  AliMUONRecoParam* recoParam(0x0);
+  
+  AliCDBEntry* e = AliCDBManager::Instance()->Get("MUON/Calib/RecoParam",runNumber);
+  if (e)
+  {
+    TObject* o = e->GetObject();
+    if ( o->IsA() == TObjArray::Class() )
+    {
+      TObjArray* a = static_cast<TObjArray*>(o);
+//      a->SetOwner(kTRUE); // FIXME: this should be done but somehow makes the reco crash at the end at cleaning stage... investigate why ?
+      TIter next(a);
+      AliMUONRecoParam* p;
+      while ( ( p = static_cast<AliMUONRecoParam*>(next()) ))
+      {
+        if ( p->IsDefault()) recoParam = p;
+      }
+    }
+    else
+    {
+      recoParam = static_cast<AliMUONRecoParam*>(o);
+    }
+  }
+  if (!recoParam)
+  {
+    AliError("Cannot get the recoParam. Failing");
+    return;
+  }
+  
+  // OK. Now get all we need and work...
+  
+  AliMUONCalibrationData calib(runNumber);
+  
+  Ctor(calib,recoParam,kFALSE);
 }
 
 //_____________________________________________________________________________
 AliMUONDigitCalibrator::AliMUONDigitCalibrator(const AliMUONCalibrationData& calib,
-                                               const char* calibMode)
+                                               const AliMUONRecoParam* recoParams)
 : TObject(),
 fLogger(new AliMUONLogger(20000)),
 fStatusMaker(0x0),
@@ -95,22 +132,44 @@ fApplyGains(0),
 fCapacitances(0x0),
 fNumberOfBadPads(0),
 fNumberOfPads(0),
-fChargeSigmaCut(0)
+fChargeSigmaCut(0),
+fMask(0)
 {
   /// ctor
   
-  Ctor(calibMode,calib,0x0);
+  Ctor(calib,recoParams);
+}
+
+//_____________________________________________________________________________
+AliMUONDigitCalibrator::AliMUONDigitCalibrator(const AliMUONCalibrationData& calib, int /*b*/)
+: TObject(),
+fLogger(new AliMUONLogger(20000)),
+fStatusMaker(0x0),
+fStatusMapMaker(0x0),
+fPedestals(0x0),
+fGains(0x0),
+fApplyGains(0),
+fCapacitances(0x0),
+fNumberOfBadPads(0),
+fNumberOfPads(0),
+fChargeSigmaCut(0),
+fMask(0)
+{
+  /// ctor
+  
+  Ctor(calib,0x0);
 }
 
 //_____________________________________________________________________________
 void
-AliMUONDigitCalibrator::Ctor(const char* calibMode,
-                             const AliMUONCalibrationData& calib,
-                             const AliMUONRecoParam* recoParams)
+AliMUONDigitCalibrator::Ctor(const AliMUONCalibrationData& calib,
+                             const AliMUONRecoParam* recoParams,
+                             Bool_t deferredInitialization)
 {
   /// designated ctor
   
-  TString cMode(calibMode);
+  TString cMode("NOGAIN");
+  if (recoParams) cMode=recoParams->GetCalibrationMode();
   cMode.ToUpper();
   
   if ( cMode == "NOGAIN" ) 
@@ -128,49 +187,56 @@ AliMUONDigitCalibrator::Ctor(const char* calibMode,
     fApplyGains = fgkGain;
     AliInfo("Will apply gain correction, with measured capacitances");
   }
+  else if ( cMode == "INJECTIONGAIN")
+       {
+               fApplyGains = fgkInjectionGain;
+    AliInfo("Will apply injection gain correction, with EMELEC factory gains");
+       }  
   else
   {
-    AliError(Form("Invalid calib mode = %s. Will use NOGAIN instead",calibMode));
+    AliError(Form("Invalid calib mode = %s. Will use NOGAIN instead",cMode.Data()));
     fApplyGains = fgkNoGain;
   }
-       
+  
+  // Load mapping manu store
+  if ( ! AliMpCDB::LoadManuStore() ) {
+    AliFatal("Could not access manu store from OCDB !");
+  }
+  
   fStatusMaker = new AliMUONPadStatusMaker(calib);
   
   // Set default values, as loose as reasonable
-
+  
   fChargeSigmaCut = 3.0;
   
-       Int_t mask(0x8080); // reject pads where ped *or* hv are missing
+       fMask = 0x8080; // reject pads where ped *or* hv are missing
        
        if ( recoParams )
        {
     // if we have reco params, we use limits and cuts from there :
     
-               fStatusMaker->SetHVSt12Limits(recoParams->HVSt12LowLimit(),recoParams->HVSt12HighLimit());
-               fStatusMaker->SetHVSt345Limits(recoParams->HVSt345LowLimit(),recoParams->HVSt345HighLimit());
-               fStatusMaker->SetPedMeanLimits(recoParams->PedMeanLowLimit(),recoParams->PedMeanHighLimit());
-               fStatusMaker->SetPedSigmaLimits(recoParams->PedSigmaLowLimit(),recoParams->PedSigmaHighLimit());
-               fStatusMaker->SetGainA1Limits(recoParams->GainA1LowLimit(),recoParams->GainA1HighLimit());
-               fStatusMaker->SetGainA2Limits(recoParams->GainA2LowLimit(),recoParams->GainA2HighLimit());
-               fStatusMaker->SetGainThresLimits(recoParams->GainThresLowLimit(),recoParams->GainThresHighLimit());
-               
-    mask = recoParams->PadGoodnessMask();
+    fStatusMaker->SetLimits(*recoParams);
+    
+    fMask = recoParams->PadGoodnessMask();
                //WARNING : getting this mask wrong is a very effective way of getting
                //no digits at all out of this class ;-)
     
     fChargeSigmaCut = recoParams->ChargeSigmaCut();
        }
+  else
+  {
+    fLogger->Log("No RecoParam available");
+    fLogger->Log(Form("SigmaCut=%e",fChargeSigmaCut));
+  }
   
-  Bool_t deferredInitialization = kTRUE;
-  
-  fStatusMapMaker = new AliMUONPadStatusMapMaker(*fStatusMaker,mask,deferredInitialization);
+  fStatusMapMaker = new AliMUONPadStatusMapMaker(*fStatusMaker,fMask,deferredInitialization);
   
   fPedestals = calib.Pedestals();
-
+  
   fGains = calib.Gains(); // we get gains whatever the calibMode is, in order
   // to get the saturation value...
-
-  if ( fApplyGains == fgkGain ) 
+  
+  if ( fApplyGains == fgkGain || fApplyGains == fgkInjectionGain ) 
   {
     fCapacitances = calib.Capacitances();
   }
@@ -180,16 +246,25 @@ AliMUONDigitCalibrator::Ctor(const char* calibMode,
 AliMUONDigitCalibrator::~AliMUONDigitCalibrator()
 {
   /// dtor.
-  delete fStatusMaker;
-  delete fStatusMapMaker;
   
-  AliInfo("Summary of messages:");
-  fLogger->Print();
+  if ( fNumberOfPads > 0 ) 
+  {
+    if ( fStatusMaker ) 
+    {
+      fStatusMaker->Report(fMask);
+    }
+    
+    AliInfo("Summary of messages:");
 
-       AliInfo(Form("We have seen %g pads, and rejected %g (%7.2f %%)",
-                                                        fNumberOfPads,fNumberOfBadPads,
-                                                        ( fNumberOfPads > 0 ) ? fNumberOfBadPads*100.0/fNumberOfPads : 0 ));
-       
+    fLogger->Print();
+    
+    AliInfo(Form("We have seen %g pads, and rejected %g (%7.2f %%)",
+                 fNumberOfPads,fNumberOfBadPads,
+                 ( fNumberOfPads > 0 ) ? fNumberOfBadPads*100.0/fNumberOfPads : 0 ));
+       }
+
+  delete fStatusMaker;
+  delete fStatusMapMaker;
   delete fLogger;
 }
 
@@ -200,8 +275,10 @@ AliMUONDigitCalibrator::Calibrate(AliMUONVDigitStore& digitStore)
   /// Calibrate the digits contained in digitStore  
   TIter next(digitStore.CreateTrackerIterator());
   AliMUONVDigit* digit;
-  Int_t detElemId(-1);
-  Double_t nsigmas = fChargeSigmaCut;
+  
+  fStatusMapMaker->RefreshRejectProbabilities(); // this will do something only for simulations
+  // (and only for those simulations where the reject list contain probabilities which are
+  // different from zero or one)
   
   AliDebug(1,Form("# of digits = %d",digitStore.GetSize()));
   
@@ -214,48 +291,31 @@ AliMUONDigitCalibrator::Calibrate(AliMUONVDigitStore& digitStore)
     }
     
     digit->Calibrated(kTRUE);
+    digit->ChargeInFC(kTRUE);
     
-    if ( digit->DetElemId() != detElemId ) 
-    {
-      // Find out occupancy of that DE
-      detElemId = digit->DetElemId();
-      AliMpDetElement* de = AliMpDDLStore::Instance()->GetDetElement(detElemId);
-      Double_t nchannels = de->NofChannels();
-      Double_t occ = digitStore.GetSize(detElemId)/nchannels;
-      if ( occ > 0.05 ) 
-      {
-        nsigmas = 10.0; // enlarge (a lot) sigma cut if occupancy is high
-        // (which probably means zero suppression was not exactly OK).
-        fLogger->Log(Form("Will use %5.0f*sigma cut for DE %04d "
-                          "due to high occupancy",nsigmas,detElemId));
-      }
-      else
-      {
-        nsigmas = fChargeSigmaCut;
-      }
-    }
-
     Float_t charge(0.0);
     Int_t statusMap;
     Bool_t isSaturated(kFALSE);
     
+    ++fNumberOfPads;
+    
     Bool_t ok = IsValidDigit(digit->DetElemId(),digit->ManuId(),digit->ManuChannel(),&statusMap);
-
-    digit->SetStatusMap(statusMap);
     
+    digit->SetStatusMap(statusMap);
+
     if (ok)
     {
-      ++fNumberOfPads;
       charge = CalibrateDigit(digit->DetElemId(),digit->ManuId(),digit->ManuChannel(),
-                              digit->ADC(),nsigmas,&isSaturated);
+                              digit->ADC(),fChargeSigmaCut,&isSaturated);
     }
     else
     {
-      ++fNumberOfBadPads;
+      ++fNumberOfBadPads;      
     }
     
     digit->SetCharge(charge);
     digit->Saturated(isSaturated);
+    
   }
 }
 
@@ -267,8 +327,15 @@ AliMUONDigitCalibrator::CalibrateDigit(Int_t detElemId, Int_t manuId, Int_t manu
 
 {
   /// Calibrate one digit
+  /// Return the digit charge, in fC
   
-  
+  if ( nsigmas < 0 ) 
+  {
+    nsigmas = fChargeSigmaCut;
+  }
+
+  fLogger->Log(Form("ChargeSigmaCut used = %e",nsigmas));
+
   AliMUONVCalibParam* pedestal = static_cast<AliMUONVCalibParam*>
   (fPedestals->FindObject(detElemId,manuId));
   
@@ -295,38 +362,48 @@ AliMUONDigitCalibrator::CalibrateDigit(Int_t detElemId, Int_t manuId, Int_t manu
   }
   
   Float_t padc = adc-pedestal->ValueAsFloat(manuChannel,0);
-  Float_t charge(0);
-  Float_t capa(1.0);
   
-  if ( fApplyGains == fgkGainConstantCapa ) 
-  {
-    capa = 0.2; // pF
-  }
-  else if ( fApplyGains == fgkGain ) 
+       // Gain (mV/fC) = 1/(a0*capa) with a0~1.25 and capa~0.2 
+  Float_t charge(0);
+  Float_t capa(AliMUONConstants::DefaultCapa()); // capa = 0.2 and a0 = 1.25
+  Float_t a0(AliMUONConstants::DefaultA0());  // is equivalent to gain = 4 mV/fC
+  Float_t a1(0);
+  Float_t adc2mv(AliMUONConstants::DefaultADC2MV()); // 1 ADC channel = 0.61 mV
+  Float_t injGain(4); // By default the gain is set to 4 mV/fC
+  //
+  // Note that the ChargeMax (for one pad) is roughly 4096 * 0.61 mV/channel / 4 mV/fC = 625 fC
+
+  if ( fApplyGains == fgkGain || fApplyGains == fgkInjectionGain ) 
   {
-    AliMpDetElement* de = AliMpDDLStore::Instance()->GetDetElement(detElemId);
-    
-    Int_t serialNumber = de->GetManuSerialFromId(manuId);
+    Int_t serialNumber 
+    = AliMpManuStore::Instance()->GetManuSerial(detElemId, manuId);
     
     AliMUONVCalibParam* param = static_cast<AliMUONVCalibParam*>(fCapacitances->FindObject(serialNumber));
     
     if ( param )
     {
-      capa = param->ValueAsFloat(manuChannel);
+      capa = param->ValueAsFloat(manuChannel,0);
+                       injGain = param->ValueAsFloat(manuChannel,1);
+      if ( injGain < 0 ) 
+      {
+        fLogger->Log(Form("injGain is %e < 0 for serialNumber=%d",injGain,serialNumber));
+        return 0.0;
+      }
     }
     else
     {
-      fLogger->Log(Form("No capa found for serialNumber=%d",serialNumber));
-      capa = 0.0;
+      // If capa not found in the OCDB we exit
+           fLogger->Log(Form("No capa (injGain) found for serialNumber=%d",serialNumber));
+                       return 0.0;
     }
   }
   
   if ( padc > nsigmas*pedestal->ValueAsFloat(manuChannel,1) ) 
   {
-    if ( fApplyGains != fgkNoGain ) 
+    if ( fApplyGains == fgkGain || fApplyGains == fgkGainConstantCapa ) 
     {
-      Float_t a0 = gain->ValueAsFloat(manuChannel,0);
-      Float_t a1 = gain->ValueAsFloat(manuChannel,1);
+      a0 = gain->ValueAsFloat(manuChannel,0);
+      a1 = gain->ValueAsFloat(manuChannel,1);
       Int_t thres = gain->ValueAsInt(manuChannel,2);
       if ( padc < thres ) 
       {
@@ -336,35 +413,39 @@ AliMUONDigitCalibrator::CalibrateDigit(Int_t detElemId, Int_t manuId, Int_t manu
       {
         charge = a0*thres + a0*(padc-thres) + a1*(padc-thres)*(padc-thres);
       }
+                       charge *= capa*adc2mv;
     }
-    else
+    else if ( fApplyGains == fgkInjectionGain ) 
     {
-      charge = padc;
+                       
+      charge = padc*adc2mv/injGain;
     }
+               else
+               {
+                       charge = a0*padc*capa*adc2mv;
+               }
   }
   
-  charge *= capa;
-  
   if ( isSaturated ) 
   {
     Int_t saturation(3000);
-  
-    if ( gain )
+    
+    if ( gain && ( fApplyGains != fgkNoGain ) )
     {
       saturation = gain->ValueAsInt(manuChannel,4);
     }
-  
+    
     if ( padc >= saturation )
     {
       *isSaturated = kTRUE;
     }
     else
     {
-      isSaturated = kFALSE;
+      *isSaturated = kFALSE;
     }
   }
   
-  return charge;
+  return ( charge > 0.0 ? charge : 0.0 );
 }
 
 //_____________________________________________________________________________
@@ -375,6 +456,9 @@ AliMUONDigitCalibrator::IsValidDigit(Int_t detElemId, Int_t manuId, Int_t manuCh
 {
   /// Check if a given pad is ok or not.
   
+  // initialize the statusmap to dead by default
+  if (statusMap) *statusMap = AliMUONPadStatusMapMaker::SelfDeadMask();
+
   // First a protection against bad input parameters
   AliMpDetElement* de = AliMpDDLStore::Instance()->GetDetElement(detElemId);
   if (!de) return kFALSE; // not existing DE
@@ -393,7 +477,7 @@ AliMUONDigitCalibrator::IsValidDigit(Int_t detElemId, Int_t manuId, Int_t manuCh
   // ok, now we have a valid channel number, so let's see if that pad
   // behaves or not ;-)
   
-  Int_t sm = fStatusMapMaker->StatusMap(detElemId,manuId,manuChannel);
+  Int_t sm = StatusMap(detElemId,manuId,manuChannel);
   
   if (statusMap) *statusMap = sm;