00001
00002
00003 #ifdef STANDALONE_FITTER
00004 #include "MuonResiduals5DOFFitter.h"
00005 #else
00006 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResiduals5DOFFitter.h"
00007 #endif
00008
00009 #include "TH2F.h"
00010 #include "TMath.h"
00011 #include "TTree.h"
00012 #include "TFile.h"
00013
00014 namespace
00015 {
00016 TMinuit *minuit;
00017
00018 double sum_of_weights;
00019 double number_of_hits;
00020 bool weight_alignment;
00021
00022 double residual_x(double delta_x, double delta_z,
00023 double delta_phix, double delta_phiy, double delta_phiz,
00024 double track_x, double track_y, double track_dxdz, double track_dydz,
00025 double alpha, double resslope)
00026 {
00027 return delta_x
00028 - track_dxdz * delta_z
00029 - track_y * track_dxdz * delta_phix
00030 + track_x * track_dxdz * delta_phiy
00031 - track_y * delta_phiz
00032 + resslope * alpha;
00033 }
00034
00035 double residual_dxdz(double delta_x, double delta_z,
00036 double delta_phix, double delta_phiy, double delta_phiz,
00037 double track_x, double track_y, double track_dxdz, double track_dydz)
00038 {
00039 return -track_dxdz * track_dydz * delta_phix
00040 + (1. + track_dxdz * track_dxdz) * delta_phiy
00041 - track_dydz * delta_phiz;
00042 }
00043
00044 Double_t residual_x_trackx_TF1(Double_t *xx, Double_t *p) { return 10.*residual_x(p[0], p[2], p[3], p[4], p[5], xx[0], p[7], p[8], p[9], p[10], p[11]); }
00045 Double_t residual_x_tracky_TF1(Double_t *xx, Double_t *p) { return 10.*residual_x(p[0], p[2], p[3], p[4], p[5], p[6], xx[0], p[8], p[9], p[10], p[11]); }
00046 Double_t residual_x_trackdxdz_TF1(Double_t *xx, Double_t *p) { return 10.*residual_x(p[0], p[2], p[3], p[4], p[5], p[6], p[7], xx[0], p[9], p[10], p[11]); }
00047 Double_t residual_x_trackdydz_TF1(Double_t *xx, Double_t *p) { return 10.*residual_x(p[0], p[2], p[3], p[4], p[5], p[6], p[7], p[8], xx[0], p[10], p[11]); }
00048
00049 Double_t residual_dxdz_trackx_TF1(Double_t *xx, Double_t *p) { return 1000.*residual_dxdz(p[0], p[2], p[3], p[4], p[5], xx[0], p[7], p[8], p[9]); }
00050 Double_t residual_dxdz_tracky_TF1(Double_t *xx, Double_t *p) { return 1000.*residual_dxdz(p[0], p[2], p[3], p[4], p[5], p[6], xx[0], p[8], p[9]); }
00051 Double_t residual_dxdz_trackdxdz_TF1(Double_t *xx, Double_t *p) { return 1000.*residual_dxdz(p[0], p[2], p[3], p[4], p[5], p[6], p[7], xx[0], p[9]); }
00052 Double_t residual_dxdz_trackdydz_TF1(Double_t *xx, Double_t *p) { return 1000.*residual_dxdz(p[0], p[2], p[3], p[4], p[5], p[6], p[7], p[8], xx[0]); }
00053 }
00054
00055
00056 void MuonResiduals5DOFFitter_FCN(int &npar, double *gin, double &fval, double *par, int iflag)
00057 {
00058 MuonResidualsFitterFitInfo *fitinfo = (MuonResidualsFitterFitInfo*)(minuit->GetObjectFit());
00059 MuonResidualsFitter *fitter = fitinfo->fitter();
00060
00061 fval = 0.;
00062 for (std::vector<double*>::const_iterator resiter = fitter->residuals_begin(); resiter != fitter->residuals_end(); ++resiter)
00063 {
00064 const double residual = (*resiter)[MuonResiduals5DOFFitter::kResid];
00065 const double resslope = (*resiter)[MuonResiduals5DOFFitter::kResSlope];
00066 const double positionX = (*resiter)[MuonResiduals5DOFFitter::kPositionX];
00067 const double positionY = (*resiter)[MuonResiduals5DOFFitter::kPositionY];
00068 const double angleX = (*resiter)[MuonResiduals5DOFFitter::kAngleX];
00069 const double angleY = (*resiter)[MuonResiduals5DOFFitter::kAngleY];
00070 const double redchi2 = (*resiter)[MuonResiduals5DOFFitter::kRedChi2];
00071
00072 const double alignx = par[MuonResiduals5DOFFitter::kAlignX];
00073 const double alignz = par[MuonResiduals5DOFFitter::kAlignZ];
00074 const double alignphix = par[MuonResiduals5DOFFitter::kAlignPhiX];
00075 const double alignphiy = par[MuonResiduals5DOFFitter::kAlignPhiY];
00076 const double alignphiz = par[MuonResiduals5DOFFitter::kAlignPhiZ];
00077 const double residsigma = par[MuonResiduals5DOFFitter::kResidSigma];
00078 const double resslopesigma = par[MuonResiduals5DOFFitter::kResSlopeSigma];
00079 const double alpha = par[MuonResiduals5DOFFitter::kAlpha];
00080 const double residgamma = par[MuonResiduals5DOFFitter::kResidGamma];
00081 const double resslopegamma = par[MuonResiduals5DOFFitter::kResSlopeGamma];
00082
00083 double coeff = alpha;
00084 if (fitter->residualsModel() == MuonResidualsFitter::kPureGaussian ||
00085 fitter->residualsModel() == MuonResidualsFitter::kPureGaussian2D) coeff = 0.;
00086 double residpeak = residual_x(alignx, alignz, alignphix, alignphiy, alignphiz, positionX, positionY, angleX, angleY, coeff, resslope);
00087 double resslopepeak = residual_dxdz(alignx, alignz, alignphix, alignphiy, alignphiz, positionX, positionY, angleX, angleY);
00088
00089 double weight = (1./redchi2) * number_of_hits / sum_of_weights;
00090 if (!weight_alignment) weight = 1.;
00091
00092 if (!weight_alignment || TMath::Prob(redchi2*8, 8) < 0.99)
00093 {
00094 if (fitter->residualsModel() == MuonResidualsFitter::kPureGaussian) {
00095 if (fitter->useRes() == MuonResidualsFitter::k1111 || fitter->useRes() == MuonResidualsFitter::k1110 || fitter->useRes() == MuonResidualsFitter::k1010) {
00096 fval += -weight * MuonResidualsFitter_logPureGaussian(residual, residpeak, residsigma);
00097 fval += -weight * MuonResidualsFitter_logPureGaussian(resslope, resslopepeak, resslopesigma);
00098 }
00099 else if (fitter->useRes() == MuonResidualsFitter::k1100) {
00100 fval += -weight * MuonResidualsFitter_logPureGaussian(residual, residpeak, residsigma);
00101 }
00102 else if (fitter->useRes() == MuonResidualsFitter::k0010) {
00103 fval += -weight * MuonResidualsFitter_logPureGaussian(resslope, resslopepeak, resslopesigma);
00104 }
00105 }
00106 else if (fitter->residualsModel() == MuonResidualsFitter::kPureGaussian2D) {
00107 if (fitter->useRes() == MuonResidualsFitter::k1111 || fitter->useRes() == MuonResidualsFitter::k1110 || fitter->useRes() == MuonResidualsFitter::k1010) {
00108 fval += -weight * MuonResidualsFitter_logPureGaussian2D(residual, resslope, residpeak, resslopepeak, residsigma, resslopesigma, alpha);
00109 }
00110 else if (fitter->useRes() == MuonResidualsFitter::k1100) {
00111 fval += -weight * MuonResidualsFitter_logPureGaussian(residual, residpeak, residsigma);
00112 }
00113 else if (fitter->useRes() == MuonResidualsFitter::k0010) {
00114 fval += -weight * MuonResidualsFitter_logPureGaussian(resslope, resslopepeak, resslopesigma);
00115 }
00116 }
00117 else if (fitter->residualsModel() == MuonResidualsFitter::kPowerLawTails) {
00118 fval += -weight * MuonResidualsFitter_logPowerLawTails(residual, residpeak, residsigma, residgamma);
00119 fval += -weight * MuonResidualsFitter_logPowerLawTails(resslope, resslopepeak, resslopesigma, resslopegamma);
00120 }
00121 else if (fitter->residualsModel() == MuonResidualsFitter::kROOTVoigt) {
00122 fval += -weight * MuonResidualsFitter_logROOTVoigt(residual, residpeak, residsigma, residgamma);
00123 fval += -weight * MuonResidualsFitter_logROOTVoigt(resslope, resslopepeak, resslopesigma, resslopegamma);
00124 }
00125 else if (fitter->residualsModel() == MuonResidualsFitter::kGaussPowerTails) {
00126 fval += -weight * MuonResidualsFitter_logGaussPowerTails(residual, residpeak, residsigma);
00127 fval += -weight * MuonResidualsFitter_logGaussPowerTails(resslope, resslopepeak, resslopesigma);
00128 }
00129 else { assert(false); }
00130 }
00131 }
00132 }
00133
00134
00135 void MuonResiduals5DOFFitter::correctBField()
00136 {
00137 MuonResidualsFitter::correctBField(kPt, kCharge);
00138 }
00139
00140
00141 void MuonResiduals5DOFFitter::inform(TMinuit *tMinuit)
00142 {
00143 minuit = tMinuit;
00144 }
00145
00146
00147 double MuonResiduals5DOFFitter::sumofweights()
00148 {
00149 sum_of_weights = 0.;
00150 number_of_hits = 0.;
00151 weight_alignment = m_weightAlignment;
00152 for (std::vector<double*>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
00153 if (m_weightAlignment) {
00154 double redchi2 = (*resiter)[MuonResiduals5DOFFitter::kRedChi2];
00155 if (TMath::Prob(redchi2*8, 8) < 0.99) {
00156 sum_of_weights += 1./redchi2;
00157 number_of_hits += 1.;
00158 }
00159 }
00160 else {
00161 sum_of_weights += 1.;
00162 number_of_hits += 1.;
00163 }
00164 }
00165 return sum_of_weights;
00166 }
00167
00168
00169 bool MuonResiduals5DOFFitter::fit(Alignable *ali)
00170 {
00171 initialize_table();
00172 sumofweights();
00173
00174 double res_std = 0.5;
00175 double resslope_std = 0.005;
00176
00177 int nums[10] = {kAlignX, kAlignZ, kAlignPhiX, kAlignPhiY, kAlignPhiZ, kResidSigma, kResSlopeSigma, kAlpha, kResidGamma, kResSlopeGamma};
00178 std::string names[10] = {"AlignX","AlignZ","AlignPhiX","AlignPhiY","AlignPhiZ", "ResidSigma","ResSlopeSigma", "Alpha", "ResidGamma","ResSlopeGamma"};
00179 double starts[10] = {0., 0., 0., 0., 0., res_std, resslope_std, 0., 0.1*res_std, 0.1*resslope_std};
00180 double steps[10] = {0.1, 0.1, 0.001, 0.001, 0.001, 0.001*res_std, 0.001*resslope_std, 0.001, 0.01*res_std, 0.01*resslope_std};
00181 double lows[10] = {0., 0., 0., 0., 0., 0., 0., -1., 0., 0.};
00182 double highs[10] = {0., 0., 0., 0., 0., 10., 0.1, 1., 0., 0.};
00183
00184 std::vector<int> num(nums, nums+5);
00185 std::vector<std::string> name(names, names+5);
00186 std::vector<double> start(starts, starts+5);
00187 std::vector<double> step(steps, steps+5);
00188 std::vector<double> low(lows, lows+5);
00189 std::vector<double> high(highs, highs+5);
00190
00191 bool add_alpha = ( residualsModel() == kPureGaussian2D);
00192 bool add_gamma = ( residualsModel() == kROOTVoigt || residualsModel() == kPowerLawTails);
00193
00194 int idx[4], ni = 0;
00195 if (useRes() == k1111 || useRes() == k1110 || useRes() == k1010) {
00196 for(ni=0; ni<2; ni++) idx[ni] = ni+5;
00197 if (add_alpha) idx[ni++] = 7;
00198 else if (add_gamma) for(; ni<4; ni++) idx[ni] = ni+6;
00199 if (!add_alpha) fix(kAlpha);
00200 }
00201 else if (useRes() == k1100) {
00202 idx[ni++] = 5;
00203 if (add_gamma) idx[ni++] = 8;
00204 fix(kResSlopeSigma);
00205 fix(kAlpha);
00206 }
00207 else if (useRes() == k0010) {
00208 idx[ni++] = 6;
00209 if (add_gamma) idx[ni++] = 9;
00210 fix(kResidSigma);
00211 fix(kAlpha);
00212 }
00213 for (int i=0; i<ni; i++){
00214 num.push_back(nums[idx[i]]);
00215 name.push_back(names[idx[i]]);
00216 start.push_back(starts[idx[i]]);
00217 step.push_back(steps[idx[i]]);
00218 low.push_back(lows[idx[i]]);
00219 high.push_back(highs[idx[i]]);
00220 }
00221
00222 return dofit(&MuonResiduals5DOFFitter_FCN, num, name, start, step, low, high);
00223 }
00224
00225
00226 double MuonResiduals5DOFFitter::plot(std::string name, TFileDirectory *dir, Alignable *ali)
00227 {
00228 sumofweights();
00229
00230 double mean_residual = 0., mean_resslope = 0.;
00231 double mean_trackx = 0., mean_tracky = 0., mean_trackdxdz = 0., mean_trackdydz = 0.;
00232 double sum_w = 0.;
00233
00234 for (std::vector<double*>::const_iterator rit = residuals_begin(); rit != residuals_end(); ++rit)
00235 {
00236 const double redchi2 = (*rit)[kRedChi2];
00237 double weight = 1./redchi2;
00238 if (!m_weightAlignment) weight = 1.;
00239
00240 if (!m_weightAlignment || TMath::Prob(redchi2*6, 6) < 0.99)
00241 {
00242 double factor_w = 1./(sum_w + weight);
00243 mean_residual = factor_w * (sum_w * mean_residual + weight * (*rit)[kResid]);
00244 mean_resslope = factor_w * (sum_w * mean_resslope + weight * (*rit)[kResSlope]);
00245 mean_trackx = factor_w * (sum_w * mean_trackx + weight * (*rit)[kPositionX]);
00246 mean_tracky = factor_w * (sum_w * mean_tracky + weight * (*rit)[kPositionY]);
00247 mean_trackdxdz = factor_w * (sum_w * mean_trackdxdz + weight * (*rit)[kAngleX]);
00248 mean_trackdydz = factor_w * (sum_w * mean_trackdydz + weight * (*rit)[kAngleY]);
00249 sum_w += weight;
00250 }
00251 }
00252
00253 std::string name_residual, name_resslope, name_residual_raw, name_resslope_raw, name_residual_cut, name_alpha;
00254 std::string name_residual_trackx, name_resslope_trackx;
00255 std::string name_residual_tracky, name_resslope_tracky;
00256 std::string name_residual_trackdxdz, name_resslope_trackdxdz;
00257 std::string name_residual_trackdydz, name_resslope_trackdydz;
00258
00259 name_residual = name + "_residual";
00260 name_resslope = name + "_resslope";
00261 name_residual_raw = name + "_residual_raw";
00262 name_resslope_raw = name + "_resslope_raw";
00263 name_residual_cut = name + "_residual_cut";
00264 name_alpha = name + "_alpha";
00265 name_residual_trackx = name + "_residual_trackx";
00266 name_resslope_trackx = name + "_resslope_trackx";
00267 name_residual_tracky = name + "_residual_tracky";
00268 name_resslope_tracky = name + "_resslope_tracky";
00269 name_residual_trackdxdz = name + "_residual_trackdxdz";
00270 name_resslope_trackdxdz = name + "_resslope_trackdxdz";
00271 name_residual_trackdydz = name + "_residual_trackdydz";
00272 name_resslope_trackdydz = name + "_resslope_trackdydz";
00273
00274 double width = ali->surface().width();
00275 double length = ali->surface().length();
00276 int bins_residual = 150, bins_resslope = 100;
00277 double min_residual = -75., max_residual = 75.;
00278 double min_resslope = -50., max_resslope = 50.;
00279 double min_trackx = -width/2., max_trackx = width/2.;
00280 double min_tracky = -length/2., max_tracky = length/2.;
00281 double min_trackdxdz = -1.5, max_trackdxdz = 1.5;
00282 double min_trackdydz = -1.5, max_trackdydz = 1.5;
00283
00284 TH1F *hist_residual = dir->make<TH1F>(name_residual.c_str(), "", bins_residual, min_residual, max_residual);
00285 TH1F *hist_resslope = dir->make<TH1F>(name_resslope.c_str(), "", bins_resslope, min_resslope, max_resslope);
00286 TH1F *hist_residual_raw = dir->make<TH1F>(name_residual_raw.c_str(), "", bins_residual, min_residual, max_residual);
00287 TH1F *hist_resslope_raw = dir->make<TH1F>(name_resslope_raw.c_str(), "", bins_resslope, min_resslope, max_resslope);
00288 TH1F *hist_residual_cut = dir->make<TH1F>(name_residual_cut.c_str(), "", bins_residual, min_residual, max_residual);
00289 TH2F *hist_alpha = dir->make<TH2F>(name_alpha.c_str(), "", 50, min_resslope, max_resslope, 50, -50., 50.);
00290
00291 TProfile *hist_residual_trackx = dir->make<TProfile>(name_residual_trackx.c_str(), "", 50, min_trackx, max_trackx, min_residual, max_residual);
00292 TProfile *hist_resslope_trackx = dir->make<TProfile>(name_resslope_trackx.c_str(), "", 50, min_trackx, max_trackx, min_resslope, max_resslope);
00293 TProfile *hist_residual_tracky = dir->make<TProfile>(name_residual_tracky.c_str(), "", 50, min_tracky, max_tracky, min_residual, max_residual);
00294 TProfile *hist_resslope_tracky = dir->make<TProfile>(name_resslope_tracky.c_str(), "", 50, min_tracky, max_tracky, min_resslope, max_resslope);
00295 TProfile *hist_residual_trackdxdz = dir->make<TProfile>(name_residual_trackdxdz.c_str(), "", 250, min_trackdxdz, max_trackdxdz, min_residual, max_residual);
00296 TProfile *hist_resslope_trackdxdz = dir->make<TProfile>(name_resslope_trackdxdz.c_str(), "", 250, min_trackdxdz, max_trackdxdz, min_resslope, max_resslope);
00297 TProfile *hist_residual_trackdydz = dir->make<TProfile>(name_residual_trackdydz.c_str(), "", 250, min_trackdydz, max_trackdydz, min_residual, max_residual);
00298 TProfile *hist_resslope_trackdydz = dir->make<TProfile>(name_resslope_trackdydz.c_str(), "", 250, min_trackdydz, max_trackdydz, min_resslope, max_resslope);
00299
00300 hist_residual_trackx->SetAxisRange(-10., 10., "Y");
00301 hist_resslope_trackx->SetAxisRange(-10., 10., "Y");
00302 hist_residual_tracky->SetAxisRange(-10., 10., "Y");
00303 hist_resslope_tracky->SetAxisRange(-10., 10., "Y");
00304 hist_residual_trackdxdz->SetAxisRange(-10., 10., "Y");
00305 hist_resslope_trackdxdz->SetAxisRange(-10., 10., "Y");
00306 hist_residual_trackdydz->SetAxisRange(-10., 10., "Y");
00307 hist_resslope_trackdydz->SetAxisRange(-10., 10., "Y");
00308
00309 name_residual += "_fit";
00310 name_resslope += "_fit";
00311 name_alpha += "_fit";
00312 name_residual_trackx += "_fit";
00313 name_resslope_trackx += "_fit";
00314 name_residual_tracky += "_fit";
00315 name_resslope_tracky += "_fit";
00316 name_residual_trackdxdz += "_fit";
00317 name_resslope_trackdxdz += "_fit";
00318 name_residual_trackdydz += "_fit";
00319 name_resslope_trackdydz += "_fit";
00320
00321 TF1 *fit_residual = NULL;
00322 TF1 *fit_resslope = NULL;
00323 if (residualsModel() == kPureGaussian || residualsModel() == kPureGaussian2D) {
00324 fit_residual = new TF1(name_residual.c_str(), MuonResidualsFitter_pureGaussian_TF1, min_residual, max_residual, 3);
00325 fit_residual->SetParameters(sum_of_weights * (max_residual - min_residual)/bins_residual, 10.*value(kAlignX), 10.*value(kResidSigma));
00326 const double er_res[3] = {0., 10.*errorerror(kAlignX), 10.*errorerror(kResidSigma)};
00327 fit_residual->SetParErrors(er_res);
00328 fit_resslope = new TF1(name_resslope.c_str(), MuonResidualsFitter_pureGaussian_TF1, min_resslope, max_resslope, 3);
00329 fit_resslope->SetParameters(sum_of_weights * (max_resslope - min_resslope)/bins_resslope, 1000.*value(kAlignPhiY), 1000.*value(kResSlopeSigma));
00330 const double er_resslope[3] = {0., 1000.*errorerror(kAlignPhiY), 1000.*errorerror(kResSlopeSigma)};
00331 fit_resslope->SetParErrors(er_resslope);
00332 }
00333 else if (residualsModel() == kPowerLawTails) {
00334 fit_residual = new TF1(name_residual.c_str(), MuonResidualsFitter_powerLawTails_TF1, min_residual, max_residual, 4);
00335 fit_residual->SetParameters(sum_of_weights * (max_residual - min_residual)/bins_residual, 10.*value(kAlignX), 10.*value(kResidSigma), 10.*value(kResidGamma));
00336 fit_resslope = new TF1(name_resslope.c_str(), MuonResidualsFitter_powerLawTails_TF1, min_resslope, max_resslope, 4);
00337 fit_resslope->SetParameters(sum_of_weights * (max_resslope - min_resslope)/bins_resslope, 1000.*value(kAlignPhiY), 1000.*value(kResSlopeSigma), 1000.*value(kResSlopeGamma));
00338 }
00339 else if (residualsModel() == kROOTVoigt) {
00340 fit_residual = new TF1(name_residual.c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_residual, max_residual, 4);
00341 fit_residual->SetParameters(sum_of_weights * (max_residual - min_residual)/bins_residual, 10.*value(kAlignX), 10.*value(kResidSigma), 10.*value(kResidGamma));
00342 fit_resslope = new TF1(name_resslope.c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_resslope, max_resslope, 4);
00343 fit_resslope->SetParameters(sum_of_weights * (max_resslope - min_resslope)/bins_resslope, 1000.*value(kAlignPhiY), 1000.*value(kResSlopeSigma), 1000.*value(kResSlopeGamma));
00344 }
00345 else if (residualsModel() == kGaussPowerTails) {
00346 fit_residual = new TF1(name_residual.c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_residual, max_residual, 3);
00347 fit_residual->SetParameters(sum_of_weights * (max_residual - min_residual)/bins_residual, 10.*value(kAlignX), 10.*value(kResidSigma));
00348 fit_resslope = new TF1(name_resslope.c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_resslope, max_resslope, 3);
00349 fit_resslope->SetParameters(sum_of_weights * (max_resslope - min_resslope)/bins_resslope, 1000.*value(kAlignPhiY), 1000.*value(kResSlopeSigma));
00350 }
00351 else { assert(false); }
00352
00353 fit_residual->SetLineColor(2); fit_residual->SetLineWidth(2); fit_residual->Write();
00354 fit_resslope->SetLineColor(2); fit_resslope->SetLineWidth(2); fit_resslope->Write();
00355
00356 TF1 *fit_alpha = new TF1(name_alpha.c_str(), "[0] + x*[1]", min_resslope, max_resslope);
00357 double a = 10.*value(kAlignX), b = 10.*value(kAlpha)/1000.;
00358 if (residualsModel() == kPureGaussian2D)
00359 {
00360 double sx = 10.*value(kResidSigma), sy = 1000.*value(kResSlopeSigma), r = value(kAlpha);
00361 a = mean_residual;
00362 b = 0.;
00363 if ( sx != 0. )
00364 {
00365 b = 1./(sy/sx*r);
00366 a = - b * mean_resslope;
00367 }
00368 }
00369 fit_alpha->SetParameters(a, b);
00370 fit_alpha->SetLineColor(2); fit_alpha->SetLineWidth(2); fit_alpha->Write();
00371
00372 TProfile *fit_residual_trackx = dir->make<TProfile>(name_residual_trackx.c_str(), "", 50, min_trackx, max_trackx);
00373 TProfile *fit_resslope_trackx = dir->make<TProfile>(name_resslope_trackx.c_str(), "", 50, min_trackx, max_trackx);
00374 TProfile *fit_residual_tracky = dir->make<TProfile>(name_residual_tracky.c_str(), "", 50, min_tracky, max_tracky);
00375 TProfile *fit_resslope_tracky = dir->make<TProfile>(name_resslope_tracky.c_str(), "", 50, min_tracky, max_tracky);
00376 TProfile *fit_residual_trackdxdz = dir->make<TProfile>(name_residual_trackdxdz.c_str(), "", 250, min_trackdxdz, max_trackdxdz);
00377 TProfile *fit_resslope_trackdxdz = dir->make<TProfile>(name_resslope_trackdxdz.c_str(), "", 250, min_trackdxdz, max_trackdxdz);
00378 TProfile *fit_residual_trackdydz = dir->make<TProfile>(name_residual_trackdydz.c_str(), "", 250, min_trackdydz, max_trackdydz);
00379 TProfile *fit_resslope_trackdydz = dir->make<TProfile>(name_resslope_trackdydz.c_str(), "", 250, min_trackdydz, max_trackdydz);
00380
00381 fit_residual_trackx->SetLineColor(2); fit_residual_trackx->SetLineWidth(2);
00382 fit_resslope_trackx->SetLineColor(2); fit_resslope_trackx->SetLineWidth(2);
00383 fit_residual_tracky->SetLineColor(2); fit_residual_tracky->SetLineWidth(2);
00384 fit_resslope_tracky->SetLineColor(2); fit_resslope_tracky->SetLineWidth(2);
00385 fit_residual_trackdxdz->SetLineColor(2); fit_residual_trackdxdz->SetLineWidth(2);
00386 fit_resslope_trackdxdz->SetLineColor(2); fit_resslope_trackdxdz->SetLineWidth(2);
00387 fit_residual_trackdydz->SetLineColor(2); fit_residual_trackdydz->SetLineWidth(2);
00388 fit_resslope_trackdydz->SetLineColor(2); fit_resslope_trackdydz->SetLineWidth(2);
00389
00390 name_residual_trackx += "line";
00391 name_resslope_trackx += "line";
00392 name_residual_tracky += "line";
00393 name_resslope_tracky += "line";
00394 name_residual_trackdxdz += "line";
00395 name_resslope_trackdxdz += "line";
00396 name_residual_trackdydz += "line";
00397 name_resslope_trackdydz += "line";
00398
00399 TF1 *fitline_residual_trackx = new TF1(name_residual_trackx.c_str(), residual_x_trackx_TF1, min_trackx, max_trackx, 12);
00400 TF1 *fitline_resslope_trackx = new TF1(name_resslope_trackx.c_str(), residual_dxdz_trackx_TF1, min_trackx, max_trackx, 12);
00401 TF1 *fitline_residual_tracky = new TF1(name_residual_tracky.c_str(), residual_x_tracky_TF1, min_tracky, max_tracky, 12);
00402 TF1 *fitline_resslope_tracky = new TF1(name_resslope_tracky.c_str(), residual_dxdz_tracky_TF1, min_tracky, max_tracky, 12);
00403 TF1 *fitline_residual_trackdxdz = new TF1(name_residual_trackdxdz.c_str(), residual_x_trackdxdz_TF1, min_trackdxdz, max_trackdxdz, 12);
00404 TF1 *fitline_resslope_trackdxdz = new TF1(name_resslope_trackdxdz.c_str(), residual_dxdz_trackdxdz_TF1, min_trackdxdz, max_trackdxdz, 12);
00405 TF1 *fitline_residual_trackdydz = new TF1(name_residual_trackdydz.c_str(), residual_x_trackdydz_TF1, min_trackdydz, max_trackdydz, 12);
00406 TF1 *fitline_resslope_trackdydz = new TF1(name_resslope_trackdydz.c_str(), residual_dxdz_trackdydz_TF1, min_trackdydz, max_trackdydz, 12);
00407
00408 std::vector<TF1*> fitlines;
00409 fitlines.push_back(fitline_residual_trackx);
00410 fitlines.push_back(fitline_resslope_trackx);
00411 fitlines.push_back(fitline_residual_tracky);
00412 fitlines.push_back(fitline_resslope_tracky);
00413 fitlines.push_back(fitline_residual_trackdxdz);
00414 fitlines.push_back(fitline_resslope_trackdxdz);
00415 fitlines.push_back(fitline_residual_trackdydz);
00416 fitlines.push_back(fitline_resslope_trackdydz);
00417
00418 double fitparameters[12] = {value(kAlignX), 0., value(kAlignZ), value(kAlignPhiX), value(kAlignPhiY), value(kAlignPhiZ),
00419 mean_trackx, mean_tracky, mean_trackdxdz, mean_trackdydz, value(kAlpha), mean_resslope};
00420 if (residualsModel() == kPureGaussian2D) fitparameters[10] = 0.;
00421
00422 for(std::vector<TF1*>::const_iterator itr = fitlines.begin(); itr != fitlines.end(); itr++)
00423 {
00424 (*itr)->SetParameters(fitparameters);
00425 (*itr)->SetLineColor(2);
00426 (*itr)->SetLineWidth(2);
00427 (*itr)->Write();
00428 }
00429
00430 for (std::vector<double*>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
00431 const double resid = (*resiter)[kResid];
00432 const double resslope = (*resiter)[kResSlope];
00433 const double positionX = (*resiter)[kPositionX];
00434 const double positionY = (*resiter)[kPositionY];
00435 const double angleX = (*resiter)[kAngleX];
00436 const double angleY = (*resiter)[kAngleY];
00437 const double redchi2 = (*resiter)[kRedChi2];
00438 double weight = 1./redchi2;
00439 if (!m_weightAlignment) weight = 1.;
00440
00441 if (!m_weightAlignment || TMath::Prob(redchi2*8, 8) < 0.99) {
00442 hist_alpha->Fill(1000.*resslope, 10.*resid);
00443
00444 double coeff = value(kAlpha);
00445 if (residualsModel() == kPureGaussian || residualsModel() == kPureGaussian2D) coeff = 0.;
00446 double geom_resid = residual_x(value(kAlignX), value(kAlignZ), value(kAlignPhiX), value(kAlignPhiY), value(kAlignPhiZ), positionX, positionY, angleX, angleY, coeff, resslope);
00447 hist_residual->Fill(10.*(resid - geom_resid + value(kAlignX)), weight);
00448 hist_residual_trackx->Fill(positionX, 10.*resid, weight);
00449 hist_residual_tracky->Fill(positionY, 10.*resid, weight);
00450 hist_residual_trackdxdz->Fill(angleX, 10.*resid, weight);
00451 hist_residual_trackdydz->Fill(angleY, 10.*resid, weight);
00452 fit_residual_trackx->Fill(positionX, 10.*geom_resid, weight);
00453 fit_residual_tracky->Fill(positionY, 10.*geom_resid, weight);
00454 fit_residual_trackdxdz->Fill(angleX, 10.*geom_resid, weight);
00455 fit_residual_trackdydz->Fill(angleY, 10.*geom_resid, weight);
00456
00457 double geom_resslope = residual_dxdz(value(kAlignX), value(kAlignZ), value(kAlignPhiX), value(kAlignPhiY), value(kAlignPhiZ), positionX, positionY, angleX, angleY);
00458 hist_resslope->Fill(1000.*(resslope - geom_resslope + value(kAlignPhiY)), weight);
00459 hist_resslope_trackx->Fill(positionX, 1000.*resslope, weight);
00460 hist_resslope_tracky->Fill(positionY, 1000.*resslope, weight);
00461 hist_resslope_trackdxdz->Fill(angleX, 1000.*resslope, weight);
00462 hist_resslope_trackdydz->Fill(angleY, 1000.*resslope, weight);
00463 fit_resslope_trackx->Fill(positionX, 1000.*geom_resslope, weight);
00464 fit_resslope_tracky->Fill(positionY, 1000.*geom_resslope, weight);
00465 fit_resslope_trackdxdz->Fill(angleX, 1000.*geom_resslope, weight);
00466 fit_resslope_trackdydz->Fill(angleY, 1000.*geom_resslope, weight);
00467 }
00468
00469 hist_residual_raw->Fill(10.*resid);
00470 hist_resslope_raw->Fill(1000.*resslope);
00471 if (fabs(resslope) < 0.005) hist_residual_cut->Fill(10.*resid);
00472 }
00473
00474 double chi2 = 0.;
00475 double ndof = 0.;
00476 for (int i = 1; i <= hist_residual->GetNbinsX(); i++) {
00477 double xi = hist_residual->GetBinCenter(i);
00478 double yi = hist_residual->GetBinContent(i);
00479 double yerri = hist_residual->GetBinError(i);
00480 double yth = fit_residual->Eval(xi);
00481 if (yerri > 0.) {
00482 chi2 += pow((yth - yi)/yerri, 2);
00483 ndof += 1.;
00484 }
00485 }
00486 for (int i = 1; i <= hist_resslope->GetNbinsX(); i++) {
00487 double xi = hist_resslope->GetBinCenter(i);
00488 double yi = hist_resslope->GetBinContent(i);
00489 double yerri = hist_resslope->GetBinError(i);
00490 double yth = fit_resslope->Eval(xi);
00491 if (yerri > 0.) {
00492 chi2 += pow((yth - yi)/yerri, 2);
00493 ndof += 1.;
00494 }
00495 }
00496 ndof -= npar();
00497
00498 return (ndof > 0. ? chi2 / ndof : -1.);
00499 }
00500
00501
00502 TTree * MuonResiduals5DOFFitter::readNtuple(std::string fname, unsigned int wheel, unsigned int station, unsigned int sector, unsigned int preselected)
00503 {
00504 TFile *f = new TFile(fname.c_str());
00505 TTree *t = (TTree*)f->Get("mual_ttree");
00506
00507
00508 TFile *tmpf = new TFile("small_tree_fit.root","recreate");
00509 assert(tmpf!=0);
00510
00511
00512 TTree *tt = t->CopyTree(Form("is_dt && ring_wheel==%d && station==%d && sector==%d && select==%d", wheel, station, sector, (bool)preselected));
00513
00514 MuonAlignmentTreeRow r;
00515 tt->SetBranchAddress("res_x", &r.res_x);
00516 tt->SetBranchAddress("res_slope_x", &r.res_slope_x);
00517 tt->SetBranchAddress("pos_x", &r.pos_x);
00518 tt->SetBranchAddress("pos_y", &r.pos_y);
00519 tt->SetBranchAddress("angle_x", &r.angle_x);
00520 tt->SetBranchAddress("angle_y", &r.angle_y);
00521 tt->SetBranchAddress("pz", &r.pz);
00522 tt->SetBranchAddress("pt", &r.pt);
00523 tt->SetBranchAddress("q", &r.q);
00524
00525 Long64_t nentries = tt->GetEntries();
00526 for (Long64_t i=0;i<nentries; i++)
00527 {
00528 tt->GetEntry(i);
00529 double *rdata = new double[MuonResiduals5DOFFitter::kNData];
00530 rdata[kResid] = r.res_x;
00531 rdata[kResSlope] = r.res_slope_x;
00532 rdata[kPositionX] = r.pos_x;
00533 rdata[kPositionY] = r.pos_y;
00534 rdata[kAngleX] = r.angle_x;
00535 rdata[kAngleY] = r.angle_y;
00536 rdata[kRedChi2] = 0.1;
00537 rdata[kPz] = r.pz;
00538 rdata[kPt] = r.pt;
00539 rdata[kCharge] = r.q;
00540 fill(rdata);
00541 }
00542 delete f;
00543
00544 return tt;
00545 }
1.7.3