CMS 3D CMS Logo

Public Member Functions | Private Member Functions | Private Attributes

CalibrationAlgorithm Class Reference

Algorithm for calibration runs. More...

#include <CalibrationAlgorithm.h>

Inheritance diagram for CalibrationAlgorithm:
CommissioningAlgorithm

List of all members.

Public Member Functions

 CalibrationAlgorithm (const edm::ParameterSet &pset, CalibrationAnalysis *const )
const Histohisto (int i) const
virtual ~CalibrationAlgorithm ()

Private Member Functions

void analyse ()
 CalibrationAlgorithm ()
void correctDistribution (TH1 *) const
virtual void extract (const std::vector< TH1 * > &)
TF1 * fitPulse (TH1 *, float rangeLow=0, float rangeHigh=-1)
float maximum (TH1 *)
float turnOn (TH1 *)

Private Attributes

CalibrationAnalysiscal_
TF1 * deconv_fitter_
Histo histo_ [32]
TF1 * peak_fitter_

Detailed Description

Algorithm for calibration runs.

Author:
C. Delaere

Definition at line 18 of file CalibrationAlgorithm.h.


Constructor & Destructor Documentation

CalibrationAlgorithm::CalibrationAlgorithm ( const edm::ParameterSet pset,
CalibrationAnalysis * const  anal 
)

Definition at line 20 of file CalibrationAlgorithm.cc.

References deconv_fitter_, fdeconv_convoluted(), fpeak_convoluted(), and peak_fitter_.

  : CommissioningAlgorithm(anal),
    deconv_fitter_(0),
    peak_fitter_(0),
    cal_(0)
{
  deconv_fitter_ = new TF1("deconv_fitter",fdeconv_convoluted,-50,50,5);
  deconv_fitter_->FixParameter(0,0);
  deconv_fitter_->SetParLimits(1,-100,0);
  deconv_fitter_->SetParLimits(2,0,200);
  deconv_fitter_->SetParLimits(3,5,100);
  deconv_fitter_->FixParameter(3,50);
  deconv_fitter_->SetParLimits(4,0,50);
  deconv_fitter_->SetParameters(0.,-10,0.96,50,20);
  peak_fitter_ = new TF1("peak_fitter",fpeak_convoluted,-50,50,5);
  peak_fitter_->FixParameter(0,0);
  peak_fitter_->SetParLimits(1,-100,0);
  peak_fitter_->SetParLimits(2,0,400);
  peak_fitter_->SetParLimits(3,5,100);
  peak_fitter_->FixParameter(3,50);
  peak_fitter_->SetParLimits(4,0,50);
  peak_fitter_->SetParameters(0.,-10,0.96,50,20);
}
virtual CalibrationAlgorithm::~CalibrationAlgorithm ( ) [inline, virtual]

Definition at line 24 of file CalibrationAlgorithm.h.

{;}
CalibrationAlgorithm::CalibrationAlgorithm ( ) [inline, private]

Definition at line 30 of file CalibrationAlgorithm.h.

{;}

Member Function Documentation

void CalibrationAlgorithm::analyse ( ) [private, virtual]

Performs histogram anaylsis.

Implements CommissioningAlgorithm.

Definition at line 101 of file CalibrationAlgorithm.cc.

References CalibrationAnalysis::amplitude_, cal_, CalibrationAnalysis::calchan_, CalibrationAnalysis::chi2_, correctDistribution(), first, fitPulse(), histo_, i, CalibrationAnalysis::isScan_, CalibrationAnalysis::max_amplitude_, CalibrationAnalysis::max_chi2_, CalibrationAnalysis::max_riseTime_, CalibrationAnalysis::max_smearing_, CalibrationAnalysis::max_tail_, CalibrationAnalysis::max_timeConstant_, maximum(), CalibrationAnalysis::mean_amplitude_, CalibrationAnalysis::mean_chi2_, CalibrationAnalysis::mean_riseTime_, CalibrationAnalysis::mean_smearing_, CalibrationAnalysis::mean_tail_, CalibrationAnalysis::mean_timeConstant_, CalibrationAnalysis::min_amplitude_, CalibrationAnalysis::min_chi2_, CalibrationAnalysis::min_riseTime_, CalibrationAnalysis::min_smearing_, CalibrationAnalysis::min_tail_, CalibrationAnalysis::min_timeConstant_, sistrip::mlCommissioning_, CalibrationAnalysis::riseTime_, Gflash::Rmax, Gflash::Rmin, CalibrationAnalysis::smearing_, CalibrationAnalysis::spread_amplitude_, CalibrationAnalysis::spread_chi2_, CalibrationAnalysis::spread_riseTime_, CalibrationAnalysis::spread_smearing_, CalibrationAnalysis::spread_tail_, CalibrationAnalysis::spread_timeConstant_, mathSSE::sqrt(), strip(), CalibrationAnalysis::tail_, CalibrationAnalysis::timeConstant_, indexGen::title, Tmax, builder_last_value_cfg::Tmin, and turnOn().

                                   { 

  if ( !cal_ ) {
    edm::LogWarning(mlCommissioning_)
      << "[CalibrationAlgorithm::" << __func__ << "]"
      << " NULL pointer to derived Analysis object!";
    return; 
  }

  float Amean[2]   = {0.,0.};
  float Amin[2]    = {2000.,2000.};
  float Amax[2]    = {0.,0.};
  float Aspread[2] = {0.,0.};
  float Tmean[2]   = {0.,0.};
  float Tmin[2]    = {2000.,2000.};
  float Tmax[2]    = {0.,0.};
  float Tspread[2] = {0.,0.};
  float Rmean[2]   = {0.,0.};
  float Rmin[2]    = {2000.,2000.};
  float Rmax[2]    = {0.,0.};
  float Rspread[2] = {0.,0.};
  float Cmean[2]   = {0.,0.};
  float Cmin[2]    = {2000.,2000.};
  float Cmax[2]    = {0.,0.};
  float Cspread[2] = {0.,0.};
  float Smean[2]   = {0.,0.};
  float Smin[2]    = {2000.,2000.};
  float Smax[2]    = {0.,0.};
  float Sspread[2] = {0.,0.};
  float Kmean[2]   = {0.,0.};
  float Kmin[2]    = {2000000.,2000000.};
  float Kmax[2]    = {0.,0.};
  float Kspread[2] = {0.,0.};
 
  unsigned int upperLimit = cal_->isScan_ ? 2 : 32;
  float nStrips = cal_->isScan_ ? 1. : 16.;
  for(unsigned int i=0;i<upperLimit;++i) {
    if ( !histo_[i].first ) {
      edm::LogWarning(mlCommissioning_) 
        << " NULL pointer to histogram " << i << "!";
      return;
    }
     
    // determine which APV and strip is being looked at.
    std::string title = histo_[i].first->GetName();
    int apv = 0;
    int strip = 0;
    if(title.find("STRIP")!=std::string::npos && title.find("Apv")!=std::string::npos) {
      strip = atoi(title.c_str()+title.find("STRIP")+6);
      apv = (atoi(title.c_str()+title.find("Apv")+3))%2;
    } else {
      strip = (atoi(title.c_str()+title.find_last_of("_")+1))%16;
      apv = (atoi(title.c_str()+title.find_last_of("_")+1))/16;
      if(title.find("Apv")!=std::string::npos) {
        apv = (atoi(title.c_str()+title.find("Apv")+3))%2;
        strip = strip*8 + cal_->calchan_;
      } else {
        edm::LogWarning(mlCommissioning_) 
          << " Malformed histogram title! Strip/APV not retreived: " 
          << title;
      }
    }
    
    // rescale the plot
    correctDistribution(histo_[i].first);
    
    // amplitude
    cal_->amplitude_[apv][strip] = histo_[i].first->GetMaximum();
    
    // rise time
    cal_->riseTime_[apv][strip] = maximum(histo_[i].first) - turnOn(histo_[i].first);
    
    // tail 125 ns after the maximum
    int lastBin = histo_[i].first->FindBin(histo_[i].first->GetBinCenter(histo_[i].first->GetMaximumBin())+125);
    if(lastBin>histo_[i].first->GetNbinsX()-4) lastBin = histo_[i].first->GetNbinsX()-4;
     if(histo_[i].first->GetMaximum()!=0)
       cal_->tail_[apv][strip] = 100*histo_[i].first->GetBinContent(lastBin)/histo_[i].first->GetMaximum();
     else
       cal_->tail_[apv][strip] = 100;
  
    // perform the fit for the next quantities
    TF1* fit = fitPulse(histo_[i].first);
    
    // time constant
    cal_->timeConstant_[apv][strip] = fit->GetParameter(3);
  
    // smearing
    cal_->smearing_[apv][strip] = fit->GetParameter(4);
  
    // chi2
    cal_->chi2_[apv][strip] = fit->GetChisquare();
    
    //compute mean, max, min, spread
    Amean[apv] += cal_->amplitude_[apv][strip]/nStrips;
    Amin[apv] = Amin[apv]<cal_->amplitude_[apv][strip] ? Amin[apv] : cal_->amplitude_[apv][strip];
    Amax[apv] = Amax[apv]>cal_->amplitude_[apv][strip] ? Amax[apv] : cal_->amplitude_[apv][strip];
    Aspread[apv] += cal_->amplitude_[apv][strip]*cal_->amplitude_[apv][strip]/nStrips;
    Tmean[apv] += cal_->tail_[apv][strip]/nStrips;
    Tmin[apv] = Tmin[apv]<cal_->tail_[apv][strip] ? Tmin[apv] : cal_->tail_[apv][strip];
    Tmax[apv] = Tmax[apv]>cal_->tail_[apv][strip] ? Tmax[apv] : cal_->tail_[apv][strip];
    Tspread[apv] += cal_->tail_[apv][strip]*cal_->tail_[apv][strip]/nStrips;
    Rmean[apv] += cal_->riseTime_[apv][strip]/nStrips;
    Rmin[apv] = Rmin[apv]<cal_->riseTime_[apv][strip] ? Rmin[apv] : cal_->riseTime_[apv][strip];
    Rmax[apv] = Rmax[apv]>cal_->riseTime_[apv][strip] ? Rmax[apv] : cal_->riseTime_[apv][strip];
    Rspread[apv] += cal_->riseTime_[apv][strip]*cal_->riseTime_[apv][strip]/nStrips;
    Cmean[apv] += cal_->timeConstant_[apv][strip]/nStrips;
    Cmin[apv] = Cmin[apv]<cal_->timeConstant_[apv][strip] ? Cmin[apv] : cal_->timeConstant_[apv][strip];
    Cmax[apv] = Cmax[apv]>cal_->timeConstant_[apv][strip] ? Cmax[apv] : cal_->timeConstant_[apv][strip];
    Cspread[apv] += cal_->timeConstant_[apv][strip]*cal_->timeConstant_[apv][strip]/nStrips;
    Smean[apv] += cal_->smearing_[apv][strip]/nStrips;
    Smin[apv] = Smin[apv]<cal_->smearing_[apv][strip] ? Smin[apv] : cal_->smearing_[apv][strip];
    Smax[apv] = Smax[apv]>cal_->smearing_[apv][strip] ? Smax[apv] : cal_->smearing_[apv][strip];
    Sspread[apv] += cal_->smearing_[apv][strip]*cal_->smearing_[apv][strip]/nStrips;
    Kmean[apv] += cal_->chi2_[apv][strip]/nStrips;
    Kmin[apv] = Kmin[apv]<cal_->chi2_[apv][strip] ? Kmin[apv] : cal_->chi2_[apv][strip];
    Kmax[apv] = Kmax[apv]>cal_->chi2_[apv][strip] ? Kmax[apv] : cal_->chi2_[apv][strip];
    Kspread[apv] += cal_->chi2_[apv][strip]*cal_->chi2_[apv][strip]/nStrips;
  }
                                
  // fill the mean, max, min, spread, ... histograms.
  for(int i=0;i<2;++i) {
    cal_->mean_amplitude_[i] = Amean[i];
    cal_->mean_tail_[i] = Tmean[i];
    cal_->mean_riseTime_[i] = Rmean[i];
    cal_->mean_timeConstant_[i] = Cmean[i];
    cal_->mean_smearing_[i] = Smean[i];
    cal_->mean_chi2_[i] = Kmean[i];
    cal_->min_amplitude_[i] = Amin[i];
    cal_->min_tail_[i] = Tmin[i];
    cal_->min_riseTime_[i] = Rmin[i];
    cal_->min_timeConstant_[i] = Cmin[i];
    cal_->min_smearing_[i] = Smin[i];
    cal_->min_chi2_[i] = Kmin[i];
    cal_->max_amplitude_[i] = Amax[i];
    cal_->max_tail_[i] = Tmax[i];
    cal_->max_riseTime_[i] = Rmax[i];
    cal_->max_timeConstant_[i] = Cmax[i];
    cal_->max_smearing_[i] = Smax[i];
    cal_->max_chi2_[i] = Kmax[i];
    cal_->spread_amplitude_[i] = sqrt(fabs(Aspread[i]-Amean[i]*Amean[i]));
    cal_->spread_tail_[i] = sqrt(fabs(Tspread[i]-Tmean[i]*Tmean[i]));
    cal_->spread_riseTime_[i] = sqrt(fabs(Rspread[i]-Rmean[i]*Rmean[i]));
    cal_->spread_timeConstant_[i] = sqrt(fabs(Cspread[i]-Cmean[i]*Cmean[i]));
    cal_->spread_smearing_[i] = sqrt(fabs(Sspread[i]-Smean[i]*Smean[i]));
    cal_->spread_chi2_[i] = sqrt(fabs(Kspread[i]-Kmean[i]*Kmean[i]));
  }
}
void CalibrationAlgorithm::correctDistribution ( TH1 *  histo) const [private]

Definition at line 251 of file CalibrationAlgorithm.cc.

References cal_, and CalibrationAnalysis::isScan_.

Referenced by analyse().

{
  // return the curve
  histo->Scale(-1);
  if ( cal_ ) { if( cal_->isScan_ ) histo->Scale(1/16.); }
}
void CalibrationAlgorithm::extract ( const std::vector< TH1 * > &  ) [private, virtual]

Extracts and organises histograms.

Implements CommissioningAlgorithm.

Definition at line 46 of file CalibrationAlgorithm.cc.

References CommissioningAnalysis::addErrorCode(), CommissioningAlgorithm::anal(), cal_, sistrip::CALIBRATION, sistrip::CALIBRATION_DECO, sistrip::CALIBRATION_SCAN, sistrip::CALIBRATION_SCAN_DECO, CommissioningAlgorithm::extractFedKey(), CommissioningAnalysis::fedKey(), histo_, CalibrationAnalysis::isScan_, sistrip::mlCommissioning_, sistrip::numberOfHistos_, indexGen::title, tmp, and sistrip::unexpectedTask_.

                                                                 {
  if ( !anal() ) {
    edm::LogWarning(mlCommissioning_)
      << "[CalibrationAlgorithm::" << __func__ << "]"
      << " NULL pointer to base Analysis object!";
    return; 
  }

  CommissioningAnalysis* tmp = const_cast<CommissioningAnalysis*>( anal() );
  cal_ = dynamic_cast<CalibrationAnalysis*>( tmp );
  if ( !cal_ ) {
    edm::LogWarning(mlCommissioning_)
      << "[CalibrationAlgorithm::" << __func__ << "]"
      << " NULL pointer to derived Analysis object!";
    return; 
  }

  // Check number of histograms
  if ( histos.size() != 32 && histos.size() !=2 ) {
    cal_->addErrorCode(sistrip::numberOfHistos_);
  }
  
  // Extract FED key from histo title
  if ( !histos.empty() ) { cal_->fedKey( extractFedKey( histos.front() ) ); }
  
  // Extract histograms
  std::vector<TH1*>::const_iterator ihis = histos.begin();
  unsigned int cnt = 0;
  for ( ; ihis != histos.end(); ihis++,cnt++ ) {
    
    // Check for NULL pointer
    if ( !(*ihis) ) { continue; }
    
    // Check name
    SiStripHistoTitle title( (*ihis)->GetName() );
    if ( title.runType() != sistrip::CALIBRATION && 
         title.runType() != sistrip::CALIBRATION_DECO &&
         title.runType() != sistrip::CALIBRATION_SCAN && 
         title.runType() != sistrip::CALIBRATION_SCAN_DECO ) {
      cal_->addErrorCode(sistrip::unexpectedTask_);
      continue;
    }
    
    cal_->isScan_ = ( title.runType() == sistrip::CALIBRATION_SCAN || 
                      title.runType() == sistrip::CALIBRATION_SCAN_DECO ); 
    
    // Extract calibration histo
    histo_[cnt].first = *ihis;
    histo_[cnt].second = (*ihis)->GetTitle();
  }
  
}
TF1 * CalibrationAlgorithm::fitPulse ( TH1 *  histo,
float  rangeLow = 0,
float  rangeHigh = -1 
) [private]

Definition at line 260 of file CalibrationAlgorithm.cc.

References cal_, CalibrationAnalysis::deconv_, deconv_fitter_, error, i, N, peak_fitter_, and mathSSE::sqrt().

Referenced by analyse(), and maximum().

{
  if(!cal_) return 0;
  if(!histo) return 0;
  float noise = 4.;
  float N = round(histo->GetMaximum()/125.);
  float error = sqrt(2*N)*noise;
  //float error = sqrt(2)*noise*N; // ?????????
  for(int i=1;i<=histo->GetNbinsX();++i) {
    histo->SetBinError(i,error);
  }
  if (cal_->deconv_) {
    if(rangeLow>rangeHigh)
      histo->Fit(deconv_fitter_,"Q");
    else 
      histo->Fit(deconv_fitter_,"0Q","",rangeLow,rangeHigh);
    return deconv_fitter_; 
  } else {
    if(rangeLow>rangeHigh)
      histo->Fit(peak_fitter_,"Q");
    else 
      histo->Fit(peak_fitter_,"0Q","",rangeLow,rangeHigh);
    return peak_fitter_; 
  } 
}
const Histo& CalibrationAlgorithm::histo ( int  i) const [inline]

Definition at line 26 of file CalibrationAlgorithm.h.

References histo_, and i.

{ return histo_[i]; }
float CalibrationAlgorithm::maximum ( TH1 *  h) [private]

Definition at line 290 of file CalibrationAlgorithm.cc.

References newFWLiteAna::bin, and fitPulse().

Referenced by analyse().

                                            {
    int bin = h->GetMaximumBin();
    // fit around the maximum with the detector response and take the max from the fit
    TF1* fit = fitPulse(h,h->GetBinCenter(bin)-25,h->GetBinCenter(bin)+25);
    return fit->GetMaximumX();
}
float CalibrationAlgorithm::turnOn ( TH1 *  h) [private]

Definition at line 299 of file CalibrationAlgorithm.cc.

References newFWLiteAna::base, newFWLiteAna::bin, end, and cond::rpcobgas::time.

Referenced by analyse().

                                           {
    // localize the rising edge
    int bin=1;
    float amplitude = h->GetMaximum();
    for(; bin<= h->GetNbinsX() && h->GetBinContent(bin)<0.4*amplitude; ++bin) {}
    float end = h->GetBinLowEdge(bin);
    // fit the rising edge with a sigmoid
    TF1* sigmoid = new TF1("sigmoid","[0]/(1+exp(-[1]*(x-[2])))+[3]",0,end);
    sigmoid->SetParLimits(0,amplitude/10.,amplitude);
    sigmoid->SetParLimits(1,0.05,0.5);
    sigmoid->SetParLimits(2,end-10,end+10);
    sigmoid->SetParLimits(3,-amplitude/10.,amplitude/10.);
    sigmoid->SetParameters(amplitude/2.,0.1,end,0.);
    h->Fit(sigmoid,"0QR");
    // return the point where the fit = 3% signal.
    float time = 0.;
    float base = sigmoid->GetMinimum(0,end);
    for(;time<end && (sigmoid->Eval(time)-base)<0.03*(amplitude-base);time += 0.1) {}
    delete sigmoid;
    return time-0.05;
}

Member Data Documentation

Definition at line 57 of file CalibrationAlgorithm.h.

Referenced by analyse(), correctDistribution(), extract(), and fitPulse().

Fitter in deconvolution mode

Definition at line 52 of file CalibrationAlgorithm.h.

Referenced by CalibrationAlgorithm(), and fitPulse().

pulse shape

Definition at line 49 of file CalibrationAlgorithm.h.

Referenced by analyse(), extract(), and histo().

Fitter in peak mode

Definition at line 55 of file CalibrationAlgorithm.h.

Referenced by CalibrationAlgorithm(), and fitPulse().