File indexing completed on 2024-04-06 11:56:49
0001 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResidualsBfieldAngleFitter.h"
0002
0003 static TMinuit *MuonResidualsBfieldAngleFitter_TMinuit;
0004
0005 void MuonResidualsBfieldAngleFitter::inform(TMinuit *tMinuit) { MuonResidualsBfieldAngleFitter_TMinuit = tMinuit; }
0006
0007 void MuonResidualsBfieldAngleFitter_FCN(int &npar, double *gin, double &fval, double *par, int iflag) {
0008 MuonResidualsFitterFitInfo *fitinfo =
0009 (MuonResidualsFitterFitInfo *)(MuonResidualsBfieldAngleFitter_TMinuit->GetObjectFit());
0010 MuonResidualsFitter *fitter = fitinfo->fitter();
0011
0012 fval = 0.;
0013 for (std::vector<double *>::const_iterator resiter = fitter->residuals_begin(); resiter != fitter->residuals_end();
0014 ++resiter) {
0015 const double residual = (*resiter)[MuonResidualsBfieldAngleFitter::kResidual];
0016 const double qoverpt = (*resiter)[MuonResidualsBfieldAngleFitter::kQoverPt];
0017 const double qoverpz = (*resiter)[MuonResidualsBfieldAngleFitter::kQoverPz];
0018
0019 double center = 0.;
0020 center += par[MuonResidualsBfieldAngleFitter::kAngle];
0021 center += par[MuonResidualsBfieldAngleFitter::kBfrompt] * qoverpt;
0022 center += par[MuonResidualsBfieldAngleFitter::kBfrompz] * qoverpz;
0023 center += par[MuonResidualsBfieldAngleFitter::kdEdx] * (1. / qoverpt / qoverpt + 1. / qoverpz / qoverpz) *
0024 (qoverpt > 0. ? 1. : -1.);
0025
0026 if (fitter->residualsModel() == MuonResidualsFitter::kPureGaussian) {
0027 fval += -MuonResidualsFitter_logPureGaussian(residual, center, par[MuonResidualsBfieldAngleFitter::kSigma]);
0028 } else if (fitter->residualsModel() == MuonResidualsFitter::kPowerLawTails) {
0029 fval += -MuonResidualsFitter_logPowerLawTails(
0030 residual, center, par[MuonResidualsBfieldAngleFitter::kSigma], par[MuonResidualsBfieldAngleFitter::kGamma]);
0031 } else if (fitter->residualsModel() == MuonResidualsFitter::kROOTVoigt) {
0032 fval += -MuonResidualsFitter_logROOTVoigt(
0033 residual, center, par[MuonResidualsBfieldAngleFitter::kSigma], par[MuonResidualsBfieldAngleFitter::kGamma]);
0034 } else if (fitter->residualsModel() == MuonResidualsFitter::kGaussPowerTails) {
0035 fval += -MuonResidualsFitter_logGaussPowerTails(residual, center, par[MuonResidualsBfieldAngleFitter::kSigma]);
0036 } else {
0037 assert(false);
0038 }
0039 }
0040 }
0041
0042 bool MuonResidualsBfieldAngleFitter::fit(Alignable *ali) {
0043 initialize_table();
0044
0045 double sum_x = 0.;
0046 double sum_xx = 0.;
0047 int N = 0;
0048
0049 for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
0050 const double residual = (*resiter)[kResidual];
0051
0052
0053 if (fabs(residual) < 0.1) {
0054 sum_x += residual;
0055 sum_xx += residual * residual;
0056 N++;
0057 }
0058 }
0059
0060 if (N < m_minHits)
0061 return false;
0062
0063
0064 double mean = sum_x / double(N);
0065 double stdev = sqrt(sum_xx / double(N) - pow(sum_x / double(N), 2));
0066
0067
0068 sum_x = 0.;
0069 sum_xx = 0.;
0070 N = 0;
0071 for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
0072 const double residual = (*resiter)[kResidual];
0073 if (mean - 1.5 * stdev < residual && residual < mean + 1.5 * stdev) {
0074 sum_x += residual;
0075 sum_xx += residual * residual;
0076 N++;
0077 }
0078 }
0079 mean = sum_x / double(N);
0080 stdev = sqrt(sum_xx / double(N) - pow(sum_x / double(N), 2));
0081
0082 sum_x = 0.;
0083 sum_xx = 0.;
0084 N = 0;
0085 for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
0086 const double residual = (*resiter)[kResidual];
0087 if (mean - 1.5 * stdev < residual && residual < mean + 1.5 * stdev) {
0088 sum_x += residual;
0089 sum_xx += residual * residual;
0090 N++;
0091 }
0092 }
0093 mean = sum_x / double(N);
0094 stdev = sqrt(sum_xx / double(N) - pow(sum_x / double(N), 2));
0095
0096 std::vector<int> parNum;
0097 std::vector<std::string> parName;
0098 std::vector<double> start;
0099 std::vector<double> step;
0100 std::vector<double> low;
0101 std::vector<double> high;
0102
0103 parNum.push_back(kAngle);
0104 parName.push_back(std::string("angle"));
0105 start.push_back(mean);
0106 step.push_back(0.1);
0107 low.push_back(0.);
0108 high.push_back(0.);
0109 parNum.push_back(kBfrompt);
0110 parName.push_back(std::string("bfrompt"));
0111 start.push_back(0.);
0112 step.push_back(0.1 * stdev / 0.05);
0113 low.push_back(0.);
0114 high.push_back(0.);
0115 parNum.push_back(kBfrompz);
0116 parName.push_back(std::string("bfrompz"));
0117 start.push_back(0.);
0118 step.push_back(0.1 * stdev / 0.05);
0119 low.push_back(0.);
0120 high.push_back(0.);
0121 parNum.push_back(kdEdx);
0122 parName.push_back(std::string("dEdx"));
0123 start.push_back(0.);
0124 step.push_back(0.1 * stdev / 0.05);
0125 low.push_back(0.);
0126 high.push_back(0.);
0127 parNum.push_back(kSigma);
0128 parName.push_back(std::string("sigma"));
0129 start.push_back(stdev);
0130 step.push_back(0.1 * stdev);
0131 low.push_back(0.);
0132 high.push_back(0.);
0133 if (residualsModel() != kPureGaussian && residualsModel() != kGaussPowerTails) {
0134 parNum.push_back(kGamma);
0135 parName.push_back(std::string("gamma"));
0136 start.push_back(stdev);
0137 step.push_back(0.1 * stdev);
0138 low.push_back(0.);
0139 high.push_back(0.);
0140 }
0141
0142 return dofit(&MuonResidualsBfieldAngleFitter_FCN, parNum, parName, start, step, low, high);
0143 }
0144
0145 double MuonResidualsBfieldAngleFitter::plot(std::string name, TFileDirectory *dir, Alignable *ali) {
0146 std::stringstream raw_name, narrowed_name, qoverpt_name, qoverpz_name, psquared_name;
0147 raw_name << name << "_raw";
0148 narrowed_name << name << "_narrowed";
0149 qoverpt_name << name << "_qoverpt";
0150 qoverpz_name << name << "_qoverpz";
0151 psquared_name << name << "_psquared";
0152
0153 TH1F *raw_hist =
0154 dir->make<TH1F>(raw_name.str().c_str(), (raw_name.str() + std::string(" (mrad)")).c_str(), 100, -100., 100.);
0155 TH1F *narrowed_hist = dir->make<TH1F>(
0156 narrowed_name.str().c_str(), (narrowed_name.str() + std::string(" (mrad)")).c_str(), 100, -100., 100.);
0157 TProfile *qoverpt_hist = dir->make<TProfile>(
0158 qoverpt_name.str().c_str(), (qoverpt_name.str() + std::string(" (mrad)")).c_str(), 100, -0.05, 0.05);
0159 TProfile *qoverpz_hist = dir->make<TProfile>(
0160 qoverpz_name.str().c_str(), (qoverpz_name.str() + std::string(" (mrad)")).c_str(), 100, -0.05, 0.05);
0161 TProfile *psquared_hist = dir->make<TProfile>(
0162 psquared_name.str().c_str(), (psquared_name.str() + std::string(" (mrad)")).c_str(), 100, -0.05, 0.05);
0163
0164 narrowed_name << "fit";
0165 qoverpt_name << "fit";
0166 qoverpz_name << "fit";
0167 psquared_name << "fit";
0168
0169 double scale_factor = double(numResiduals()) * (100. - -100.) / 100;
0170
0171 TF1 *narrowed_fit = nullptr;
0172 if (residualsModel() == kPureGaussian) {
0173 narrowed_fit = new TF1(narrowed_name.str().c_str(), MuonResidualsFitter_pureGaussian_TF1, -100., 100., 3);
0174 narrowed_fit->SetParameters(scale_factor, value(kAngle) * 1000., value(kSigma) * 1000.);
0175 narrowed_fit->Write();
0176 } else if (residualsModel() == kPowerLawTails) {
0177 narrowed_fit = new TF1(narrowed_name.str().c_str(), MuonResidualsFitter_powerLawTails_TF1, -100., 100., 4);
0178 narrowed_fit->SetParameters(scale_factor, value(kAngle) * 1000., value(kSigma) * 1000., value(kGamma) * 1000.);
0179 narrowed_fit->Write();
0180 } else if (residualsModel() == kROOTVoigt) {
0181 narrowed_fit = new TF1(narrowed_name.str().c_str(), MuonResidualsFitter_ROOTVoigt_TF1, -100., 100., 4);
0182 narrowed_fit->SetParameters(scale_factor, value(kAngle) * 1000., value(kSigma) * 1000., value(kGamma) * 1000.);
0183 narrowed_fit->Write();
0184 } else if (residualsModel() == kGaussPowerTails) {
0185 narrowed_fit = new TF1(narrowed_name.str().c_str(), MuonResidualsFitter_GaussPowerTails_TF1, -100., 100., 3);
0186 narrowed_fit->SetParameters(scale_factor, value(kAngle) * 1000., value(kSigma) * 1000.);
0187 narrowed_fit->Write();
0188 }
0189
0190 TF1 *qoverpt_fit = new TF1(qoverpt_name.str().c_str(), "[0]+x*[1]", -0.05, 0.05);
0191 qoverpt_fit->SetParameters(value(kAngle) * 1000., value(kBfrompt) * 1000.);
0192 qoverpt_fit->Write();
0193
0194 TF1 *qoverpz_fit = new TF1(qoverpz_name.str().c_str(), "[0]+x*[1]", -0.05, 0.05);
0195 qoverpz_fit->SetParameters(value(kAngle) * 1000., value(kBfrompz) * 1000.);
0196 qoverpz_fit->Write();
0197
0198 TF1 *psquared_fit = new TF1(psquared_name.str().c_str(), "[0]+[1]*x**2", -0.05, 0.05);
0199 psquared_fit->SetParameters(value(kAngle) * 1000., value(kdEdx) * 1000.);
0200 psquared_fit->Write();
0201
0202 for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
0203 const double raw_residual = (*resiter)[kResidual];
0204 const double qoverpt = (*resiter)[kQoverPt];
0205 const double qoverpz = (*resiter)[kQoverPz];
0206 const double psquared = (1. / qoverpt / qoverpt + 1. / qoverpz / qoverpz) * (qoverpt > 0. ? 1. : -1.);
0207
0208 double qoverpt_correction = value(kBfrompt) * qoverpt;
0209 double qoverpz_correction = value(kBfrompz) * qoverpz;
0210 double dEdx_correction = value(kdEdx) * psquared;
0211 double corrected_residual = raw_residual - qoverpt_correction - qoverpz_correction - dEdx_correction;
0212
0213 raw_hist->Fill(raw_residual * 1000.);
0214 narrowed_hist->Fill(corrected_residual * 1000.);
0215
0216 qoverpt_hist->Fill(qoverpt, (raw_residual - qoverpz_correction - dEdx_correction) * 1000.);
0217 qoverpz_hist->Fill(qoverpz, (raw_residual - qoverpt_correction - dEdx_correction) * 1000.);
0218 psquared_hist->Fill(psquared, (raw_residual - qoverpt_correction - qoverpz_correction) * 1000.);
0219 }
0220
0221 return 0.;
0222 }