SiStripCommissioningAnalysis/src/CalibrationAlgorithm.cc

0001 #include "DQM/SiStripCommissioningAnalysis/interface/CalibrationAlgorithm.h"
0002 #include "CondFormats/SiStripObjects/interface/CalibrationAnalysis.h"
0003 #include "DataFormats/SiStripCommon/interface/SiStripHistoTitle.h"
0004 #include "DataFormats/SiStripCommon/interface/SiStripEnumsAndStrings.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "DQM/SiStripCommissioningAnalysis/interface/SiStripPulseShape.h"
0007 #include "TProfile.h"
0008 #include "TF1.h"
0009 #include "TH1.h"
0010 #include "TVirtualFitter.h"
0011 #include "TFitResultPtr.h"
0012 #include "TFitResult.h"
0013 #include <iostream>
0014 #include <sstream>
0015 #include <iomanip>
0016 #include <cmath>
0017 #include "Math/MinimizerOptions.h"
0018
0019 using namespace sistrip;
0020
0021 // ----------------------------------------------------------------------------
0022 //
0023 CalibrationAlgorithm::CalibrationAlgorithm(const edm::ParameterSet& pset, CalibrationAnalysis* const anal)
0024     : CommissioningAlgorithm(anal), cal_(nullptr) {}
0025
0026 // ----------------------------------------------------------------------------
0027 //
0028 void CalibrationAlgorithm::extract(const std::vector<TH1*>& histos) {
0029   // extract analysis object which should be already created
0030   if (!anal()) {
0031     edm::LogWarning(mlCommissioning_) << "[CalibrationAlgorithm::" << __func__ << "]"
0032                                       << " NULL pointer to base Analysis object!";
0033     return;
0034   }
0035
0036   CommissioningAnalysis* tmp = const_cast<CommissioningAnalysis*>(anal());
0037   cal_ = dynamic_cast<CalibrationAnalysis*>(tmp);
0038
0039   if (!cal_) {
0040     edm::LogWarning(mlCommissioning_) << "[CalibrationAlgorithm::" << __func__ << "]"
0041                                       << " NULL pointer to derived Analysis object!";
0042     return;
0043   }
0044
0045   // Extract FED key from histo title
0046   if (!histos.empty()) {
0047     cal_->fedKey(extractFedKey(histos.front()));
0048   }
0049
0050   // Extract histograms
0051   std::vector<TH1*>::const_iterator ihis = histos.begin();
0052   unsigned int cnt = 0;
0053   for (; ihis != histos.end(); ihis++, cnt++) {
0054     // Check for NULL pointer
0055     if (!(*ihis)) {
0056       continue;
0057     }
0058
0059     // Check name
0060     SiStripHistoTitle title((*ihis)->GetName());
0061     if (title.runType() != sistrip::CALIBRATION && title.runType() != sistrip::CALIBRATION_DECO) {
0062       cal_->addErrorCode(sistrip::unexpectedTask_);
0063       continue;
0064     }
0065
0066     /// extract isha, vfs and calchan values, as well as filling the histogram objects
0067     std::vector<std::string> tokens;
0068     std::string token;
0069     std::istringstream tokenStream(title.extraInfo());
0070     while (std::getline(tokenStream, token, '_')) {
0071       tokens.push_back(token);
0072     }
0073
0074     ////////
0075     Histo histo_temp;
0076     histo_temp.first = *ihis;
0077     histo_temp.second = (*ihis)->GetTitle();
0078     histo_temp.first->Sumw2();
0079     histo_.push_back(histo_temp);
0080     apvId_.push_back(title.channel() % 2);
0081     stripId_.push_back(std::stoi(tokens.at(1)) * 16 + std::stoi(tokens.at(3)));
0082     calChan_.push_back(std::stoi(tokens.at(1)));
0083   }
0084 }
0085
0086 // ----------------------------------------------------------------------------
0087 //
0088 void CalibrationAlgorithm::analyse() {
0089   ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2", "Migrad");
0090   ROOT::Math::MinimizerOptions::SetDefaultStrategy(0);
0091
0092   if (!cal_) {
0093     edm::LogWarning(mlCommissioning_) << "[CalibrationAlgorithm::" << __func__ << "]"
0094                                       << " NULL pointer to derived Analysis object!";
0095     return;
0096   }
0097
0098   float Amean[2] = {-1., -1.};
0099   float Amin[2] = {-1., -1.};
0100   float Amax[2] = {-1., -1.};
0101   float Aspread[2] = {-1., -1.};
0102   float Tmean[2] = {-1., -1.};
0103   float Tmin[2] = {-1., -1.};
0104   float Tmax[2] = {-1., -1.};
0105   float Tspread[2] = {-1., -1.};
0106   float Rmean[2] = {-1., -1.};
0107   float Rmin[2] = {-1., -1.};
0108   float Rmax[2] = {-1., -1.};
0109   float Rspread[2] = {-1., -1.};
0110   float Cmean[2] = {-1., -1.};
0111   float Cmin[2] = {-1., -1.};
0112   float Cmax[2] = {-1., -1.};
0113   float Cspread[2] = {-1., -1.};
0114   float Smean[2] = {-1., -1.};
0115   float Smin[2] = {-1., -1.};
0116   float Smax[2] = {-1., -1.};
0117   float Sspread[2] = {-1., -1.};
0118   float Kmean[2] = {-1., -1.};
0119   float Kmin[2] = {-1., -1.};
0120   float Kmax[2] = {-1., -1.};
0121   float Kspread[2] = {-1., -1.};
0122   // turnOn
0123   float Omean[2] = {-1., -1.};
0124   float Omin[2] = {-1., -1.};
0125   float Omax[2] = {-1., -1.};
0126   float Ospread[2] = {-1., -1.};
0127   // maximum
0128   float Mmean[2] = {-1., -1.};
0129   float Mmin[2] = {-1., -1.};
0130   float Mmax[2] = {-1., -1.};
0131   float Mspread[2] = {-1., -1.};
0132   // undershoot
0133   float Umean[2] = {-1., -1.};
0134   float Umin[2] = {-1., -1.};
0135   float Umax[2] = {-1., -1.};
0136   float Uspread[2] = {-1., -1.};
0137   // baseline
0138   float Bmean[2] = {-1., -1.};
0139   float Bmin[2] = {-1., -1.};
0140   float Bmax[2] = {-1., -1.};
0141   float Bspread[2] = {-1., -1.};
0142
0143   ////////
0144   TFitResultPtr fit_result;
0145   TF1* fit_function = nullptr;
0146   if (cal_->deconv_) {
0147     fit_function = new TF1("fit_function_deco", fdeconv, 0, 400, 7);
0148     fit_function->SetParameters(4, 25, 25, 50, 250, 25, 0.75);
0149   } else {
0150     fit_function = new TF1("fit_function_peak", fpeak, 0, 400, 6);
0151     fit_function->SetParameters(4, 50, 50, 70, 250, 20);
0152   }
0153
0154   //////////
0155   std::vector<unsigned int> nStrips(2, 0.);
0156
0157   for (size_t ihist = 0; ihist < histo_.size(); ihist++) {
0158     if (!histo_[ihist].first) {
0159       edm::LogWarning(mlCommissioning_) << " NULL pointer to histogram for: " << histo_[ihist].second << " !";
0160       return;
0161     }
0162
0163     cal_->amplitude_[apvId_[ihist]][stripId_[ihist]] = 0;
0164     cal_->baseline_[apvId_[ihist]][stripId_[ihist]] = 0;
0165     cal_->riseTime_[apvId_[ihist]][stripId_[ihist]] = 0;
0166     cal_->turnOn_[apvId_[ihist]][stripId_[ihist]] = 0;
0167     cal_->peakTime_[apvId_[ihist]][stripId_[ihist]] = 0;
0168     cal_->undershoot_[apvId_[ihist]][stripId_[ihist]] = 0;
0169     cal_->tail_[apvId_[ihist]][stripId_[ihist]] = 0;
0170     cal_->decayTime_[apvId_[ihist]][stripId_[ihist]] = 0;
0171     cal_->smearing_[apvId_[ihist]][stripId_[ihist]] = 0;
0172     cal_->chi2_[apvId_[ihist]][stripId_[ihist]] = 0;
0173     cal_->isvalid_[apvId_[ihist]][stripId_[ihist]] = true;
0174
0175     if (histo_[ihist].first->Integral() == 0) {
0176       cal_->isvalid_[apvId_[ihist]][stripId_[ihist]] = false;
0177       continue;
0178     }
0179
0180     // rescale the plot and set reasonable errors
0181     correctDistribution(histo_[ihist].first);
0182
0183     // from NOTE2009_021 : The charge injection provided by the calibration circuit is known with a precision of 5%
0184     float error = histo_[ihist].first->GetMaximum() * 0.05;
0185     for (int i = 1; i <= histo_[ihist].first->GetNbinsX(); ++i)
0186       histo_[ihist].first->SetBinError(i, error);
0187
0188     // set intial par
0189     if (cal_->deconv_)
0190       fit_function->SetParameters(10, 15, 30, 10, 350, 50, 0.75);
0191     else
0192       fit_function->SetParameters(6, 40, 40, 70, 350, 20);
0193
0194     fit_result = histo_[ihist].first->Fit(fit_function, "QS");
0195
0196     // fit-result should exist and have a resonably good status
0197     if (not fit_result.Get())
0198       continue;
0199
0200     float maximum_ampl = fit_function->GetMaximum();
0201     float peak_time = fit_function->GetMaximumX();
0202     float baseline = baseLine(fit_function);
0203     float turn_on_time = turnOn(fit_function, baseline);
0204     float rise_time = peak_time - turn_on_time;
0205
0206     // start filling info
0207     cal_->amplitude_[apvId_[ihist]][stripId_[ihist]] = maximum_ampl - baseline;
0208     cal_->baseline_[apvId_[ihist]][stripId_[ihist]] = baseline;
0209     cal_->riseTime_[apvId_[ihist]][stripId_[ihist]] = rise_time;
0210     cal_->turnOn_[apvId_[ihist]][stripId_[ihist]] = turn_on_time;
0211     cal_->peakTime_[apvId_[ihist]][stripId_[ihist]] = peak_time;
0212
0213     if (cal_->deconv_ and fit_function->GetMinimumX() > peak_time)  // make sure the minimum is after the peak-time
0214       cal_->undershoot_[apvId_[ihist]][stripId_[ihist]] =
0215           100 * (fit_function->GetMinimum() - baseline) / (maximum_ampl - baseline);
0216     else
0217       cal_->undershoot_[apvId_[ihist]][stripId_[ihist]] = 0;
0218
0219     // Bin related to peak + 125 ns
0220     int lastBin = histo_[ihist].first->FindBin(peak_time + 125);
0221     if (lastBin > histo_[ihist].first->GetNbinsX() - 4)
0222       lastBin = histo_[ihist].first->GetNbinsX() - 4;
0223
0224     // tail is the amplitude at 5 bx from the maximum
0225     cal_->tail_[apvId_[ihist]][stripId_[ihist]] =
0226         100 * (histo_[ihist].first->GetBinContent(lastBin) - baseline) / (maximum_ampl - baseline);
0227
0228     // reaches 1/e of the peak amplitude
0229     cal_->decayTime_[apvId_[ihist]][stripId_[ihist]] = decayTime(fit_function) - peak_time;
0230     cal_->smearing_[apvId_[ihist]][stripId_[ihist]] = 0;
0231     cal_->chi2_[apvId_[ihist]][stripId_[ihist]] =
0232         fit_function->GetChisquare() / (histo_[ihist].first->GetNbinsX() - fit_function->GetNpar());
0233
0234     // calibration channel
0235     cal_->calChan_ = calChan_[ihist];
0236
0237     // apply quality requirements
0238     bool isvalid = true;
0239     if (not cal_->deconv_) {  // peak-mode
0240
0241       if (cal_->amplitude_[apvId_[ihist]][stripId_[ihist]] < CalibrationAnalysis::minAmplitudeThreshold_)
0242         isvalid = false;
0243       if (cal_->baseline_[apvId_[ihist]][stripId_[ihist]] < CalibrationAnalysis::minBaselineThreshold_)
0244         isvalid = false;
0245       else if (cal_->baseline_[apvId_[ihist]][stripId_[ihist]] > CalibrationAnalysis::maxBaselineThreshold_)
0246         isvalid = false;
0247       if (cal_->decayTime_[apvId_[ihist]][stripId_[ihist]] < CalibrationAnalysis::minDecayTimeThreshold_)
0248         isvalid = false;
0249       else if (cal_->decayTime_[apvId_[ihist]][stripId_[ihist]] > CalibrationAnalysis::maxDecayTimeThreshold_)
0250         isvalid = false;
0251       if (cal_->peakTime_[apvId_[ihist]][stripId_[ihist]] < CalibrationAnalysis::minPeakTimeThreshold_)
0252         isvalid = false;
0253       else if (cal_->peakTime_[apvId_[ihist]][stripId_[ihist]] > CalibrationAnalysis::maxPeakTimeThreshold_)
0254         isvalid = false;
0255       if (cal_->riseTime_[apvId_[ihist]][stripId_[ihist]] < CalibrationAnalysis::minRiseTimeThreshold_)
0256         isvalid = false;
0257       else if (cal_->riseTime_[apvId_[ihist]][stripId_[ihist]] > CalibrationAnalysis::maxRiseTimeThreshold_)
0258         isvalid = false;
0259       if (cal_->turnOn_[apvId_[ihist]][stripId_[ihist]] < CalibrationAnalysis::minTurnOnThreshold_)
0260         isvalid = false;
0261       else if (cal_->turnOn_[apvId_[ihist]][stripId_[ihist]] > CalibrationAnalysis::maxTurnOnThreshold_)
0262         isvalid = false;
0263       if (cal_->chi2_[apvId_[ihist]][stripId_[ihist]] > CalibrationAnalysis::maxChi2Threshold_)
0264         isvalid = false;
0265
0266     } else {
0267       if (fit_function->GetMinimumX() < peak_time)
0268         isvalid = false;
0269       if (cal_->amplitude_[apvId_[ihist]][stripId_[ihist]] < CalibrationAnalysis::minAmplitudeThreshold_)
0270         isvalid = false;
0271       if (cal_->baseline_[apvId_[ihist]][stripId_[ihist]] < CalibrationAnalysis::minBaselineThreshold_)
0272         isvalid = false;
0273       if (cal_->baseline_[apvId_[ihist]][stripId_[ihist]] < CalibrationAnalysis::minBaselineThreshold_)
0274         isvalid = false;
0275       else if (cal_->baseline_[apvId_[ihist]][stripId_[ihist]] > CalibrationAnalysis::maxBaselineThreshold_)
0276         isvalid = false;
0277       if (cal_->chi2_[apvId_[ihist]][stripId_[ihist]] > CalibrationAnalysis::maxChi2Threshold_)
0278         isvalid = false;
0279       if (cal_->turnOn_[apvId_[ihist]][stripId_[ihist]] < CalibrationAnalysis::minTurnOnThresholdDeco_)
0280         isvalid = false;
0281       else if (cal_->turnOn_[apvId_[ihist]][stripId_[ihist]] > CalibrationAnalysis::maxTurnOnThresholdDeco_)
0282         isvalid = false;
0283       if (cal_->decayTime_[apvId_[ihist]][stripId_[ihist]] < CalibrationAnalysis::minDecayTimeThresholdDeco_)
0284         isvalid = false;
0285       else if (cal_->decayTime_[apvId_[ihist]][stripId_[ihist]] > CalibrationAnalysis::maxDecayTimeThresholdDeco_)
0286         isvalid = false;
0287       if (cal_->peakTime_[apvId_[ihist]][stripId_[ihist]] < CalibrationAnalysis::minPeakTimeThresholdDeco_)
0288         isvalid = false;
0289       else if (cal_->peakTime_[apvId_[ihist]][stripId_[ihist]] > CalibrationAnalysis::maxPeakTimeThresholdDeco_)
0290         isvalid = false;
0291       if (cal_->riseTime_[apvId_[ihist]][stripId_[ihist]] < CalibrationAnalysis::minRiseTimeThresholdDeco_)
0292         isvalid = false;
0293       else if (cal_->riseTime_[apvId_[ihist]][stripId_[ihist]] > CalibrationAnalysis::maxRiseTimeThresholdDeco_)
0294         isvalid = false;
0295     }
0296
0297     if (not isvalid) {  // not valid set default to zero for all quantities
0298
0299       cal_->amplitude_[apvId_[ihist]][stripId_[ihist]] = 0;
0300       cal_->baseline_[apvId_[ihist]][stripId_[ihist]] = 0;
0301       cal_->riseTime_[apvId_[ihist]][stripId_[ihist]] = 0;
0302       cal_->turnOn_[apvId_[ihist]][stripId_[ihist]] = 0;
0303       cal_->peakTime_[apvId_[ihist]][stripId_[ihist]] = 0;
0304       cal_->undershoot_[apvId_[ihist]][stripId_[ihist]] = 0;
0305       cal_->tail_[apvId_[ihist]][stripId_[ihist]] = 0;
0306       cal_->decayTime_[apvId_[ihist]][stripId_[ihist]] = 0;
0307       cal_->smearing_[apvId_[ihist]][stripId_[ihist]] = 0;
0308       cal_->chi2_[apvId_[ihist]][stripId_[ihist]] = 0;
0309       cal_->isvalid_[apvId_[ihist]][stripId_[ihist]] = false;
0310       continue;
0311     }
0312
0313     // in case is valid
0314     nStrips[apvId_[ihist]]++;
0315
0316     //compute mean, max, min, spread only for valid strips
0317     Amean[apvId_[ihist]] += cal_->amplitude_[apvId_[ihist]][stripId_[ihist]];
0318     Amin[apvId_[ihist]] = Amin[apvId_[ihist]] < cal_->amplitude_[apvId_[ihist]][stripId_[ihist]]
0319                               ? Amin[apvId_[ihist]]
0320                               : cal_->amplitude_[apvId_[ihist]][stripId_[ihist]];
0321     Amax[apvId_[ihist]] = Amax[apvId_[ihist]] > cal_->amplitude_[apvId_[ihist]][stripId_[ihist]]
0322                               ? Amax[apvId_[ihist]]
0323                               : cal_->amplitude_[apvId_[ihist]][stripId_[ihist]];
0324     Aspread[apvId_[ihist]] +=
0325         cal_->amplitude_[apvId_[ihist]][stripId_[ihist]] * cal_->amplitude_[apvId_[ihist]][stripId_[ihist]];
0326
0327     Tmean[apvId_[ihist]] += cal_->tail_[apvId_[ihist]][stripId_[ihist]];
0328     Tmin[apvId_[ihist]] = Tmin[apvId_[ihist]] < cal_->tail_[apvId_[ihist]][stripId_[ihist]]
0329                               ? Tmin[apvId_[ihist]]
0330                               : cal_->tail_[apvId_[ihist]][stripId_[ihist]];
0331     Tmax[apvId_[ihist]] = Tmax[apvId_[ihist]] > cal_->tail_[apvId_[ihist]][stripId_[ihist]]
0332                               ? Tmax[apvId_[ihist]]
0333                               : cal_->tail_[apvId_[ihist]][stripId_[ihist]];
0334     Tspread[apvId_[ihist]] += cal_->tail_[apvId_[ihist]][stripId_[ihist]] * cal_->tail_[apvId_[ihist]][stripId_[ihist]];
0335
0336     Rmean[apvId_[ihist]] += cal_->riseTime_[apvId_[ihist]][stripId_[ihist]];
0337     Rmin[apvId_[ihist]] = Rmin[apvId_[ihist]] < cal_->riseTime_[apvId_[ihist]][stripId_[ihist]]
0338                               ? Rmin[apvId_[ihist]]
0339                               : cal_->riseTime_[apvId_[ihist]][stripId_[ihist]];
0340     Rmax[apvId_[ihist]] = Rmax[apvId_[ihist]] > cal_->riseTime_[apvId_[ihist]][stripId_[ihist]]
0341                               ? Rmax[apvId_[ihist]]
0342                               : cal_->riseTime_[apvId_[ihist]][stripId_[ihist]];
0343     Rspread[apvId_[ihist]] +=
0344         cal_->riseTime_[apvId_[ihist]][stripId_[ihist]] * cal_->riseTime_[apvId_[ihist]][stripId_[ihist]];
0345
0346     Cmean[apvId_[ihist]] += cal_->decayTime_[apvId_[ihist]][stripId_[ihist]];
0347     Cmin[apvId_[ihist]] = Cmin[apvId_[ihist]] < cal_->decayTime_[apvId_[ihist]][stripId_[ihist]]
0348                               ? Cmin[apvId_[ihist]]
0349                               : cal_->decayTime_[apvId_[ihist]][stripId_[ihist]];
0350     Cmax[apvId_[ihist]] = Cmax[apvId_[ihist]] > cal_->decayTime_[apvId_[ihist]][stripId_[ihist]]
0351                               ? Cmax[apvId_[ihist]]
0352                               : cal_->decayTime_[apvId_[ihist]][stripId_[ihist]];
0353     Cspread[apvId_[ihist]] +=
0354         cal_->decayTime_[apvId_[ihist]][stripId_[ihist]] * cal_->decayTime_[apvId_[ihist]][stripId_[ihist]];
0355
0356     Smean[apvId_[ihist]] += cal_->smearing_[apvId_[ihist]][stripId_[ihist]];
0357     Smin[apvId_[ihist]] = Smin[apvId_[ihist]] < cal_->smearing_[apvId_[ihist]][stripId_[ihist]]
0358                               ? Smin[apvId_[ihist]]
0359                               : cal_->smearing_[apvId_[ihist]][stripId_[ihist]];
0360     Smax[apvId_[ihist]] = Smax[apvId_[ihist]] > cal_->smearing_[apvId_[ihist]][stripId_[ihist]]
0361                               ? Smax[apvId_[ihist]]
0362                               : cal_->smearing_[apvId_[ihist]][stripId_[ihist]];
0363     Sspread[apvId_[ihist]] +=
0364         cal_->smearing_[apvId_[ihist]][stripId_[ihist]] * cal_->smearing_[apvId_[ihist]][stripId_[ihist]];
0365
0366     Kmean[apvId_[ihist]] += cal_->chi2_[apvId_[ihist]][stripId_[ihist]];
0367     Kmin[apvId_[ihist]] = Kmin[apvId_[ihist]] < cal_->chi2_[apvId_[ihist]][stripId_[ihist]]
0368                               ? Kmin[apvId_[ihist]]
0369                               : cal_->chi2_[apvId_[ihist]][stripId_[ihist]];
0370     Kmax[apvId_[ihist]] = Kmax[apvId_[ihist]] > cal_->chi2_[apvId_[ihist]][stripId_[ihist]]
0371                               ? Kmax[apvId_[ihist]]
0372                               : cal_->chi2_[apvId_[ihist]][stripId_[ihist]];
0373     Kspread[apvId_[ihist]] += cal_->chi2_[apvId_[ihist]][stripId_[ihist]] * cal_->chi2_[apvId_[ihist]][stripId_[ihist]];
0374
0375     Omean[apvId_[ihist]] += cal_->turnOn_[apvId_[ihist]][stripId_[ihist]];
0376     Omin[apvId_[ihist]] = Omin[apvId_[ihist]] < cal_->turnOn_[apvId_[ihist]][stripId_[ihist]]
0377                               ? Omin[apvId_[ihist]]
0378                               : cal_->turnOn_[apvId_[ihist]][stripId_[ihist]];
0379     Omax[apvId_[ihist]] = Omax[apvId_[ihist]] > cal_->turnOn_[apvId_[ihist]][stripId_[ihist]]
0380                               ? Omax[apvId_[ihist]]
0381                               : cal_->turnOn_[apvId_[ihist]][stripId_[ihist]];
0382     Ospread[apvId_[ihist]] +=
0383         cal_->turnOn_[apvId_[ihist]][stripId_[ihist]] * cal_->turnOn_[apvId_[ihist]][stripId_[ihist]];
0384
0385     Mmean[apvId_[ihist]] += cal_->peakTime_[apvId_[ihist]][stripId_[ihist]];
0386     Mmin[apvId_[ihist]] = Mmin[apvId_[ihist]] < cal_->peakTime_[apvId_[ihist]][stripId_[ihist]]
0387                               ? Mmin[apvId_[ihist]]
0388                               : cal_->peakTime_[apvId_[ihist]][stripId_[ihist]];
0389     Mmax[apvId_[ihist]] = Mmax[apvId_[ihist]] > cal_->peakTime_[apvId_[ihist]][stripId_[ihist]]
0390                               ? Mmax[apvId_[ihist]]
0391                               : cal_->peakTime_[apvId_[ihist]][stripId_[ihist]];
0392     Mspread[apvId_[ihist]] +=
0393         cal_->peakTime_[apvId_[ihist]][stripId_[ihist]] * cal_->peakTime_[apvId_[ihist]][stripId_[ihist]];
0394
0395     Umean[apvId_[ihist]] += cal_->undershoot_[apvId_[ihist]][stripId_[ihist]];
0396     Umin[apvId_[ihist]] = Umin[apvId_[ihist]] < cal_->undershoot_[apvId_[ihist]][stripId_[ihist]]
0397                               ? Umin[apvId_[ihist]]
0398                               : cal_->undershoot_[apvId_[ihist]][stripId_[ihist]];
0399     Umax[apvId_[ihist]] = Umax[apvId_[ihist]] > cal_->undershoot_[apvId_[ihist]][stripId_[ihist]]
0400                               ? Umax[apvId_[ihist]]
0401                               : cal_->undershoot_[apvId_[ihist]][stripId_[ihist]];
0402     Uspread[apvId_[ihist]] +=
0403         cal_->undershoot_[apvId_[ihist]][stripId_[ihist]] * cal_->undershoot_[apvId_[ihist]][stripId_[ihist]];
0404
0405     Bmean[apvId_[ihist]] += cal_->baseline_[apvId_[ihist]][stripId_[ihist]];
0406     Bmin[apvId_[ihist]] = Bmin[apvId_[ihist]] < cal_->baseline_[apvId_[ihist]][stripId_[ihist]]
0407                               ? Bmin[apvId_[ihist]]
0408                               : cal_->baseline_[apvId_[ihist]][stripId_[ihist]];
0409     Bmax[apvId_[ihist]] = Bmax[apvId_[ihist]] > cal_->baseline_[apvId_[ihist]][stripId_[ihist]]
0410                               ? Bmax[apvId_[ihist]]
0411                               : cal_->baseline_[apvId_[ihist]][stripId_[ihist]];
0412     Bspread[apvId_[ihist]] +=
0413         cal_->baseline_[apvId_[ihist]][stripId_[ihist]] * cal_->baseline_[apvId_[ihist]][stripId_[ihist]];
0414   }
0415
0416   // make mean values
0417   for (int i = 0; i < 2; i++) {
0418     if (nStrips[i] != 0) {
0419       Amean[i] = Amean[i] / nStrips[i];
0420       Tmean[i] = Tmean[i] / nStrips[i];
0421       Rmean[i] = Rmean[i] / nStrips[i];
0422       Cmean[i] = Cmean[i] / nStrips[i];
0423       Omean[i] = Omean[i] / nStrips[i];
0424       Mmean[i] = Mmean[i] / nStrips[i];
0425       Umean[i] = Umean[i] / nStrips[i];
0426       Bmean[i] = Bmean[i] / nStrips[i];
0427       Smean[i] = Smean[i] / nStrips[i];
0428       Kmean[i] = Kmean[i] / nStrips[i];
0429
0430       Aspread[i] = Aspread[i] / nStrips[i];
0431       Tspread[i] = Tspread[i] / nStrips[i];
0432       Rspread[i] = Rspread[i] / nStrips[i];
0433       Cspread[i] = Cspread[i] / nStrips[i];
0434       Ospread[i] = Ospread[i] / nStrips[i];
0435       Mspread[i] = Mspread[i] / nStrips[i];
0436       Uspread[i] = Uspread[i] / nStrips[i];
0437       Bspread[i] = Bspread[i] / nStrips[i];
0438       Sspread[i] = Sspread[i] / nStrips[i];
0439       Kspread[i] = Kspread[i] / nStrips[i];
0440     }
0441   }
0442
0443   // fill the mean, max, min, spread, ... histograms.
0444   for (int i = 0; i < 2; ++i) {
0445     cal_->mean_amplitude_[i] = Amean[i];
0446     cal_->mean_tail_[i] = Tmean[i];
0447     cal_->mean_riseTime_[i] = Rmean[i];
0448     cal_->mean_decayTime_[i] = Cmean[i];
0449     cal_->mean_turnOn_[i] = Omean[i];
0450     cal_->mean_peakTime_[i] = Mmean[i];
0451     cal_->mean_undershoot_[i] = Umean[i];
0452     cal_->mean_baseline_[i] = Bmean[i];
0453     cal_->mean_smearing_[i] = Smean[i];
0454     cal_->mean_chi2_[i] = Kmean[i];
0455
0456     cal_->min_amplitude_[i] = Amin[i];
0457     cal_->min_tail_[i] = Tmin[i];
0458     cal_->min_riseTime_[i] = Rmin[i];
0459     cal_->min_decayTime_[i] = Cmin[i];
0460     cal_->min_turnOn_[i] = Omin[i];
0461     cal_->min_peakTime_[i] = Mmin[i];
0462     cal_->min_undershoot_[i] = Umin[i];
0463     cal_->min_baseline_[i] = Bmin[i];
0464     cal_->min_smearing_[i] = Smin[i];
0465     cal_->min_chi2_[i] = Kmin[i];
0466
0467     cal_->max_amplitude_[i] = Amax[i];
0468     cal_->max_tail_[i] = Tmax[i];
0469     cal_->max_riseTime_[i] = Rmax[i];
0470     cal_->max_decayTime_[i] = Cmax[i];
0471     cal_->max_turnOn_[i] = Omax[i];
0472     cal_->max_peakTime_[i] = Mmax[i];
0473     cal_->max_undershoot_[i] = Umax[i];
0474     cal_->max_baseline_[i] = Bmax[i];
0475     cal_->max_smearing_[i] = Smax[i];
0476     cal_->max_chi2_[i] = Kmax[i];
0477
0478     cal_->spread_amplitude_[i] = sqrt(fabs(Aspread[i] - Amean[i] * Amean[i]));
0479     cal_->spread_tail_[i] = sqrt(fabs(Tspread[i] - Tmean[i] * Tmean[i]));
0480     cal_->spread_riseTime_[i] = sqrt(fabs(Rspread[i] - Rmean[i] * Rmean[i]));
0481     cal_->spread_decayTime_[i] = sqrt(fabs(Cspread[i] - Cmean[i] * Cmean[i]));
0482     cal_->spread_turnOn_[i] = sqrt(fabs(Ospread[i] - Omean[i] * Omean[i]));
0483     cal_->spread_peakTime_[i] = sqrt(fabs(Mspread[i] - Mmean[i] * Mmean[i]));
0484     cal_->spread_undershoot_[i] = sqrt(fabs(Uspread[i] - Umean[i] * Umean[i]));
0485     cal_->spread_baseline_[i] = sqrt(fabs(Bspread[i] - Bmean[i] * Bmean[i]));
0486     cal_->spread_smearing_[i] = sqrt(fabs(Sspread[i] - Smean[i] * Smean[i]));
0487     cal_->spread_chi2_[i] = sqrt(fabs(Kspread[i] - Kmean[i] * Kmean[i]));
0488   }
0489
0490   if (fit_function)
0491     delete fit_function;
0492 }
0493
0494 // ------
0495 void CalibrationAlgorithm::correctDistribution(TH1* histo) const {
0496   // 20 events per point in the TM loop  --> divide by 20 to have the amplitude of a single event readout
0497   for (int iBin = 0; iBin < histo->GetNbinsX(); iBin++) {
0498     histo->SetBinContent(iBin + 1, -histo->GetBinContent(iBin + 1) / 20.);
0499   }
0500 }
0501
0502 // ----------------------------------------------------------------------------
0503 float CalibrationAlgorithm::baseLine(TF1* f) {
0504   float xmax = 10;
0505   float baseline = 0;
0506   int npoints = 0;
0507   float x = f->GetXmin();
0508   for (; x < xmax; x += 0.1) {
0509     baseline += f->Eval(x);
0510     npoints++;
0511   }
0512   return baseline / npoints;
0513 }
0514
0515 // ----------------------------------------------------------------------------
0516 float CalibrationAlgorithm::turnOn(TF1* f,
0517                                    const float& baseline) {  // should happen within 100 ns in both deco and peak modes
0518   float max_amplitude = f->GetMaximum();
0519   float time = 10.;
0520   for (; time < 100 && (f->Eval(time) - baseline) < 0.05 * (max_amplitude - baseline); time += 0.1) {
0521   }  // flucutation higher than 5% of the pulse height
0522   return time;
0523 }
0524
0525 // ----------------------------------------------------------------------------
0526 float CalibrationAlgorithm::decayTime(
0527     TF1* f) {  // if we approximate the decay to an exp(-t/tau), in one constant unit, the amplited is reduced by e^{-1}
0528   float xval = f->GetMaximumX();
0529   float max_amplitude = f->GetMaximum();
0530   float x = xval;
0531   for (; x < 1000; x = x + 0.1) {
0532     if (f->Eval(x) < max_amplitude * exp(-1))
0533       break;
0534   }
0535   return x;
0536 }