00001 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResiduals5DOFFitter.h"
00002 #include "TH2F.h"
00003 #include "TMath.h"
00004
00005 static TMinuit *MuonResiduals5DOFFitter_TMinuit;
00006 static double MuonResiduals5DOFFitter_sum_of_weights;
00007 static double MuonResiduals5DOFFitter_number_of_hits;
00008 static bool MuonResiduals5DOFFitter_weightAlignment;
00009
00010 void MuonResiduals5DOFFitter::inform(TMinuit *tMinuit) {
00011 MuonResiduals5DOFFitter_TMinuit = tMinuit;
00012 }
00013
00014 double MuonResiduals5DOFFitter_residual(double delta_x, double delta_z, double delta_phix, double delta_phiy, double delta_phiz, double track_x, double track_y, double track_dxdz, double track_dydz, double alpha, double resslope) {
00015 return delta_x - track_dxdz * delta_z - track_y * track_dxdz * delta_phix + track_x * track_dxdz * delta_phiy - track_y * delta_phiz + resslope * alpha;
00016 }
00017
00018 double MuonResiduals5DOFFitter_resslope(double delta_x, double delta_z, double delta_phix, double delta_phiy, double delta_phiz, double track_x, double track_y, double track_dxdz, double track_dydz) {
00019 return -track_dxdz * track_dydz * delta_phix + (1. + track_dxdz * track_dxdz) * delta_phiy - track_dydz * delta_phiz;
00020 }
00021
00022 Double_t MuonResiduals5DOFFitter_residual_trackx_TF1(Double_t *xvec, Double_t *par) { return 10.*MuonResiduals5DOFFitter_residual(par[0], par[2], par[3], par[4], par[5], xvec[0], par[7], par[8], par[9], par[10], par[11]); }
00023 Double_t MuonResiduals5DOFFitter_residual_tracky_TF1(Double_t *xvec, Double_t *par) { return 10.*MuonResiduals5DOFFitter_residual(par[0], par[2], par[3], par[4], par[5], par[6], xvec[0], par[8], par[9], par[10], par[11]); }
00024 Double_t MuonResiduals5DOFFitter_residual_trackdxdz_TF1(Double_t *xvec, Double_t *par) { return 10.*MuonResiduals5DOFFitter_residual(par[0], par[2], par[3], par[4], par[5], par[6], par[7], xvec[0], par[9], par[10], par[11]); }
00025 Double_t MuonResiduals5DOFFitter_residual_trackdydz_TF1(Double_t *xvec, Double_t *par) { return 10.*MuonResiduals5DOFFitter_residual(par[0], par[2], par[3], par[4], par[5], par[6], par[7], par[8], xvec[0], par[10], par[11]); }
00026
00027 Double_t MuonResiduals5DOFFitter_resslope_trackx_TF1(Double_t *xvec, Double_t *par) { return 1000.*MuonResiduals5DOFFitter_resslope(par[0], par[2], par[3], par[4], par[5], xvec[0], par[7], par[8], par[9]); }
00028 Double_t MuonResiduals5DOFFitter_resslope_tracky_TF1(Double_t *xvec, Double_t *par) { return 1000.*MuonResiduals5DOFFitter_resslope(par[0], par[2], par[3], par[4], par[5], par[6], xvec[0], par[8], par[9]); }
00029 Double_t MuonResiduals5DOFFitter_resslope_trackdxdz_TF1(Double_t *xvec, Double_t *par) { return 1000.*MuonResiduals5DOFFitter_resslope(par[0], par[2], par[3], par[4], par[5], par[6], par[7], xvec[0], par[9]); }
00030 Double_t MuonResiduals5DOFFitter_resslope_trackdydz_TF1(Double_t *xvec, Double_t *par) { return 1000.*MuonResiduals5DOFFitter_resslope(par[0], par[2], par[3], par[4], par[5], par[6], par[7], par[8], xvec[0]); }
00031
00032 void MuonResiduals5DOFFitter_FCN(int &npar, double *gin, double &fval, double *par, int iflag) {
00033 MuonResidualsFitterFitInfo *fitinfo = (MuonResidualsFitterFitInfo*)(MuonResiduals5DOFFitter_TMinuit->GetObjectFit());
00034 MuonResidualsFitter *fitter = fitinfo->fitter();
00035
00036 fval = 0.;
00037 for (std::vector<double*>::const_iterator resiter = fitter->residuals_begin(); resiter != fitter->residuals_end(); ++resiter) {
00038 const double residual = (*resiter)[MuonResiduals5DOFFitter::kResid];
00039 const double resslope = (*resiter)[MuonResiduals5DOFFitter::kResSlope];
00040 const double positionX = (*resiter)[MuonResiduals5DOFFitter::kPositionX];
00041 const double positionY = (*resiter)[MuonResiduals5DOFFitter::kPositionY];
00042 const double angleX = (*resiter)[MuonResiduals5DOFFitter::kAngleX];
00043 const double angleY = (*resiter)[MuonResiduals5DOFFitter::kAngleY];
00044 const double redchi2 = (*resiter)[MuonResiduals5DOFFitter::kRedChi2];
00045
00046 const double alignx = par[MuonResiduals5DOFFitter::kAlignX];
00047 const double alignz = par[MuonResiduals5DOFFitter::kAlignZ];
00048 const double alignphix = par[MuonResiduals5DOFFitter::kAlignPhiX];
00049 const double alignphiy = par[MuonResiduals5DOFFitter::kAlignPhiY];
00050 const double alignphiz = par[MuonResiduals5DOFFitter::kAlignPhiZ];
00051 const double residsigma = par[MuonResiduals5DOFFitter::kResidSigma];
00052 const double resslopesigma = par[MuonResiduals5DOFFitter::kResSlopeSigma];
00053 const double alpha = par[MuonResiduals5DOFFitter::kAlpha];
00054 const double residgamma = par[MuonResiduals5DOFFitter::kResidGamma];
00055 const double resslopegamma = par[MuonResiduals5DOFFitter::kResSlopeGamma];
00056
00057 double residpeak = MuonResiduals5DOFFitter_residual(alignx, alignz, alignphix, alignphiy, alignphiz, positionX, positionY, angleX, angleY, alpha, resslope);
00058 double resslopepeak = MuonResiduals5DOFFitter_resslope(alignx, alignz, alignphix, alignphiy, alignphiz, positionX, positionY, angleX, angleY);
00059
00060 double weight = (1./redchi2) * MuonResiduals5DOFFitter_number_of_hits / MuonResiduals5DOFFitter_sum_of_weights;
00061 if (!MuonResiduals5DOFFitter_weightAlignment) weight = 1.;
00062
00063 if (!MuonResiduals5DOFFitter_weightAlignment || TMath::Prob(redchi2*8, 8) < 0.99) {
00064
00065 if (fitter->residualsModel() == MuonResidualsFitter::kPureGaussian) {
00066 fval += -weight * MuonResidualsFitter_logPureGaussian(residual, residpeak, residsigma);
00067 fval += -weight * MuonResidualsFitter_logPureGaussian(resslope, resslopepeak, resslopesigma);
00068 }
00069 else if (fitter->residualsModel() == MuonResidualsFitter::kPowerLawTails) {
00070 fval += -weight * MuonResidualsFitter_logPowerLawTails(residual, residpeak, residsigma, residgamma);
00071 fval += -weight * MuonResidualsFitter_logPowerLawTails(resslope, resslopepeak, resslopesigma, resslopegamma);
00072 }
00073 else if (fitter->residualsModel() == MuonResidualsFitter::kROOTVoigt) {
00074 fval += -weight * MuonResidualsFitter_logROOTVoigt(residual, residpeak, residsigma, residgamma);
00075 fval += -weight * MuonResidualsFitter_logROOTVoigt(resslope, resslopepeak, resslopesigma, resslopegamma);
00076 }
00077 else if (fitter->residualsModel() == MuonResidualsFitter::kGaussPowerTails) {
00078 fval += -weight * MuonResidualsFitter_logGaussPowerTails(residual, residpeak, residsigma);
00079 fval += -weight * MuonResidualsFitter_logGaussPowerTails(resslope, resslopepeak, resslopesigma);
00080 }
00081 else { assert(false); }
00082
00083 }
00084 }
00085 }
00086
00087 double MuonResiduals5DOFFitter::sumofweights() {
00088 MuonResiduals5DOFFitter_sum_of_weights = 0.;
00089 MuonResiduals5DOFFitter_number_of_hits = 0.;
00090 MuonResiduals5DOFFitter_weightAlignment = m_weightAlignment;
00091 for (std::vector<double*>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
00092 if (m_weightAlignment) {
00093 double redchi2 = (*resiter)[MuonResiduals5DOFFitter::kRedChi2];
00094 if (TMath::Prob(redchi2*8, 8) < 0.99) {
00095 MuonResiduals5DOFFitter_sum_of_weights += 1./redchi2;
00096 MuonResiduals5DOFFitter_number_of_hits += 1.;
00097 }
00098 }
00099 else {
00100 MuonResiduals5DOFFitter_sum_of_weights += 1.;
00101 MuonResiduals5DOFFitter_number_of_hits += 1.;
00102 }
00103 }
00104 return MuonResiduals5DOFFitter_sum_of_weights;
00105 }
00106
00107 bool MuonResiduals5DOFFitter::fit(Alignable *ali) {
00108 initialize_table();
00109 sumofweights();
00110
00111 double resid_mean = 0;
00112 double resslope_mean = 0;
00113 double resid_stdev = 0.5;
00114 double resslope_stdev = 0.005;
00115 double alpha_estimate = 0;
00116
00117 std::vector<int> num;
00118 std::vector<std::string> name;
00119 std::vector<double> start;
00120 std::vector<double> step;
00121 std::vector<double> low;
00122 std::vector<double> high;
00123
00124 if (fixed(kAlignX)) {
00125 num.push_back(kAlignX); name.push_back(std::string("AlignX")); start.push_back(0.); step.push_back(0.1); low.push_back(0.); high.push_back(0.);
00126 } else {
00127 num.push_back(kAlignX); name.push_back(std::string("AlignX")); start.push_back(resid_mean); step.push_back(0.1); low.push_back(0.); high.push_back(0.);
00128 }
00129 num.push_back(kAlignZ); name.push_back(std::string("AlignZ")); start.push_back(0.); step.push_back(0.1); low.push_back(0.); high.push_back(0.);
00130 num.push_back(kAlignPhiX); name.push_back(std::string("AlignPhiX")); start.push_back(0.); step.push_back(0.001); low.push_back(0.); high.push_back(0.);
00131 if (fixed(kAlignPhiY)) {
00132 num.push_back(kAlignPhiY); name.push_back(std::string("AlignPhiY")); start.push_back(0.); step.push_back(0.001); low.push_back(0.); high.push_back(0.);
00133 } else {
00134 num.push_back(kAlignPhiY); name.push_back(std::string("AlignPhiY")); start.push_back(resslope_mean); step.push_back(0.001); low.push_back(0.); high.push_back(0.);
00135 }
00136 num.push_back(kAlignPhiZ); name.push_back(std::string("AlignPhiZ")); start.push_back(0.); step.push_back(0.001); low.push_back(0.); high.push_back(0.);
00137 num.push_back(kResidSigma); name.push_back(std::string("ResidSigma")); start.push_back(resid_stdev); step.push_back(0.01*resid_stdev); low.push_back(0.); high.push_back(0.);
00138 num.push_back(kResSlopeSigma); name.push_back(std::string("ResSlopeSigma")); start.push_back(resslope_stdev); step.push_back(0.01*resslope_stdev); low.push_back(0.); high.push_back(0.);
00139 num.push_back(kAlpha); name.push_back(std::string("Alpha")); start.push_back(alpha_estimate); step.push_back(0.01*resslope_stdev); low.push_back(0.); high.push_back(0.);
00140 if (residualsModel() != kPureGaussian && residualsModel() != kGaussPowerTails) {
00141 num.push_back(kResidGamma); name.push_back(std::string("ResidGamma")); start.push_back(0.1*resid_stdev); step.push_back(0.01*resid_stdev); low.push_back(0.); high.push_back(0.);
00142 num.push_back(kResSlopeGamma); name.push_back(std::string("ResSlopeGamma")); start.push_back(0.1*resslope_stdev); step.push_back(0.01*resslope_stdev); low.push_back(0.); high.push_back(0.);
00143 }
00144
00145 return dofit(&MuonResiduals5DOFFitter_FCN, num, name, start, step, low, high);
00146 }
00147
00148 double MuonResiduals5DOFFitter::plot(std::string name, TFileDirectory *dir, Alignable *ali) {
00149 sumofweights();
00150
00151 std::stringstream name_residual, name_resslope, name_residual_raw, name_resslope_raw, name_residual_cut, name_alpha;
00152 std::stringstream name_residual_trackx, name_resslope_trackx;
00153 std::stringstream name_residual_tracky, name_resslope_tracky;
00154 std::stringstream name_residual_trackdxdz, name_resslope_trackdxdz;
00155 std::stringstream name_residual_trackdydz, name_resslope_trackdydz;
00156
00157 name_residual << name << "_residual";
00158 name_resslope << name << "_resslope";
00159 name_residual_raw << name << "_residual_raw";
00160 name_resslope_raw << name << "_resslope_raw";
00161 name_residual_cut << name << "_residual_cut";
00162 name_alpha << name << "_alpha";
00163 name_residual_trackx << name << "_residual_trackx";
00164 name_resslope_trackx << name << "_resslope_trackx";
00165 name_residual_tracky << name << "_residual_tracky";
00166 name_resslope_tracky << name << "_resslope_tracky";
00167 name_residual_trackdxdz << name << "_residual_trackdxdz";
00168 name_resslope_trackdxdz << name << "_resslope_trackdxdz";
00169 name_residual_trackdydz << name << "_residual_trackdydz";
00170 name_resslope_trackdydz << name << "_resslope_trackdydz";
00171
00172 double width = ali->surface().width();
00173 double length = ali->surface().length();
00174 double min_residual = -100.; double max_residual = 100.;
00175 double min_resslope = -100.; double max_resslope = 100.;
00176 double min_trackx = -width/2.; double max_trackx = width/2.;
00177 double min_tracky = -length/2.; double max_tracky = length/2.;
00178 double min_trackdxdz = -1.5; double max_trackdxdz = 1.5;
00179 double min_trackdydz = -1.5; double max_trackdydz = 1.5;
00180
00181 TH1F *hist_residual = dir->make<TH1F>(name_residual.str().c_str(), "", 100, min_residual, max_residual);
00182 TH1F *hist_resslope = dir->make<TH1F>(name_resslope.str().c_str(), "", 100, min_resslope, max_resslope);
00183 TH1F *hist_residual_raw = dir->make<TH1F>(name_residual_raw.str().c_str(), "", 100, min_residual, max_residual);
00184 TH1F *hist_resslope_raw = dir->make<TH1F>(name_resslope_raw.str().c_str(), "", 100, min_resslope, max_resslope);
00185 TH1F *hist_residual_cut = dir->make<TH1F>(name_residual_cut.str().c_str(), "", 100, min_residual, max_residual);
00186 TH2F *hist_alpha = dir->make<TH2F>(name_alpha.str().c_str(), "", 40, min_resslope, max_resslope, 40, -20., 20.);
00187 TProfile *hist_residual_trackx = dir->make<TProfile>(name_residual_trackx.str().c_str(), "", 100, min_trackx, max_trackx, min_residual, max_residual);
00188 TProfile *hist_resslope_trackx = dir->make<TProfile>(name_resslope_trackx.str().c_str(), "", 100, min_trackx, max_trackx, min_resslope, max_resslope);
00189 TProfile *hist_residual_tracky = dir->make<TProfile>(name_residual_tracky.str().c_str(), "", 100, min_tracky, max_tracky, min_residual, max_residual);
00190 TProfile *hist_resslope_tracky = dir->make<TProfile>(name_resslope_tracky.str().c_str(), "", 100, min_tracky, max_tracky, min_resslope, max_resslope);
00191 TProfile *hist_residual_trackdxdz = dir->make<TProfile>(name_residual_trackdxdz.str().c_str(), "", 500, min_trackdxdz, max_trackdxdz, min_residual, max_residual);
00192 TProfile *hist_resslope_trackdxdz = dir->make<TProfile>(name_resslope_trackdxdz.str().c_str(), "", 500, min_trackdxdz, max_trackdxdz, min_resslope, max_resslope);
00193 TProfile *hist_residual_trackdydz = dir->make<TProfile>(name_residual_trackdydz.str().c_str(), "", 500, min_trackdydz, max_trackdydz, min_residual, max_residual);
00194 TProfile *hist_resslope_trackdydz = dir->make<TProfile>(name_resslope_trackdydz.str().c_str(), "", 500, min_trackdydz, max_trackdydz, min_resslope, max_resslope);
00195
00196 hist_residual_trackx->SetAxisRange(-10., 10., "Y");
00197 hist_resslope_trackx->SetAxisRange(-10., 10., "Y");
00198 hist_residual_tracky->SetAxisRange(-10., 10., "Y");
00199 hist_resslope_tracky->SetAxisRange(-10., 10., "Y");
00200 hist_residual_trackdxdz->SetAxisRange(-10., 10., "Y");
00201 hist_resslope_trackdxdz->SetAxisRange(-10., 10., "Y");
00202 hist_residual_trackdydz->SetAxisRange(-10., 10., "Y");
00203 hist_resslope_trackdydz->SetAxisRange(-10., 10., "Y");
00204
00205 name_residual << "_fit";
00206 name_resslope << "_fit";
00207 name_alpha << "_fit";
00208 name_residual_trackx << "_fit";
00209 name_resslope_trackx << "_fit";
00210 name_residual_tracky << "_fit";
00211 name_resslope_tracky << "_fit";
00212 name_residual_trackdxdz << "_fit";
00213 name_resslope_trackdxdz << "_fit";
00214 name_residual_trackdydz << "_fit";
00215 name_resslope_trackdydz << "_fit";
00216
00217 TF1 *fit_residual = NULL;
00218 TF1 *fit_resslope = NULL;
00219 if (residualsModel() == kPureGaussian) {
00220 fit_residual = new TF1(name_residual.str().c_str(), MuonResidualsFitter_pureGaussian_TF1, min_residual, max_residual, 3);
00221 fit_residual->SetParameters(MuonResiduals5DOFFitter_sum_of_weights * (max_residual - min_residual)/100., 10.*value(kAlignX), 10.*value(kResidSigma));
00222 fit_resslope = new TF1(name_resslope.str().c_str(), MuonResidualsFitter_pureGaussian_TF1, min_resslope, max_resslope, 3);
00223 fit_resslope->SetParameters(MuonResiduals5DOFFitter_sum_of_weights * (max_resslope - min_resslope)/100., 1000.*value(kAlignPhiY), 1000.*value(kResSlopeSigma));
00224 }
00225 else if (residualsModel() == kPowerLawTails) {
00226 fit_residual = new TF1(name_residual.str().c_str(), MuonResidualsFitter_powerLawTails_TF1, min_residual, max_residual, 4);
00227 fit_residual->SetParameters(MuonResiduals5DOFFitter_sum_of_weights * (max_residual - min_residual)/100., 10.*value(kAlignX), 10.*value(kResidSigma), 10.*value(kResidGamma));
00228 fit_resslope = new TF1(name_resslope.str().c_str(), MuonResidualsFitter_powerLawTails_TF1, min_resslope, max_resslope, 4);
00229 fit_resslope->SetParameters(MuonResiduals5DOFFitter_sum_of_weights * (max_resslope - min_resslope)/100., 1000.*value(kAlignPhiY), 1000.*value(kResSlopeSigma), 1000.*value(kResSlopeGamma));
00230 }
00231 else if (residualsModel() == kROOTVoigt) {
00232 fit_residual = new TF1(name_residual.str().c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_residual, max_residual, 4);
00233 fit_residual->SetParameters(MuonResiduals5DOFFitter_sum_of_weights * (max_residual - min_residual)/100., 10.*value(kAlignX), 10.*value(kResidSigma), 10.*value(kResidGamma));
00234 fit_resslope = new TF1(name_resslope.str().c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_resslope, max_resslope, 4);
00235 fit_resslope->SetParameters(MuonResiduals5DOFFitter_sum_of_weights * (max_resslope - min_resslope)/100., 1000.*value(kAlignPhiY), 1000.*value(kResSlopeSigma), 1000.*value(kResSlopeGamma));
00236 }
00237 else if (residualsModel() == kGaussPowerTails) {
00238 fit_residual = new TF1(name_residual.str().c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_residual, max_residual, 3);
00239 fit_residual->SetParameters(MuonResiduals5DOFFitter_sum_of_weights * (max_residual - min_residual)/100., 10.*value(kAlignX), 10.*value(kResidSigma));
00240 fit_resslope = new TF1(name_resslope.str().c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_resslope, max_resslope, 3);
00241 fit_resslope->SetParameters(MuonResiduals5DOFFitter_sum_of_weights * (max_resslope - min_resslope)/100., 1000.*value(kAlignPhiY), 1000.*value(kResSlopeSigma));
00242 }
00243 else { assert(false); }
00244
00245 fit_residual->SetLineColor(2); fit_residual->SetLineWidth(2);
00246 fit_resslope->SetLineColor(2); fit_resslope->SetLineWidth(2);
00247 fit_residual->Write();
00248 fit_resslope->Write();
00249
00250 TF1 *fit_alpha = new TF1(name_alpha.str().c_str(), "[0] + x*[1]", min_resslope, max_resslope);
00251 fit_alpha->SetParameters(10.*value(kAlignX), 10.*value(kAlpha)/1000.);
00252 fit_alpha->SetLineColor(2); fit_alpha->SetLineWidth(2);
00253 fit_alpha->Write();
00254
00255 TProfile *fit_residual_trackx = dir->make<TProfile>(name_residual_trackx.str().c_str(), "", 100, min_trackx, max_trackx);
00256 TProfile *fit_resslope_trackx = dir->make<TProfile>(name_resslope_trackx.str().c_str(), "", 100, min_trackx, max_trackx);
00257 TProfile *fit_residual_tracky = dir->make<TProfile>(name_residual_tracky.str().c_str(), "", 100, min_tracky, max_tracky);
00258 TProfile *fit_resslope_tracky = dir->make<TProfile>(name_resslope_tracky.str().c_str(), "", 100, min_tracky, max_tracky);
00259 TProfile *fit_residual_trackdxdz = dir->make<TProfile>(name_residual_trackdxdz.str().c_str(), "", 500, min_trackdxdz, max_trackdxdz);
00260 TProfile *fit_resslope_trackdxdz = dir->make<TProfile>(name_resslope_trackdxdz.str().c_str(), "", 500, min_trackdxdz, max_trackdxdz);
00261 TProfile *fit_residual_trackdydz = dir->make<TProfile>(name_residual_trackdydz.str().c_str(), "", 500, min_trackdydz, max_trackdydz);
00262 TProfile *fit_resslope_trackdydz = dir->make<TProfile>(name_resslope_trackdydz.str().c_str(), "", 500, min_trackdydz, max_trackdydz);
00263
00264 fit_residual_trackx->SetLineColor(2); fit_residual_trackx->SetLineWidth(2);
00265 fit_resslope_trackx->SetLineColor(2); fit_resslope_trackx->SetLineWidth(2);
00266 fit_residual_tracky->SetLineColor(2); fit_residual_tracky->SetLineWidth(2);
00267 fit_resslope_tracky->SetLineColor(2); fit_resslope_tracky->SetLineWidth(2);
00268 fit_residual_trackdxdz->SetLineColor(2); fit_residual_trackdxdz->SetLineWidth(2);
00269 fit_resslope_trackdxdz->SetLineColor(2); fit_resslope_trackdxdz->SetLineWidth(2);
00270 fit_residual_trackdydz->SetLineColor(2); fit_residual_trackdydz->SetLineWidth(2);
00271 fit_resslope_trackdydz->SetLineColor(2); fit_resslope_trackdydz->SetLineWidth(2);
00272
00273 name_residual_trackx << "line";
00274 name_resslope_trackx << "line";
00275 name_residual_tracky << "line";
00276 name_resslope_tracky << "line";
00277 name_residual_trackdxdz << "line";
00278 name_resslope_trackdxdz << "line";
00279 name_residual_trackdydz << "line";
00280 name_resslope_trackdydz << "line";
00281
00282 TF1 *fitline_residual_trackx = new TF1(name_residual_trackx.str().c_str(), MuonResiduals5DOFFitter_residual_trackx_TF1, min_trackx, max_trackx, 12);
00283 TF1 *fitline_resslope_trackx = new TF1(name_resslope_trackx.str().c_str(), MuonResiduals5DOFFitter_resslope_trackx_TF1, min_trackx, max_trackx, 12);
00284 TF1 *fitline_residual_tracky = new TF1(name_residual_tracky.str().c_str(), MuonResiduals5DOFFitter_residual_tracky_TF1, min_tracky, max_tracky, 12);
00285 TF1 *fitline_resslope_tracky = new TF1(name_resslope_tracky.str().c_str(), MuonResiduals5DOFFitter_resslope_tracky_TF1, min_tracky, max_tracky, 12);
00286 TF1 *fitline_residual_trackdxdz = new TF1(name_residual_trackdxdz.str().c_str(), MuonResiduals5DOFFitter_residual_trackdxdz_TF1, min_trackdxdz, max_trackdxdz, 12);
00287 TF1 *fitline_resslope_trackdxdz = new TF1(name_resslope_trackdxdz.str().c_str(), MuonResiduals5DOFFitter_resslope_trackdxdz_TF1, min_trackdxdz, max_trackdxdz, 12);
00288 TF1 *fitline_residual_trackdydz = new TF1(name_residual_trackdydz.str().c_str(), MuonResiduals5DOFFitter_residual_trackdydz_TF1, min_trackdydz, max_trackdydz, 12);
00289 TF1 *fitline_resslope_trackdydz = new TF1(name_resslope_trackdydz.str().c_str(), MuonResiduals5DOFFitter_resslope_trackdydz_TF1, min_trackdydz, max_trackdydz, 12);
00290
00291 double sum_resslope = 0.;
00292 double sum_trackx = 0.;
00293 double sum_tracky = 0.;
00294 double sum_trackdxdz = 0.;
00295 double sum_trackdydz = 0.;
00296 for (std::vector<double*>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
00297 const double resslope = (*resiter)[MuonResiduals5DOFFitter::kResSlope];
00298 const double positionX = (*resiter)[MuonResiduals5DOFFitter::kPositionX];
00299 const double positionY = (*resiter)[MuonResiduals5DOFFitter::kPositionY];
00300 const double angleX = (*resiter)[MuonResiduals5DOFFitter::kAngleX];
00301 const double angleY = (*resiter)[MuonResiduals5DOFFitter::kAngleY];
00302 const double redchi2 = (*resiter)[MuonResiduals5DOFFitter::kRedChi2];
00303 double weight = 1./redchi2;
00304 if (!m_weightAlignment) weight = 1.;
00305
00306 if (!m_weightAlignment || TMath::Prob(redchi2*8, 8) < 0.99) {
00307 sum_resslope += weight * resslope;
00308 sum_trackx += weight * positionX;
00309 sum_tracky += weight * positionY;
00310 sum_trackdxdz += weight * angleX;
00311 sum_trackdydz += weight * angleY;
00312 }
00313 }
00314 double mean_resslope = sum_resslope / MuonResiduals5DOFFitter_sum_of_weights;
00315 double mean_trackx = sum_trackx / MuonResiduals5DOFFitter_sum_of_weights;
00316 double mean_tracky = sum_tracky / MuonResiduals5DOFFitter_sum_of_weights;
00317 double mean_trackdxdz = sum_trackdxdz / MuonResiduals5DOFFitter_sum_of_weights;
00318 double mean_trackdydz = sum_trackdydz / MuonResiduals5DOFFitter_sum_of_weights;
00319
00320 double fitparameters[12];
00321 fitparameters[0] = value(kAlignX);
00322 fitparameters[1] = 0.;
00323 fitparameters[2] = value(kAlignZ);
00324 fitparameters[3] = value(kAlignPhiX);
00325 fitparameters[4] = value(kAlignPhiY);
00326 fitparameters[5] = value(kAlignPhiZ);
00327 fitparameters[6] = mean_trackx;
00328 fitparameters[7] = mean_tracky;
00329 fitparameters[8] = mean_trackdxdz;
00330 fitparameters[9] = mean_trackdydz;
00331 fitparameters[10] = value(kAlpha);
00332 fitparameters[11] = mean_resslope;
00333
00334 fitline_residual_trackx->SetParameters(fitparameters);
00335 fitline_resslope_trackx->SetParameters(fitparameters);
00336 fitline_residual_tracky->SetParameters(fitparameters);
00337 fitline_resslope_tracky->SetParameters(fitparameters);
00338 fitline_residual_trackdxdz->SetParameters(fitparameters);
00339 fitline_resslope_trackdxdz->SetParameters(fitparameters);
00340 fitline_residual_trackdydz->SetParameters(fitparameters);
00341 fitline_resslope_trackdydz->SetParameters(fitparameters);
00342
00343 fitline_residual_trackx->SetLineColor(2); fitline_residual_trackx->SetLineWidth(2);
00344 fitline_resslope_trackx->SetLineColor(2); fitline_resslope_trackx->SetLineWidth(2);
00345 fitline_residual_tracky->SetLineColor(2); fitline_residual_tracky->SetLineWidth(2);
00346 fitline_resslope_tracky->SetLineColor(2); fitline_resslope_tracky->SetLineWidth(2);
00347 fitline_residual_trackdxdz->SetLineColor(2); fitline_residual_trackdxdz->SetLineWidth(2);
00348 fitline_resslope_trackdxdz->SetLineColor(2); fitline_resslope_trackdxdz->SetLineWidth(2);
00349 fitline_residual_trackdydz->SetLineColor(2); fitline_residual_trackdydz->SetLineWidth(2);
00350 fitline_resslope_trackdydz->SetLineColor(2); fitline_resslope_trackdydz->SetLineWidth(2);
00351
00352 fitline_residual_trackx->Write();
00353 fitline_resslope_trackx->Write();
00354 fitline_residual_tracky->Write();
00355 fitline_resslope_tracky->Write();
00356 fitline_residual_trackdxdz->Write();
00357 fitline_resslope_trackdxdz->Write();
00358 fitline_residual_trackdydz->Write();
00359 fitline_resslope_trackdydz->Write();
00360
00361 for (std::vector<double*>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
00362 const double resid = (*resiter)[MuonResiduals5DOFFitter::kResid];
00363 const double resslope = (*resiter)[MuonResiduals5DOFFitter::kResSlope];
00364 const double positionX = (*resiter)[MuonResiduals5DOFFitter::kPositionX];
00365 const double positionY = (*resiter)[MuonResiduals5DOFFitter::kPositionY];
00366 const double angleX = (*resiter)[MuonResiduals5DOFFitter::kAngleX];
00367 const double angleY = (*resiter)[MuonResiduals5DOFFitter::kAngleY];
00368 const double redchi2 = (*resiter)[MuonResiduals5DOFFitter::kRedChi2];
00369 double weight = 1./redchi2;
00370 if (!m_weightAlignment) weight = 1.;
00371
00372 if (!m_weightAlignment || TMath::Prob(redchi2*8, 8) < 0.99) {
00373 hist_alpha->Fill(1000.*resslope, 10.*resid);
00374
00375 double geom_resid = MuonResiduals5DOFFitter_residual(value(kAlignX), value(kAlignZ), value(kAlignPhiX), value(kAlignPhiY), value(kAlignPhiZ), positionX, positionY, angleX, angleY, value(kAlpha), resslope);
00376 hist_residual->Fill(10.*(resid - geom_resid + value(kAlignX)), weight);
00377 hist_residual_trackx->Fill(positionX, 10.*resid, weight);
00378 hist_residual_tracky->Fill(positionY, 10.*resid, weight);
00379 hist_residual_trackdxdz->Fill(angleX, 10.*resid, weight);
00380 hist_residual_trackdydz->Fill(angleY, 10.*resid, weight);
00381 fit_residual_trackx->Fill(positionX, 10.*geom_resid, weight);
00382 fit_residual_tracky->Fill(positionY, 10.*geom_resid, weight);
00383 fit_residual_trackdxdz->Fill(angleX, 10.*geom_resid, weight);
00384 fit_residual_trackdydz->Fill(angleY, 10.*geom_resid, weight);
00385
00386 double geom_resslope = MuonResiduals5DOFFitter_resslope(value(kAlignX), value(kAlignZ), value(kAlignPhiX), value(kAlignPhiY), value(kAlignPhiZ), positionX, positionY, angleX, angleY);
00387 hist_resslope->Fill(1000.*(resslope - geom_resslope + value(kAlignPhiY)), weight);
00388 hist_resslope_trackx->Fill(positionX, 1000.*resslope, weight);
00389 hist_resslope_tracky->Fill(positionY, 1000.*resslope, weight);
00390 hist_resslope_trackdxdz->Fill(angleX, 1000.*resslope, weight);
00391 hist_resslope_trackdydz->Fill(angleY, 1000.*resslope, weight);
00392 fit_resslope_trackx->Fill(positionX, 1000.*geom_resslope, weight);
00393 fit_resslope_tracky->Fill(positionY, 1000.*geom_resslope, weight);
00394 fit_resslope_trackdxdz->Fill(angleX, 1000.*geom_resslope, weight);
00395 fit_resslope_trackdydz->Fill(angleY, 1000.*geom_resslope, weight);
00396 }
00397
00398 hist_residual_raw->Fill(10.*resid);
00399 hist_resslope_raw->Fill(1000.*resslope);
00400 if (fabs(resslope) < 0.005) hist_residual_cut->Fill(10.*resid);
00401 }
00402
00403 double chi2 = 0.;
00404 double ndof = 0.;
00405 for (int i = 1; i <= hist_residual->GetNbinsX(); i++) {
00406 double xi = hist_residual->GetBinCenter(i);
00407 double yi = hist_residual->GetBinContent(i);
00408 double yerri = hist_residual->GetBinError(i);
00409 double yth = fit_residual->Eval(xi);
00410 if (yerri > 0.) {
00411 chi2 += pow((yth - yi)/yerri, 2);
00412 ndof += 1.;
00413 }
00414 }
00415 for (int i = 1; i <= hist_resslope->GetNbinsX(); i++) {
00416 double xi = hist_resslope->GetBinCenter(i);
00417 double yi = hist_resslope->GetBinContent(i);
00418 double yerri = hist_resslope->GetBinError(i);
00419 double yth = fit_resslope->Eval(xi);
00420 if (yerri > 0.) {
00421 chi2 += pow((yth - yi)/yerri, 2);
00422 ndof += 1.;
00423 }
00424 }
00425 ndof -= npar();
00426
00427 return (ndof > 0. ? chi2 / ndof : -1.);
00428 }
1.7.3