Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 10:56:20 

0001 #include "DQM/SiStripCommissioningAnalysis/interface/CalibrationScanAlgorithm.h"
0002 #include "CondFormats/SiStripObjects/interface/CalibrationScanAnalysis.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 "TFitResult.h"
0012 #include "TMath.h"
0013 #include "TGraph2D.h"
0014 #include "TH2F.h"
0015 #include "TCanvas.h"
0016 #include "TROOT.h"
0017 #include <iostream>
0018 #include <sstream>
0019 #include <iomanip>
0020 #include <cmath>
0021 #include "Math/MinimizerOptions.h"
0022 
0023 using namespace sistrip;
0024 
0025 // ----------------------------------------------------------------------------
0026 //
0027 CalibrationScanAlgorithm::CalibrationScanAlgorithm(const edm::ParameterSet& pset, CalibrationScanAnalysis* const anal)
0028     : CommissioningAlgorithm(anal), cal_(nullptr) {}
0029 
0030 // ----------------------------------------------------------------------------
0031 //
0032 void CalibrationScanAlgorithm::extract(const std::vector<TH1*>& histos) {
0033   if (!anal()) {
0034     edm::LogWarning(mlCommissioning_) << "[CalibrationScanAlgorithm::" << __func__ << "]"
0035                                       << " NULL pointer to base Analysis object!";
0036     return;
0037   }
0038 
0039   CommissioningAnalysis* tmp = const_cast<CommissioningAnalysis*>(anal());
0040   cal_ = dynamic_cast<CalibrationScanAnalysis*>(tmp);
0041   if (!cal_) {
0042     edm::LogWarning(mlCommissioning_) << "[CalibrationScanAlgorithm::" << __func__ << "]"
0043                                       << " NULL pointer to derived Analysis object!";
0044     return;
0045   }
0046 
0047   // Extract FED key from histo title
0048   if (!histos.empty()) {
0049     cal_->fedKey(extractFedKey(histos.front()));
0050   }
0051 
0052   // Extract histograms
0053   std::vector<TH1*>::const_iterator ihis = histos.begin();
0054   unsigned int cnt = 0;
0055   for (; ihis != histos.end(); ihis++, cnt++) {
0056     // Check for NULL pointer
0057     if (!(*ihis)) {
0058       continue;
0059     }
0060 
0061     // Check name
0062     SiStripHistoTitle title((*ihis)->GetName());
0063     if (title.runType() != sistrip::CALIBRATION_SCAN && title.runType() != sistrip::CALIBRATION_SCAN_DECO) {
0064       cal_->addErrorCode(sistrip::unexpectedTask_);
0065       continue;
0066     }
0067 
0068     /// extract isha, vfs and calchan values, as well as filling the histogram objects
0069     Histo histo_temp;
0070     histo_temp.first = *ihis;
0071     histo_temp.first->Sumw2();
0072     histo_temp.second = (*ihis)->GetTitle();
0073     histo_[title.extraInfo()].resize(2);
0074     if (title.channel() % 2 == 0)
0075       histo_[title.extraInfo()][0] = histo_temp;
0076     else
0077       histo_[title.extraInfo()][1] = histo_temp;
0078   }
0079 }
0080 // ----------------------------------------------------------------------------
0081 //
0082 void CalibrationScanAlgorithm::analyse() {
0083   ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2", "Migrad");
0084   ROOT::Math::MinimizerOptions::SetDefaultStrategy(0);
0085 
0086   if (!cal_) {
0087     edm::LogWarning(mlCommissioning_) << "[CalibrationScanAlgorithm::" << __func__ << "]"
0088                                       << " NULL pointer to derived Analysis object!";
0089     return;
0090   }
0091 
0092   ////////
0093   TFitResultPtr fit_result;
0094   TF1* fit_function_turnOn = nullptr;
0095   TF1* fit_function_decay = nullptr;
0096   TF1* fit_function_deco = nullptr;
0097   if (cal_->deconv_) {
0098     fit_function_deco = new TF1("fit_function_deco", fdeconv, 0, 400, 7);
0099     fit_function_deco->SetParameters(4, 25, 25, 50, 250, 25, 0.75);
0100   } else {
0101     fit_function_turnOn = new TF1("fit_function_turnOn", fturnOn, 0, 400, 4);
0102     fit_function_decay = new TF1("fit_function_decay", fdecay, 0, 400, 3);
0103     fit_function_turnOn->SetParameters(50, 50, 40, 20);
0104     fit_function_decay->SetParameters(-150, -0.01, -0.1);
0105   }
0106 
0107   /// loop over histograms for this fiber
0108   for (auto map_element : histo_) {
0109     // add to the analysis result
0110     cal_->addOneCalibrationPoint(map_element.first);
0111 
0112     // stored as integer the scanned isha and vfs values
0113     std::vector<std::string> tokens;
0114     std::string token;
0115     std::istringstream tokenStream(map_element.first);
0116     while (std::getline(tokenStream, token, '_')) {
0117       tokens.push_back(token);
0118     }
0119 
0120     scanned_isha_.push_back(std::stoi(tokens.at(1)));
0121     scanned_vfs_.push_back(std::stoi(tokens.at(3)));
0122 
0123     // loop on APVs
0124     for (size_t iapv = 0; iapv < 2; iapv++) {
0125       if (!map_element.second[iapv].first) {
0126         edm::LogWarning(mlCommissioning_)
0127             << " NULL pointer to histogram for: " << map_element.second[iapv].second << " !";
0128         return;
0129       }
0130 
0131       cal_->amplitude_[map_element.first][iapv] = 0;
0132       cal_->baseline_[map_element.first][iapv] = 0;
0133       cal_->riseTime_[map_element.first][iapv] = 0;
0134       cal_->turnOn_[map_element.first][iapv] = 0;
0135       cal_->peakTime_[map_element.first][iapv] = 0;
0136       cal_->undershoot_[map_element.first][iapv] = 0;
0137       cal_->tail_[map_element.first][iapv] = 0;
0138       cal_->decayTime_[map_element.first][iapv] = 0;
0139       cal_->smearing_[map_element.first][iapv] = 0;
0140       cal_->chi2_[map_element.first][iapv] = 0;
0141       cal_->isvalid_[map_element.first][iapv] = true;
0142 
0143       if (map_element.second[iapv].first->Integral() == 0) {
0144         cal_->isvalid_[map_element.first][iapv] = false;
0145         continue;
0146       }
0147 
0148       // rescale the plot
0149       correctDistribution(map_element.second[iapv].first, false);
0150 
0151       // from NOTE2009_021 : The charge injection provided by the calibration circuit is known with a precision of 5%;
0152       float error = (map_element.second[iapv].first->GetMaximum() * 0.05);
0153       for (int i = 1; i <= map_element.second[iapv].first->GetNbinsX(); ++i)
0154         map_element.second[iapv].first->SetBinError(i, error);
0155 
0156       /////
0157       if (cal_->deconv_) {  // deconvolution mode
0158         fit_function_deco->SetParameters(4, 25, 25, 50, 250, 25, 0.75);
0159         fit_result = map_element.second[iapv].first->Fit(fit_function_deco, "QRS");
0160 
0161         if (not fit_result.Get()) {
0162           cal_->isvalid_[map_element.first][iapv] = false;
0163           continue;
0164         }
0165 
0166         /// make the fit
0167         float maximum_ampl = fit_function_deco->GetMaximum();
0168         float peak_time = fit_function_deco->GetMaximumX();
0169         float baseline = baseLine(fit_function_deco);
0170         float turn_on_time = turnOn(fit_function_deco, baseline);
0171         float rise_time = peak_time - turn_on_time;
0172 
0173         // start filling info
0174         cal_->amplitude_[map_element.first][iapv] = maximum_ampl - baseline;
0175         cal_->baseline_[map_element.first][iapv] = baseline;
0176         cal_->riseTime_[map_element.first][iapv] = rise_time;
0177         cal_->turnOn_[map_element.first][iapv] = turn_on_time;
0178         cal_->peakTime_[map_element.first][iapv] = peak_time;
0179         if (fit_function_deco->GetMinimumX() > rise_time)
0180           cal_->undershoot_[map_element.first][iapv] =
0181               100 * (fit_function_deco->GetMinimum() - baseline) / (maximum_ampl - baseline);
0182         else
0183           cal_->undershoot_[map_element.first][iapv] = 0;
0184 
0185         // Bin related to peak + 125 ns
0186         int lastBin = map_element.second[iapv].first->FindBin(peak_time + 125);
0187         if (lastBin > map_element.second[iapv].first->GetNbinsX() - 4)
0188           lastBin = map_element.second[iapv].first->GetNbinsX() - 4;
0189 
0190         // tail is the amplitude at 5 bx from the maximum
0191         cal_->tail_[map_element.first][iapv] =
0192             100 * (map_element.second[iapv].first->GetBinContent(lastBin) - baseline) / (maximum_ampl - baseline);
0193 
0194         // reaches 1/e of the peak amplitude
0195         cal_->decayTime_[map_element.first][iapv] = decayTime(fit_function_deco) - peak_time;
0196         cal_->smearing_[map_element.first][iapv] = 0;
0197         cal_->chi2_[map_element.first][iapv] =
0198             fit_function_deco->GetChisquare() /
0199             (map_element.second[iapv].first->GetNbinsX() - fit_function_deco->GetNpar());
0200 
0201       } else {
0202         // peak mode
0203         fit_function_turnOn->SetParameters(50, 50, 40, 20);
0204         fit_function_turnOn->SetRange(fit_function_turnOn->GetXmin(),
0205                                       map_element.second[iapv].first->GetBinCenter(
0206                                           map_element.second[iapv].first->GetMaximumBin()));  // up to the maximum
0207         fit_result = map_element.second[iapv].first->Fit(fit_function_turnOn, "QSR");
0208         if (not fit_result.Get()) {
0209           cal_->isvalid_[map_element.first][iapv] = false;
0210           continue;
0211         }
0212 
0213         /// make the fit
0214         float maximum_ampl = fit_function_turnOn->GetMaximum();
0215         float peak_time = fit_function_turnOn->GetMaximumX();
0216         float baseline = baseLine(fit_function_turnOn);
0217         float turn_on_time = turnOn(fit_function_turnOn, baseline);
0218         float rise_time = peak_time - turn_on_time;
0219 
0220         // start filling info
0221         cal_->amplitude_[map_element.first][iapv] = maximum_ampl - baseline;
0222         cal_->baseline_[map_element.first][iapv] = baseline;
0223         cal_->riseTime_[map_element.first][iapv] = rise_time;
0224         cal_->turnOn_[map_element.first][iapv] = turn_on_time;
0225         cal_->peakTime_[map_element.first][iapv] = peak_time;
0226 
0227         fit_function_decay->SetParameters(-150, -0.01, -0.1);
0228         fit_function_decay->SetRange(
0229             map_element.second[iapv].first->GetBinCenter(map_element.second[iapv].first->GetMaximumBin()) + 10.,
0230             fit_function_decay->GetXmax());  // up to the maximum
0231         fit_result = map_element.second[iapv].first->Fit(fit_function_decay, "QSR+");
0232 
0233         if (fit_result.Get() and fit_result->Status() >= 4) {
0234           cal_->isvalid_[map_element.first][iapv] = false;
0235           continue;
0236         }
0237 
0238         cal_->undershoot_[map_element.first][iapv] = 0;
0239 
0240         // Bin related to peak + 125 ns
0241         int lastBin = map_element.second[iapv].first->FindBin(peak_time + 125);
0242         if (lastBin > map_element.second[iapv].first->GetNbinsX() - 4)
0243           lastBin = map_element.second[iapv].first->GetNbinsX() - 4;
0244 
0245         // tail is the amplitude at 5 bx from the maximum
0246         cal_->tail_[map_element.first][iapv] =
0247             100 * (map_element.second[iapv].first->GetBinContent(lastBin) - baseline) / (maximum_ampl - baseline);
0248 
0249         // reaches 1/e of the peak amplitude
0250         cal_->decayTime_[map_element.first][iapv] = decayTime(fit_function_decay) - peak_time;
0251         cal_->smearing_[map_element.first][iapv] = 0;
0252         cal_->chi2_[map_element.first][iapv] =
0253             (fit_function_turnOn->GetChisquare() + fit_function_decay->GetChisquare()) /
0254             (map_element.second[iapv].first->GetNbinsX() - fit_function_turnOn->GetNpar() -
0255              fit_function_decay->GetNpar());
0256 
0257         // apply quality requirements
0258         bool isvalid = true;
0259         if (cal_->amplitude_[map_element.first][iapv] < CalibrationScanAnalysis::minAmplitudeThreshold_)
0260           isvalid = false;
0261         if (cal_->baseline_[map_element.first][iapv] < CalibrationScanAnalysis::minBaselineThreshold_)
0262           isvalid = false;
0263         else if (cal_->baseline_[map_element.first][iapv] > CalibrationScanAnalysis::maxBaselineThreshold_)
0264           isvalid = false;
0265         if (cal_->decayTime_[map_element.first][iapv] < CalibrationScanAnalysis::minDecayTimeThreshold_)
0266           isvalid = false;
0267         else if (cal_->decayTime_[map_element.first][iapv] > CalibrationScanAnalysis::maxDecayTimeThreshold_)
0268           isvalid = false;
0269         if (cal_->peakTime_[map_element.first][iapv] < CalibrationScanAnalysis::minPeakTimeThreshold_)
0270           isvalid = false;
0271         else if (cal_->peakTime_[map_element.first][iapv] > CalibrationScanAnalysis::maxPeakTimeThreshold_)
0272           isvalid = false;
0273         if (cal_->riseTime_[map_element.first][iapv] < CalibrationScanAnalysis::minRiseTimeThreshold_)
0274           isvalid = false;
0275         else if (cal_->riseTime_[map_element.first][iapv] > CalibrationScanAnalysis::maxRiseTimeThreshold_)
0276           isvalid = false;
0277         if (cal_->turnOn_[map_element.first][iapv] < CalibrationScanAnalysis::minTurnOnThreshold_)
0278           isvalid = false;
0279         else if (cal_->turnOn_[map_element.first][iapv] > CalibrationScanAnalysis::maxTurnOnThreshold_)
0280           isvalid = false;
0281         if (cal_->chi2_[map_element.first][iapv] > CalibrationScanAnalysis::maxChi2Threshold_)
0282           isvalid = false;
0283 
0284         if (not isvalid) {
0285           cal_->amplitude_[map_element.first][iapv] = 0;
0286           cal_->baseline_[map_element.first][iapv] = 0;
0287           cal_->riseTime_[map_element.first][iapv] = 0;
0288           cal_->turnOn_[map_element.first][iapv] = 0;
0289           cal_->peakTime_[map_element.first][iapv] = 0;
0290           cal_->undershoot_[map_element.first][iapv] = 0;
0291           cal_->tail_[map_element.first][iapv] = 0;
0292           cal_->decayTime_[map_element.first][iapv] = 0;
0293           cal_->smearing_[map_element.first][iapv] = 0;
0294           cal_->chi2_[map_element.first][iapv] = 0;
0295           cal_->isvalid_[map_element.first][iapv] = false;
0296         }
0297       }
0298     }
0299   }
0300 
0301   if (fit_function_deco)
0302     delete fit_function_deco;
0303   if (fit_function_decay)
0304     delete fit_function_decay;
0305   if (fit_function_turnOn)
0306     delete fit_function_turnOn;
0307 }
0308 
0309 // ------
0310 void CalibrationScanAlgorithm::correctDistribution(TH1* histo, const bool& isShape) const {
0311   // 5 events per point in the TM loop
0312   // total signal is obtained by summing 16 strips of the same calChan
0313   if (not isShape) {
0314     for (int iBin = 0; iBin < histo->GetNbinsX(); iBin++) {
0315       histo->SetBinContent(iBin + 1, -histo->GetBinContent(iBin + 1) / 16.);
0316       histo->SetBinContent(iBin + 1, histo->GetBinContent(iBin + 1) / 5.);
0317     }
0318   } else
0319     histo->Scale(1. / histo->Integral());
0320 }
0321 
0322 // ----------------------------------------------------------------------------
0323 float CalibrationScanAlgorithm::baseLine(TF1* f) {
0324   float x = f->GetXmin();
0325   float xmax = 10;
0326   float baseline = 0;
0327   int npoints = 0;
0328   for (; x < xmax; x += 0.1) {
0329     baseline += f->Eval(x);
0330     npoints++;
0331   }
0332   return baseline / npoints;
0333 }
0334 
0335 // ----------------------------------------------------------------------------
0336 float CalibrationScanAlgorithm::turnOn(
0337     TF1* f, const float& baseline) {  // should happen within 100 ns in both deco and peak modes
0338   float max_amplitude = f->GetMaximum();
0339   float time = 10.;
0340   for (; time < 100 && (f->Eval(time) - baseline) < 0.05 * (max_amplitude - baseline); time += 0.1) {
0341   }  // flucutation higher than 5% of the pulse height
0342   return time;
0343 }
0344 
0345 // ----------------------------------------------------------------------------
0346 float CalibrationScanAlgorithm::decayTime(
0347     TF1* f) {  // if we approximate the decay to an exp(-t/tau), in one constant unit, the amplited is reduced by e^{-1}
0348   float xval = std::max(f->GetXmin(), f->GetMaximumX());
0349   float max_amplitude = f->GetMaximum();
0350   float x = xval;
0351   for (; x < 1000;
0352        x = x +
0353            0.1) {  // 1000 is a reasoable large bound to compute the decay time .. in case the function is bad it is useful to break the loop
0354     if (f->Eval(x) < max_amplitude * exp(-1))
0355       break;
0356   }
0357   return x;
0358 }
0359 
0360 // --- function to extract the VFS value corresponding to decay time of 125ns, then ISHA close to 50 ns
0361 void CalibrationScanAlgorithm::tuneIndependently(const int& iapv,
0362                                                  const float& targetRiseTime,
0363                                                  const float& targetDecayTime) {
0364   std::map<int, std::vector<float> > decayTime_vs_vfs;
0365   TString name;
0366   int imap = 0;
0367 
0368   for (auto map_element : histo_) {
0369     // only consider isha values in the middle of the scanned range
0370     if (scanned_isha_.at(imap) <= CalibrationScanAnalysis::minISHAforVFSTune_ or
0371         scanned_isha_.at(imap) >= CalibrationScanAnalysis::maxISHAforVFSTune_) {
0372       imap++;
0373       continue;
0374     }
0375 
0376     if (cal_->isValid(map_element.first)[iapv])
0377       decayTime_vs_vfs[scanned_vfs_.at(imap)].push_back(cal_->decayTime(map_element.first)[iapv]);
0378 
0379     if (name == "") {  // store the base name
0380       name = Form("%s", map_element.second[iapv].first->GetName());
0381       name.ReplaceAll("_" + map_element.first, "");
0382     }
0383     imap++;
0384   }
0385 
0386   // sort before taking the median
0387   for (auto iter : decayTime_vs_vfs)
0388     sort(iter.second.begin(), iter.second.end());
0389 
0390   name.ReplaceAll("ExpertHisto_", "");
0391 
0392   // transform the dependance vs vfs in graph
0393   cal_->decayTime_vs_vfs_.push_back(new TGraph());
0394   cal_->decayTime_vs_vfs_.back()->SetName(Form("decayTime_%s", name.Data()));
0395 
0396   // transform the dependance vs isha in graph
0397   cal_->riseTime_vs_isha_.push_back(new TGraph());
0398   cal_->riseTime_vs_isha_.back()->SetName(Form("riseTime_%s", name.Data()));
0399 
0400   if (!decayTime_vs_vfs.empty()) {
0401     int ipoint = 0;
0402     for (auto map_element : decayTime_vs_vfs) {
0403       if (!map_element.second.empty()) {
0404         cal_->decayTime_vs_vfs_.at(iapv)->SetPoint(
0405             ipoint, map_element.second.at(round(map_element.second.size() / 2)), map_element.first);
0406         ipoint++;
0407       }
0408     }
0409 
0410     double max_apv =
0411         TMath::MaxElement(cal_->decayTime_vs_vfs_.at(iapv)->GetN(), cal_->decayTime_vs_vfs_.at(iapv)->GetY());
0412     double min_apv =
0413         TMath::MinElement(cal_->decayTime_vs_vfs_.at(iapv)->GetN(), cal_->decayTime_vs_vfs_.at(iapv)->GetY());
0414 
0415     cal_->vfs_[iapv] = cal_->decayTime_vs_vfs_.at(iapv)->Eval(targetDecayTime);
0416 
0417     // avoid extrapolations
0418     if (cal_->vfs_[iapv] < min_apv)
0419       cal_->vfs_[iapv] = min_apv;
0420     else if (cal_->vfs_[iapv] > max_apv)
0421       cal_->vfs_[iapv] = max_apv;
0422 
0423     // value for each isha but different ISHA
0424     std::map<int, std::vector<float> > riseTime_vs_isha;
0425     imap = 0;
0426     // store for each isha value all rise time (changing isha)
0427     for (auto map_element : histo_) {
0428       if (fabs(scanned_vfs_.at(imap) - cal_->vfs_[iapv]) < CalibrationScanAnalysis::VFSrange_ and
0429           cal_->isValid(map_element.first)[iapv])  //around chosen VFS by \pm 20
0430         riseTime_vs_isha[scanned_isha_.at(imap)].push_back(cal_->riseTime(map_element.first)[iapv]);
0431       if (name == "") {
0432         name = Form("%s", map_element.second[iapv].first->GetName());
0433         name.ReplaceAll("_" + map_element.first, "");
0434       }
0435       imap++;
0436     }
0437 
0438     // sort before taking the median
0439     for (auto iter : riseTime_vs_isha)
0440       sort(iter.second.begin(), iter.second.end());
0441     name.ReplaceAll("ExpertHisto_", "");
0442 
0443     ////
0444     if (!riseTime_vs_isha.empty()) {
0445       int ipoint = 0;
0446       for (auto map_element : riseTime_vs_isha) {
0447         if (!map_element.second.empty()) {
0448           cal_->riseTime_vs_isha_.at(iapv)->SetPoint(
0449               ipoint, map_element.second.at(round(map_element.second.size() / 2)), map_element.first);
0450           ipoint++;
0451         }
0452       }
0453 
0454       double max_apv =
0455           TMath::MaxElement(cal_->riseTime_vs_isha_.at(iapv)->GetN(), cal_->riseTime_vs_isha_.at(iapv)->GetY());
0456       double min_apv =
0457           TMath::MinElement(cal_->riseTime_vs_isha_.at(iapv)->GetN(), cal_->riseTime_vs_isha_.at(iapv)->GetY());
0458 
0459       cal_->isha_[iapv] = cal_->riseTime_vs_isha_.at(iapv)->Eval(targetRiseTime);
0460 
0461       if (cal_->isha_[iapv] < min_apv)
0462         cal_->isha_[iapv] = min_apv;
0463       else if (cal_->isha_[iapv] > max_apv)
0464         cal_->isha_[iapv] = max_apv;
0465     } else
0466       cal_->isha_[iapv] = -1;
0467   }
0468 }
0469 
0470 ////////////// Simultaneously tune isha and vfs
0471 void CalibrationScanAlgorithm::tuneSimultaneously(const int& iapv,
0472                                                   const float& targetRiseTime,
0473                                                   const float& targetDecayTime) {
0474   // Build 2D graph for each APV with rise and decay time trend vs ISHA and VFS
0475   cal_->decayTime_vs_isha_vfs_.push_back(new TGraph2D());
0476   cal_->riseTime_vs_isha_vfs_.push_back(new TGraph2D());
0477 
0478   // store for each vfs value all decay time (changing vfs)
0479   TString name_apv;
0480   int ipoint_apv = 0;
0481   int imap = 0;
0482 
0483   for (auto map_element : histo_) {
0484     if (cal_->isValid(map_element.first)[iapv]) {
0485       cal_->decayTime_vs_isha_vfs_.at(iapv)->SetPoint(
0486           ipoint_apv, scanned_isha_.at(imap), scanned_vfs_.at(imap), cal_->decayTime(map_element.first)[iapv]);
0487       cal_->riseTime_vs_isha_vfs_.at(iapv)->SetPoint(
0488           ipoint_apv, scanned_isha_.at(imap), scanned_vfs_.at(imap), cal_->riseTime(map_element.first)[iapv]);
0489       ipoint_apv++;
0490     }
0491     if (name_apv == "") {  // store the base name
0492       name_apv = Form("%s", map_element.second[iapv].first->GetName());
0493       name_apv.ReplaceAll("_" + map_element.first, "");
0494     }
0495     imap++;
0496   }
0497 
0498   name_apv.ReplaceAll("ExpertHisto_", "");
0499 
0500   cal_->decayTime_vs_isha_vfs_.at(iapv)->SetName(Form("decayTime_%s", name_apv.Data()));
0501   cal_->riseTime_vs_isha_vfs_.at(iapv)->SetName(Form("riseTime_%s", name_apv.Data()));
0502 
0503   // Define 2D histogram for the distance between values and target
0504   TH2F* hist_decay_apv = new TH2F("hist_decay_apv",
0505                                   "hist_decay_apv",
0506                                   500,
0507                                   *min_element(scanned_isha_.begin(), scanned_isha_.end()),
0508                                   *max_element(scanned_isha_.begin(), scanned_isha_.end()),
0509                                   500,
0510                                   *min_element(scanned_vfs_.begin(), scanned_vfs_.end()),
0511                                   *max_element(scanned_vfs_.begin(), scanned_vfs_.end()));
0512 
0513   TH2F* hist_rise_apv = (TH2F*)hist_decay_apv->Clone();
0514   hist_rise_apv->SetName("hist_rise_apv");
0515   hist_rise_apv->Reset();
0516 
0517   TH2F* hist_distance = (TH2F*)hist_decay_apv->Clone();
0518   hist_distance->SetName("hist_distance");
0519   hist_distance->Reset();
0520 
0521   for (int iBin = 1; iBin <= hist_decay_apv->GetNbinsX(); iBin++) {
0522     for (int jBin = 1; jBin <= hist_decay_apv->GetNbinsY(); jBin++) {
0523       if (ipoint_apv != 0) {
0524         if (cal_->decayTime_vs_isha_vfs_.at(iapv)->GetN() > 10)  // to make sure the interpolation can work
0525           hist_decay_apv->SetBinContent(
0526               iBin,
0527               jBin,
0528               cal_->decayTime_vs_isha_vfs_.at(iapv)->Interpolate(hist_decay_apv->GetXaxis()->GetBinCenter(iBin),
0529                                                                  hist_decay_apv->GetYaxis()->GetBinCenter(jBin)));
0530         if (cal_->riseTime_vs_isha_vfs_.at(iapv)->GetN() > 10)
0531           hist_rise_apv->SetBinContent(
0532               iBin,
0533               jBin,
0534               cal_->riseTime_vs_isha_vfs_.at(iapv)->Interpolate(hist_rise_apv->GetXaxis()->GetBinCenter(iBin),
0535                                                                 hist_rise_apv->GetYaxis()->GetBinCenter(jBin)));
0536       }
0537     }
0538   }
0539 
0540   // further smoothing --> a smooth behaviour is indeed expected
0541   hist_decay_apv->Smooth();
0542   hist_rise_apv->Smooth();
0543 
0544   for (int iBin = 1; iBin <= hist_decay_apv->GetNbinsX(); iBin++) {
0545     for (int jBin = 1; jBin <= hist_decay_apv->GetNbinsY(); jBin++) {
0546       hist_distance->SetBinContent(
0547           iBin,
0548           jBin,
0549           sqrt(pow((hist_decay_apv->GetBinContent(iBin, jBin) - targetDecayTime) / targetDecayTime, 2) +
0550                pow((hist_rise_apv->GetBinContent(iBin, jBin) - targetRiseTime) / targetRiseTime, 2)));
0551     }
0552   }
0553 
0554   int minx, miny, minz;
0555   hist_distance->GetMinimumBin(minx, miny, minz);
0556 
0557   cal_->isha_[iapv] = round(hist_distance->GetXaxis()->GetBinCenter(minx));
0558   cal_->vfs_[iapv] = round(hist_distance->GetYaxis()->GetBinCenter(miny));
0559 
0560   delete hist_decay_apv;
0561   delete hist_rise_apv;
0562   delete hist_distance;
0563 }
0564 
0565 void CalibrationScanAlgorithm::fillTunedObservables(const int& apvid) {
0566   // find the closest isha and vfs for each APV
0567   int distance_apv = 10000;
0568 
0569   // find close by ISHA
0570   for (size_t i = 0; i < scanned_isha_.size(); i++) {
0571     if (fabs(scanned_isha_.at(i) - cal_->bestISHA().at(apvid)) < distance_apv) {
0572       distance_apv = fabs(scanned_isha_.at(i) - cal_->bestISHA().at(apvid));
0573       cal_->tunedISHA_.at(apvid) = scanned_isha_.at(i);
0574     }
0575   }
0576 
0577   distance_apv = 10000;
0578 
0579   // find close by VFS
0580   for (size_t i = 0; i < scanned_vfs_.size(); i++) {
0581     if (fabs(scanned_vfs_.at(i) - cal_->bestVFS().at(apvid)) < distance_apv) {
0582       distance_apv = fabs(scanned_vfs_.at(i) - cal_->bestVFS().at(apvid));
0583       cal_->tunedVFS_.at(apvid) = scanned_vfs_.at(i);
0584     }
0585   }
0586 
0587   ///
0588   std::string key_apv = std::string(Form("isha_%d_vfs_%d", cal_->tunedISHA().at(apvid), cal_->tunedVFS().at(apvid)));
0589   if (!cal_->amplitude(key_apv).empty()) {
0590     cal_->tunedAmplitude_[apvid] = cal_->amplitude(key_apv)[apvid];
0591     cal_->tunedTail_[apvid] = cal_->tail(key_apv)[apvid];
0592     cal_->tunedRiseTime_[apvid] = cal_->riseTime(key_apv)[apvid];
0593     cal_->tunedDecayTime_[apvid] = cal_->decayTime(key_apv)[apvid];
0594     cal_->tunedTurnOn_[apvid] = cal_->turnOn(key_apv)[apvid];
0595     cal_->tunedPeakTime_[apvid] = cal_->peakTime(key_apv)[apvid];
0596     cal_->tunedUndershoot_[apvid] = cal_->undershoot(key_apv)[apvid];
0597     cal_->tunedBaseline_[apvid] = cal_->baseline(key_apv)[apvid];
0598     cal_->tunedSmearing_[apvid] = cal_->smearing(key_apv)[apvid];
0599     cal_->tunedChi2_[apvid] = cal_->chi2(key_apv)[apvid];
0600   } else {
0601     cal_->tunedAmplitude_[apvid] = 0;
0602     cal_->tunedTail_[apvid] = 0;
0603     cal_->tunedRiseTime_[apvid] = 0;
0604     cal_->tunedDecayTime_[apvid] = 0;
0605     cal_->tunedTurnOn_[apvid] = 0;
0606     cal_->tunedPeakTime_[apvid] = 0;
0607     cal_->tunedUndershoot_[apvid] = 0;
0608     cal_->tunedBaseline_[apvid] = 0;
0609     cal_->tunedSmearing_[apvid] = 0;
0610     cal_->tunedChi2_[apvid] = 0;
0611   }
0612 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!