Being more realistic with gain calibration (Laurent)
authorivana <ivana@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 28 Jun 2007 15:37:34 +0000 (15:37 +0000)
committerivana <ivana@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 28 Jun 2007 15:37:34 +0000 (15:37 +0000)
MUON/AliMUONCDB.cxx
MUON/AliMUONDigitCalibrator.cxx
MUON/AliMUONDigitizerV3.cxx
MUON/AliMUONGainSubprocessor.cxx

index e1ac639f93ffd441737059d2fccaa4dd5a0ff37e..608530de8100e7a8d8144df1fdd622d484fc2245 100644 (file)
@@ -65,6 +65,7 @@
 #include <TRandom.h>
 #include <TStopwatch.h>
 #include <TSystem.h>
+#include <TMath.h>
 
 /// \cond CLASSIMP
 ClassImp(AliMUONCDB)
@@ -112,6 +113,24 @@ void getBoundaries(const AliMUONVStore& store,
   }  
 }
 
+//_____________________________________________________________________________
+Double_t GetRandom(Double_t mean, Double_t sigma, Bool_t mustBePositive)
+{
+  Double_t x(-1);
+  if ( mustBePositive ) 
+  {
+    while ( x < 0 ) 
+    {
+      x = gRandom->Gaus(mean,sigma);
+    }
+  }
+  else
+  {
+    x = gRandom->Gaus(mean,sigma);
+  }
+  return x;
+}
+
 }
 
 //_____________________________________________________________________________
@@ -345,7 +364,7 @@ AliMUONCDB::MakeHVStore(TMap& aliasMap, Bool_t defaultValues)
       for ( UInt_t timeStamp = 0; timeStamp < 60*15; timeStamp += 120 )
       {
         Float_t value = 1500;
-        if (!defaultValues) value = gRandom->Gaus(1750,62.5);
+        if (!defaultValues) value = GetRandom(1750,62.5,true);
         AliDCSValue* dcsValue = new AliDCSValue(value,timeStamp);
         valueSet->Add(dcsValue);
       }
@@ -411,16 +430,9 @@ AliMUONCDB::MakePedestalStore(AliMUONVStore& pedestalStore, Bool_t defaultValues
       }
       else
       {
-        meanPedestal = -1;
-        while ( meanPedestal < 0 )
-        {
-          meanPedestal = gRandom->Gaus(kPedestalMeanMean,kPedestalMeanSigma);
-        }
-        sigmaPedestal = -1;
-        while ( sigmaPedestal < 0 )
-        {
-          sigmaPedestal = gRandom->Gaus(kPedestalSigmaMean,kPedestalSigmaSigma);
-        }
+        Bool_t positive(kTRUE);
+        meanPedestal = GetRandom(kPedestalMeanMean,kPedestalMeanSigma,positive);
+        sigmaPedestal = GetRandom(kPedestalSigmaMean,kPedestalSigmaSigma,positive);
       }
       ped->SetValueAsFloat(manuChannel,0,meanPedestal);
       ped->SetValueAsFloat(manuChannel,1,sigmaPedestal);
@@ -501,11 +513,7 @@ AliMUONCDB::MakeCapacitanceStore(AliMUONVStore& capaStore, Bool_t defaultValues)
       }
       else
       {
-        capaValue = -1;
-        while ( capaValue < 0 )
-        {
-          capaValue = gRandom->Gaus(kCapaMean,kCapaSigma);
-        }
+        capaValue = GetRandom(kCapaMean,kCapaSigma,kTRUE);
       }
       capa->SetValueAsFloat(manuChannel,0,capaValue);
     }
@@ -528,9 +536,11 @@ AliMUONCDB::MakeCapacitanceStore(AliMUONVStore& capaStore, Bool_t defaultValues)
 Int_t 
 AliMUONCDB::MakeGainStore(AliMUONVStore& gainStore, Bool_t defaultValues)
 {  
-  /// Create a gain store. if defaultValues=true, all gain are 1.0,
-  /// otherwise they are from a gaussian with parameters defined in the
-  /// kGain* constants below.
+  /// Create a gain store. if defaultValues=true, all gains set so that
+  /// charge = (adc-ped)
+  ///
+  /// otherwise parameters are taken from gaussians with parameters 
+  /// defined in the k* constants below.
   
   TIter next(ManuList());
   
@@ -539,9 +549,15 @@ AliMUONCDB::MakeGainStore(AliMUONVStore& gainStore, Bool_t defaultValues)
   Int_t nchannels(0);
   Int_t nmanus(0);
     
-  const Double_t kSaturation(3000);
-  const Double_t kGainMean(1.0);
-  const Double_t kGainSigma(0.05);
+  const Int_t kSaturation(3000);
+  const Double_t kA0Mean(1.2);
+  const Double_t kA0Sigma(0.1);
+  const Double_t kA1Mean(1E-5);
+  const Double_t kA1Sigma(1E-6);
+  const Double_t kQualMean(0xFF);
+  const Double_t kQualSigma(0x10);
+  const Int_t kThresMean(1600);
+  const Int_t kThresSigma(100);
   
   while ( ( p = (AliMpIntPair*)next() ) )
   {
@@ -551,7 +567,7 @@ AliMUONCDB::MakeGainStore(AliMUONVStore& gainStore, Bool_t defaultValues)
     Int_t manuId = p->GetSecond();
 
     AliMUONVCalibParam* gain = 
-      new AliMUONCalibParamNF(2,AliMpConstants::ManuNofChannels(),
+      new AliMUONCalibParamNF(5,AliMpConstants::ManuNofChannels(),
                               detElemId,
                               manuId,
                               AliMUONVCalibParam::InvalidFloatValue());
@@ -567,23 +583,23 @@ AliMUONCDB::MakeGainStore(AliMUONVStore& gainStore, Bool_t defaultValues)
       
       ++nchannels;
       
-      Float_t meanGain;
-      Float_t saturation(kSaturation);
-    
       if ( defaultValues ) 
       {
-        meanGain = 1.0;
+        gain->SetValueAsFloat(manuChannel,0,1.0);
+        gain->SetValueAsFloat(manuChannel,1,0.0);
+        gain->SetValueAsInt(manuChannel,2,4095); 
+        gain->SetValueAsInt(manuChannel,3,1);
+        gain->SetValueAsInt(manuChannel,4,kSaturation);
       }
       else
       {
-        meanGain = -1;
-        while ( meanGain < 0 )
-        {
-          meanGain = gRandom->Gaus(kGainMean,kGainSigma);
-        }
+        Bool_t positive(kTRUE);
+        gain->SetValueAsFloat(manuChannel,0,GetRandom(kA0Mean,kA0Sigma,positive));
+        gain->SetValueAsFloat(manuChannel,1,GetRandom(kA1Mean,kA1Sigma,!positive));
+        gain->SetValueAsInt(manuChannel,2,(Int_t)TMath::Nint(GetRandom(kThresMean,kThresSigma,positive)));
+        gain->SetValueAsInt(manuChannel,3,(Int_t)TMath::Nint(GetRandom(kQualMean,kQualSigma,positive)));
+        gain->SetValueAsInt(manuChannel,4,kSaturation);        
       }
-      gain->SetValueAsFloat(manuChannel,0,meanGain);
-      gain->SetValueAsFloat(manuChannel,1,saturation);
       
     }
     Bool_t ok = gainStore.Add(gain);
index b0c2ae4330e8a662abf89f7ad1fc7f0e76dc170f..7ff75144ab80b0bb702d49f8104942f6feeaebc5 100644 (file)
@@ -129,6 +129,12 @@ AliMUONDigitCalibrator::CalibrateDigit(AliMUONVDigit& digit)
 {
   /// Calibrate one digit
   
+  if ( digit.IsCalibrated() ) 
+  {
+    fLogger->Log("ERROR : trying to calibrate a digit twice");
+    return;
+  }
+  
   AliMUONVCalibParam* deadmap = static_cast<AliMUONVCalibParam*>
   (fStatusMap->FindObject(digit.DetElemId(),digit.ManuId()));
   Int_t statusMap = deadmap->ValueAsInt(digit.ManuChannel());
@@ -168,13 +174,23 @@ AliMUONDigitCalibrator::CalibrateDigit(AliMUONVDigit& digit)
     Int_t manuChannel = digit.ManuChannel();
     Float_t adc = digit.ADC();
     Float_t padc = adc-pedestal->ValueAsFloat(manuChannel,0);
-    if ( padc < 3.0*pedestal->ValueAsFloat(manuChannel,1) ) 
+    Float_t charge(0);
+    if ( padc > 3.0*pedestal->ValueAsFloat(manuChannel,1) ) 
     {
-      padc = 0.0;
+      Float_t a0 = gain->ValueAsFloat(manuChannel,0);
+      Float_t a1 = gain->ValueAsFloat(manuChannel,1);
+      Int_t thres = gain->ValueAsInt(manuChannel,2);
+      if ( padc < thres ) 
+      {
+        charge = a0*padc;
+      }
+      else
+      {
+        charge = a0*thres + a0*padc + a1*padc*padc;
+      }
     }
-    Float_t charge = padc*gain->ValueAsFloat(manuChannel,0);
     digit.SetCharge(charge);
-    Int_t saturation = gain->ValueAsInt(manuChannel,1);
+    Int_t saturation = gain->ValueAsInt(manuChannel,4);
     if ( charge >= saturation )
     {
       digit.Saturated(kTRUE);
index 2075c76b001530dc04a9d971031a9223e112f2e0..d391b4002abb38edf2e1e0fc4a21b7b7748fd36c 100644 (file)
@@ -150,17 +150,21 @@ AliMUONDigitizerV3::ApplyResponseToTrackerDigit(AliMUONVDigit& digit, Bool_t add
   /// For tracking digits, starting from an ideal digit's charge, we :
   ///
   /// - add some noise (thus leading to a realistic charge), if requested to do so
-  /// - divide by a gain (thus decalibrating the digit)
+  /// - "divide" by a gain (thus decalibrating the digit)
   /// - add a pedestal (thus decalibrating the digit)
   /// - sets the signal to zero if below 3*sigma of the noise
 
   static const Int_t kMaxADC = (1<<12)-1; // We code the charge on a 12 bits ADC.
   
-  Float_t signal = digit.Charge();
+  Float_t charge = digit.Charge();
+  
+  // We set the charge to 0, as the only relevant piece of information
+  // after Digitization is the ADC value.  
+  digit.SetCharge(0);
   
   if ( !addNoise )
   {
-    digit.SetADC(TMath::Min(kMaxADC,TMath::Nint(signal)));
+    digit.SetADC(TMath::Min(kMaxADC,TMath::Nint(charge)));
     return;
   }
   
@@ -175,7 +179,6 @@ AliMUONDigitizerV3::ApplyResponseToTrackerDigit(AliMUONVDigit& digit, Bool_t add
     fLogger->Log(Form("%s:%d:Could not get pedestal for DE=%4d manuId=%4d. Disabling.",
                       __FILE__,__LINE__,
                       detElemId,manuId));
-    digit.SetCharge(0);
     digit.SetADC(0);
     return;    
   }
@@ -188,32 +191,40 @@ AliMUONDigitizerV3::ApplyResponseToTrackerDigit(AliMUONVDigit& digit, Bool_t add
     fLogger->Log(Form("%s:%d:Could not get gain for DE=%4d manuId=%4d. Disabling.",
                       __FILE__,__LINE__,
                       detElemId,manuId));
-    digit.SetCharge(0);
     digit.SetADC(0);
     return;        
   }    
-  Float_t gainMean = gain->ValueAsFloat(manuChannel,0);
-  
-  Float_t adcNoise = gRandom->Gaus(0.0,pedestalSigma);
+
+  Float_t a0 = gain->ValueAsFloat(manuChannel,0);
+  Float_t a1 = gain->ValueAsFloat(manuChannel,1);
+
+  Int_t thres = gain->ValueAsInt(manuChannel,2);
+  Float_t chargeThres = a0*thres;
   
-  Int_t adc;
+  Float_t padc(0); // (adc - ped) value
   
-  if ( gainMean < 1E-6 )
-  {
-    AliError(Form("Got a too small gain %e for DE=%d manuId=%d manuChannel=%d. "
-                  "Setting signal to 0.",
-                  gainMean,detElemId,manuId,manuChannel));
-    adc = 0;
-  }
-  else
+  if ( charge <= chargeThres || TMath::Abs(a1) < 1E-12 ) 
   {
-    adc = TMath::Nint( signal / gainMean + pedestalMean + adcNoise);///
+    // linear part only
     
-    if ( adc <= pedestalMean + fgkNSigmas*pedestalSigma ) 
+    if ( TMath::Abs(a0) > 1E-12 ) 
     {
-      adc = 0;
+      padc = charge/a0;    
     }
   }
+  else // charge > chargeThres && a1 not zero
+  {
+    // parabolic part
+    padc = TMath::Sqrt((chargeThres-charge)/a1) + thres;
+  }
+
+  Float_t adcNoise = gRandom->Gaus(0.0,pedestalSigma);
+  Int_t adc(0);
+  
+  if ( padc > 0 ) 
+  {
+    adc = TMath::Nint(padc + pedestalMean + adcNoise);
+  }
   
   // be sure we stick to 12 bits.
   if ( adc > kMaxADC )
@@ -221,7 +232,14 @@ AliMUONDigitizerV3::ApplyResponseToTrackerDigit(AliMUONVDigit& digit, Bool_t add
     adc = kMaxADC;
   }
   
-  digit.SetCharge(adc);
+  AliDebug(3,Form("DE %4d Manu %4d Ch %2d Charge %e A0 %e A1 %e Thres %d padc %e ADC %4d",
+                  detElemId,manuId,manuChannel,charge,
+                  gain->ValueAsFloat(manuChannel,0),
+                  gain->ValueAsFloat(manuChannel,1),
+                  gain->ValueAsInt(manuChannel,2),
+                  padc,
+                  adc));
+                  
   digit.SetADC(adc);
 }
 
@@ -309,7 +327,7 @@ AliMUONDigitizerV3::ApplyResponse(const AliMUONVDigitStore& store,
     {
       ApplyResponseToTriggerDigit(store,*digit);
     }
-    if ( digit->Charge() > 0  )
+    if ( digit->ADC() > 0  || digit->Charge() > 0 )
     {
       filteredStore.Add(*digit,AliMUONVDigitStore::kIgnore);
     }
@@ -558,7 +576,7 @@ AliMUONDigitizerV3::GenerateNoisyDigitsForOneCathode(AliMUONVDigitStore& digitSt
     d->SetCharge(TMath::Nint(ped+pedestalMean+0.5));
     d->NoiseOnly(kTRUE);
     ApplyResponseToTrackerDigit(*d,kFALSE);
-    if ( d->Charge() > 0 )
+    if ( d->ADC() > 0 )
     {
       Bool_t ok = digitStore.Add(*d,AliMUONVDigitStore::kDeny);
       // this can happen (that we randomly chose a digit that is
index 285f9f8b2b5e0de59fe2908f692c8814741ac252..c3150024260a57a1fd2d18f5b84e9fdceb91a7e5 100644 (file)
 ///
 /// Implementation of AliMUONVSubprocessor class to deal with MUON TRK Gains.
 ///
-/// Gains are read in from an ascii file, with the format :                   \n            
-///---------------------------------------------------------------------------\n
-///BUS_PATCH   MANU   CHANNEL    Ped.     a0        a1         a2         xlim        P(chi2)    P(chi2)_2  \n 
-///---------------------------------------------------------------------------\n
+/// Gains are read in from an ascii file, with the format :                               
+///
+///---------------------------------------------------------------------------
+///
+///BUS_PATCH   MANU   CHANNEL    a0   a1   thres Qual
+///
+///---------------------------------------------------------------------------
 ///
 /// \author L. Aphecetche
 ///
@@ -164,13 +167,13 @@ AliMUONGainSubprocessor::Process(TMap* /*dcsAliasMap*/)
 Int_t
 AliMUONGainSubprocessor::ReadFile(const char* filename)
 {
-  /// Read the Gains from an ASCII file.                                  \n
-  /// Format of that file is one line per channel :                           \n
-  ///-------------------------------------------------------------------------\n
-  ///BUS_PATCH   MANU   CHANNEL    Ped.     a0        a1         a2         xlim        P(chi2)    P(chi2)_2  \n 
-  ///-------------------------------------------------------------------------\n
-  ///                                                                         \n
-  /// Return kFALSE if reading was not successfull.                           \n
+  /// Read the Gains from an ASCII file.                                  
+  /// Format of that file is one line per channel :                       
+  ///-------------------------------------------------------------------------
+  ///BUS_PATCH   MANU   CHANNEL  a0  a1 thres Qual
+  ///-------------------------------------------------------------------------
+  ///                                                                         
+  /// Return kFALSE if reading was not successfull.                           
   ///
   
   TString sFilename(gSystem->ExpandPathName(filename));
@@ -184,10 +187,10 @@ AliMUONGainSubprocessor::ReadFile(const char* filename)
   }
   char line[256];
   Int_t busPatchID, manuID, manuChannel;
-  Float_t ped;
-  Float_t a0, a1, a2;
-  Float_t xlim;
-  Float_t chi2, chi22;
+  Float_t a0, a1;
+  Int_t thres;
+  Int_t qual;
+  const Int_t kSaturation(3000); // FIXME: how to get this number ?
   
   static const Int_t kNchannels(AliMpConstants::ManuNofChannels());
   Int_t n(0);
@@ -196,29 +199,29 @@ AliMUONGainSubprocessor::ReadFile(const char* filename)
   {
     if ( strlen(line) < 10 ) continue;
     if ( line[0] == '/' && line[1] == '/' ) continue;
-    std::istringstream sin(line);
+
+    sscanf(line,"%d %d %d %f %f %d %x",&busPatchID,&manuID,&manuChannel,
+           &a0,&a1,&thres,&qual); 
     AliDebug(3,Form("line=%s",line));
-    sin >> busPatchID >> manuID >> manuChannel >> ped >> a0 >> a1 >> a2 >> xlim >> chi2 >> chi22;
     Int_t detElemID = AliMpDDLStore::Instance()->GetDEfromBus(busPatchID);
-    AliDebug(3,Form("BUSPATCH %3d DETELEMID %4d MANU %3d CH %3d PED %7.2f A0 %7.2f A1 %7.2f A2 %7.2f"
-                    " XLIM %7.2f CHI2 %7.2f CHI22 %7.2f",
-                    busPatchID,detElemID,manuID,manuChannel,ped,a0,a1,a2,xlim,chi2,chi22));
-    if ( a0==a1 && a1==a2 && a0==-2) continue;
+    AliDebug(3,Form("BUSPATCH %3d DETELEMID %4d MANU %3d CH %3d A0 %7.2f "
+                    "A1 %e THRES %5d QUAL %x",
+                    busPatchID,detElemID,manuID,manuChannel,a0,a1,thres,qual));
+    if ( qual == 0 ) continue;
     
     AliMUONVCalibParam* gain = 
       static_cast<AliMUONVCalibParam*>(fGains->FindObject(detElemID,manuID));
     
     if (!gain) 
     {
-      gain = new AliMUONCalibParamNF(6,detElemID,manuID,kNchannels,0);
+      gain = new AliMUONCalibParamNF(5,kNchannels,detElemID,manuID,0);
       fGains->Add(gain);
     }
     gain->SetValueAsFloat(manuChannel,0,a0);
     gain->SetValueAsFloat(manuChannel,1,a1);
-    gain->SetValueAsFloat(manuChannel,2,a2);
-    gain->SetValueAsFloat(manuChannel,3,xlim);
-    gain->SetValueAsFloat(manuChannel,4,chi2);
-    gain->SetValueAsFloat(manuChannel,5,chi22);
+    gain->SetValueAsInt(manuChannel,2,thres);
+    gain->SetValueAsInt(manuChannel,3,qual);
+    gain->SetValueAsInt(manuChannel,4,kSaturation);
     ++n;
   }
   in.close();