00001
00002
00003 #ifdef STANDALONE_FITTER
00004 #include "MuonResiduals6DOFFitter.h"
00005 #else
00006 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResiduals6DOFFitter.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_y, 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 alphax, double residual_dxdz)
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 + residual_dxdz * alphax;
00033 }
00034
00035 double residual_y(double delta_x, double delta_y, 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 double alphay, double residual_dydz)
00039 {
00040 return delta_y
00041 - track_dydz * delta_z
00042 - track_y * track_dydz * delta_phix
00043 + track_x * track_dydz * delta_phiy
00044 + track_x * delta_phiz
00045 + residual_dydz * alphay;
00046 }
00047
00048 double residual_dxdz(double delta_x, double delta_y, double delta_z,
00049 double delta_phix, double delta_phiy, double delta_phiz,
00050 double track_x, double track_y, double track_dxdz, double track_dydz)
00051 {
00052 return -track_dxdz * track_dydz * delta_phix
00053 + (1. + track_dxdz * track_dxdz) * delta_phiy
00054 - track_dydz * delta_phiz;
00055 }
00056
00057 double residual_dydz(double delta_x, double delta_y, double delta_z,
00058 double delta_phix, double delta_phiy, double delta_phiz,
00059 double track_x, double track_y, double track_dxdz, double track_dydz)
00060 {
00061 return -(1. + track_dydz * track_dydz) * delta_phix
00062 + track_dxdz * track_dydz * delta_phiy
00063 + track_dxdz * delta_phiz;
00064 }
00065
00066 Double_t residual_x_trackx_TF1(Double_t *xx, Double_t *p) { return 10.*residual_x(p[0], p[1], p[2], p[3], p[4], p[5], xx[0], p[7], p[8], p[9], p[10], p[11]); }
00067 Double_t residual_x_tracky_TF1(Double_t *xx, Double_t *p) { return 10.*residual_x(p[0], p[1], p[2], p[3], p[4], p[5], p[6], xx[0], p[8], p[9], p[10], p[11]); }
00068 Double_t residual_x_trackdxdz_TF1(Double_t *xx, Double_t *p) { return 10.*residual_x(p[0], p[1], p[2], p[3], p[4], p[5], p[6], p[7], xx[0], p[9], p[10], p[11]); }
00069 Double_t residual_x_trackdydz_TF1(Double_t *xx, Double_t *p) { return 10.*residual_x(p[0], p[1], p[2], p[3], p[4], p[5], p[6], p[7], p[8], xx[0], p[10], p[11]); }
00070
00071 Double_t residual_y_trackx_TF1(Double_t *xx, Double_t *p) { return 10.*residual_y(p[0], p[1], p[2], p[3], p[4], p[5], xx[0], p[7], p[8], p[9], p[12], p[13]); }
00072 Double_t residual_y_tracky_TF1(Double_t *xx, Double_t *p) { return 10.*residual_y(p[0], p[1], p[2], p[3], p[4], p[5], p[6], xx[0], p[8], p[9], p[12], p[13]); }
00073 Double_t residual_y_trackdxdz_TF1(Double_t *xx, Double_t *p) { return 10.*residual_y(p[0], p[1], p[2], p[3], p[4], p[5], p[6], p[7], xx[0], p[9], p[12], p[13]); }
00074 Double_t residual_y_trackdydz_TF1(Double_t *xx, Double_t *p) { return 10.*residual_y(p[0], p[1], p[2], p[3], p[4], p[5], p[6], p[7], p[8], xx[0], p[12], p[13]); }
00075
00076 Double_t residual_dxdz_trackx_TF1(Double_t *xx, Double_t *p) { return 1000.*residual_dxdz(p[0], p[1], p[2], p[3], p[4], p[5], xx[0], p[7], p[8], p[9]); }
00077 Double_t residual_dxdz_tracky_TF1(Double_t *xx, Double_t *p) { return 1000.*residual_dxdz(p[0], p[1], p[2], p[3], p[4], p[5], p[6], xx[0], p[8], p[9]); }
00078 Double_t residual_dxdz_trackdxdz_TF1(Double_t *xx, Double_t *p) { return 1000.*residual_dxdz(p[0], p[1], p[2], p[3], p[4], p[5], p[6], p[7], xx[0], p[9]); }
00079 Double_t residual_dxdz_trackdydz_TF1(Double_t *xx, Double_t *p) { return 1000.*residual_dxdz(p[0], p[1], p[2], p[3], p[4], p[5], p[6], p[7], p[8], xx[0]); }
00080
00081 Double_t residual_dydz_trackx_TF1(Double_t *xx, Double_t *p) { return 1000.*residual_dydz(p[0], p[1], p[2], p[3], p[4], p[5], xx[0], p[7], p[8], p[9]); }
00082 Double_t residual_dydz_tracky_TF1(Double_t *xx, Double_t *p) { return 1000.*residual_dydz(p[0], p[1], p[2], p[3], p[4], p[5], p[6], xx[0], p[8], p[9]); }
00083 Double_t residual_dydz_trackdxdz_TF1(Double_t *xx, Double_t *p) { return 1000.*residual_dydz(p[0], p[1], p[2], p[3], p[4], p[5], p[6], p[7], xx[0], p[9]); }
00084 Double_t residual_dydz_trackdydz_TF1(Double_t *xx, Double_t *p) { return 1000.*residual_dydz(p[0], p[1], p[2], p[3], p[4], p[5], p[6], p[7], p[8], xx[0]); }
00085 }
00086
00087
00088 void MuonResiduals6DOFFitter_FCN(int &npar, double *gin, double &fval, double *par, int iflag)
00089 {
00090 MuonResidualsFitterFitInfo *fitinfo = (MuonResidualsFitterFitInfo*)(minuit->GetObjectFit());
00091 MuonResidualsFitter *fitter = fitinfo->fitter();
00092
00093 fval = 0.;
00094 for (std::vector<double*>::const_iterator resiter = fitter->residuals_begin(); resiter != fitter->residuals_end(); ++resiter)
00095 {
00096 const double residX = (*resiter)[MuonResiduals6DOFFitter::kResidX];
00097 const double residY = (*resiter)[MuonResiduals6DOFFitter::kResidY];
00098 const double resslopeX = (*resiter)[MuonResiduals6DOFFitter::kResSlopeX];
00099 const double resslopeY = (*resiter)[MuonResiduals6DOFFitter::kResSlopeY];
00100 const double positionX = (*resiter)[MuonResiduals6DOFFitter::kPositionX];
00101 const double positionY = (*resiter)[MuonResiduals6DOFFitter::kPositionY];
00102 const double angleX = (*resiter)[MuonResiduals6DOFFitter::kAngleX];
00103 const double angleY = (*resiter)[MuonResiduals6DOFFitter::kAngleY];
00104 const double redchi2 = (*resiter)[MuonResiduals6DOFFitter::kRedChi2];
00105
00106 const double alignx = par[MuonResiduals6DOFFitter::kAlignX];
00107 const double aligny = par[MuonResiduals6DOFFitter::kAlignY];
00108 const double alignz = par[MuonResiduals6DOFFitter::kAlignZ];
00109 const double alignphix = par[MuonResiduals6DOFFitter::kAlignPhiX];
00110 const double alignphiy = par[MuonResiduals6DOFFitter::kAlignPhiY];
00111 const double alignphiz = par[MuonResiduals6DOFFitter::kAlignPhiZ];
00112 const double resXsigma = par[MuonResiduals6DOFFitter::kResidXSigma];
00113 const double resYsigma = par[MuonResiduals6DOFFitter::kResidYSigma];
00114 const double slopeXsigma = par[MuonResiduals6DOFFitter::kResSlopeXSigma];
00115 const double slopeYsigma = par[MuonResiduals6DOFFitter::kResSlopeYSigma];
00116 const double alphax = par[MuonResiduals6DOFFitter::kAlphaX];
00117 const double alphay = par[MuonResiduals6DOFFitter::kAlphaY];
00118 const double resXgamma = par[MuonResiduals6DOFFitter::kResidXGamma];
00119 const double resYgamma = par[MuonResiduals6DOFFitter::kResidYGamma];
00120 const double slopeXgamma = par[MuonResiduals6DOFFitter::kResSlopeXGamma];
00121 const double slopeYgamma = par[MuonResiduals6DOFFitter::kResSlopeYGamma];
00122
00123 double coefX = alphax, coefY = alphay;
00124 if (fitter->residualsModel() == MuonResidualsFitter::kPureGaussian ||
00125 fitter->residualsModel() == MuonResidualsFitter::kPureGaussian2D) coefX = coefY = 0.;
00126 double residXpeak = residual_x(alignx, aligny, alignz, alignphix, alignphiy, alignphiz, positionX, positionY, angleX, angleY, coefX, resslopeX);
00127 double residYpeak = residual_y(alignx, aligny, alignz, alignphix, alignphiy, alignphiz, positionX, positionY, angleX, angleY, coefY, resslopeY);
00128 double slopeXpeak = residual_dxdz(alignx, aligny, alignz, alignphix, alignphiy, alignphiz, positionX, positionY, angleX, angleY);
00129 double slopeYpeak = residual_dydz(alignx, aligny, alignz, alignphix, alignphiy, alignphiz, positionX, positionY, angleX, angleY);
00130
00131 double weight = (1./redchi2) * number_of_hits / sum_of_weights;
00132 if (!weight_alignment) weight = 1.;
00133
00134 if (!weight_alignment || TMath::Prob(redchi2*12, 12) < 0.99)
00135 {
00136 if (fitter->residualsModel() == MuonResidualsFitter::kPureGaussian) {
00137 if (fitter->useRes() == MuonResidualsFitter::k1111) {
00138 fval += -weight * MuonResidualsFitter_logPureGaussian(residX, residXpeak, resXsigma);
00139 fval += -weight * MuonResidualsFitter_logPureGaussian(residY, residYpeak, resYsigma);
00140 fval += -weight * MuonResidualsFitter_logPureGaussian(resslopeX, slopeXpeak, slopeXsigma);
00141 fval += -weight * MuonResidualsFitter_logPureGaussian(resslopeY, slopeYpeak, slopeYsigma);
00142 }
00143 else if (fitter->useRes() == MuonResidualsFitter::k1110) {
00144 fval += -weight * MuonResidualsFitter_logPureGaussian(residX, residXpeak, resXsigma);
00145 fval += -weight * MuonResidualsFitter_logPureGaussian(residY, residYpeak, resYsigma);
00146 fval += -weight * MuonResidualsFitter_logPureGaussian(resslopeX, slopeXpeak, slopeXsigma);
00147 }
00148 else if (fitter->useRes() == MuonResidualsFitter::k1100) {
00149 fval += -weight * MuonResidualsFitter_logPureGaussian(residX, residXpeak, resXsigma);
00150 fval += -weight * MuonResidualsFitter_logPureGaussian(residY, residYpeak, resYsigma);
00151 }
00152 else if (fitter->useRes() == MuonResidualsFitter::k1010) {
00153 fval += -weight * MuonResidualsFitter_logPureGaussian(residX, residXpeak, resXsigma);
00154 fval += -weight * MuonResidualsFitter_logPureGaussian(resslopeX, slopeXpeak, slopeXsigma);
00155 }
00156 else if (fitter->useRes() == MuonResidualsFitter::k0010) {
00157 fval += -weight * MuonResidualsFitter_logPureGaussian(resslopeX, slopeXpeak, slopeXsigma);
00158 }
00159
00160
00161
00162
00163 }
00164 else if (fitter->residualsModel() == MuonResidualsFitter::kPureGaussian2D) {
00165 if (fitter->useRes() == MuonResidualsFitter::k1111) {
00166 fval += -weight * MuonResidualsFitter_logPureGaussian2D(residX, resslopeX, residXpeak, slopeXpeak, resXsigma, slopeXsigma, alphax);
00167 fval += -weight * MuonResidualsFitter_logPureGaussian2D(residY, resslopeY, residYpeak, slopeYpeak, resYsigma, slopeYsigma, alphay);
00168 }
00169 else if (fitter->useRes() == MuonResidualsFitter::k1110) {
00170 fval += -weight * MuonResidualsFitter_logPureGaussian2D(residX, resslopeX, residXpeak, slopeXpeak, resXsigma, slopeXsigma, alphax);
00171 fval += -weight * MuonResidualsFitter_logPureGaussian(residY, residYpeak, resYsigma);
00172 }
00173 else if (fitter->useRes() == MuonResidualsFitter::k1100) {
00174 fval += -weight * MuonResidualsFitter_logPureGaussian(residX, residXpeak, resXsigma);
00175 fval += -weight * MuonResidualsFitter_logPureGaussian(residY, residYpeak, resYsigma);
00176 }
00177 else if (fitter->useRes() == MuonResidualsFitter::k1010) {
00178 fval += -weight * MuonResidualsFitter_logPureGaussian2D(residX, resslopeX, residXpeak, slopeXpeak, resXsigma, slopeXsigma, alphax);
00179 }
00180 else if (fitter->useRes() == MuonResidualsFitter::k0010) {
00181 fval += -weight * MuonResidualsFitter_logPureGaussian(resslopeX, slopeXpeak, slopeXsigma);
00182 }
00183
00184
00185 }
00186 else if (fitter->residualsModel() == MuonResidualsFitter::kPowerLawTails) {
00187 fval += -weight * MuonResidualsFitter_logPowerLawTails(residX, residXpeak, resXsigma, resXgamma);
00188 fval += -weight * MuonResidualsFitter_logPowerLawTails(residY, residYpeak, resYsigma, resYgamma);
00189 fval += -weight * MuonResidualsFitter_logPowerLawTails(resslopeX, slopeXpeak, slopeXsigma, slopeXgamma);
00190 fval += -weight * MuonResidualsFitter_logPowerLawTails(resslopeY, slopeYpeak, slopeYsigma, slopeYgamma);
00191 }
00192 else if (fitter->residualsModel() == MuonResidualsFitter::kROOTVoigt) {
00193 fval += -weight * MuonResidualsFitter_logROOTVoigt(residX, residXpeak, resXsigma, resXgamma);
00194 fval += -weight * MuonResidualsFitter_logROOTVoigt(residY, residYpeak, resYsigma, resYgamma);
00195 fval += -weight * MuonResidualsFitter_logROOTVoigt(resslopeX, slopeXpeak, slopeXsigma, slopeXgamma);
00196 fval += -weight * MuonResidualsFitter_logROOTVoigt(resslopeY, slopeYpeak, slopeYsigma, slopeYgamma);
00197 }
00198 else if (fitter->residualsModel() == MuonResidualsFitter::kGaussPowerTails) {
00199 fval += -weight * MuonResidualsFitter_logGaussPowerTails(residX, residXpeak, resXsigma);
00200 fval += -weight * MuonResidualsFitter_logGaussPowerTails(residY, residYpeak, resYsigma);
00201 fval += -weight * MuonResidualsFitter_logGaussPowerTails(resslopeX, slopeXpeak, slopeXsigma);
00202 fval += -weight * MuonResidualsFitter_logGaussPowerTails(resslopeY, slopeYpeak, slopeYsigma);
00203 }
00204 else { assert(false); }
00205 }
00206 }
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216 }
00217
00218
00219 void MuonResiduals6DOFFitter::inform(TMinuit *tMinuit)
00220 {
00221 minuit = tMinuit;
00222 }
00223
00224
00225 void MuonResiduals6DOFFitter::correctBField()
00226 {
00227 MuonResidualsFitter::correctBField(kPt, kCharge);
00228 }
00229
00230
00231 double MuonResiduals6DOFFitter::sumofweights()
00232 {
00233 sum_of_weights = 0.;
00234 number_of_hits = 0.;
00235 weight_alignment = m_weightAlignment;
00236 for (std::vector<double*>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
00237 if (m_weightAlignment) {
00238 double redchi2 = (*resiter)[kRedChi2];
00239 if (TMath::Prob(redchi2*12, 12) < 0.99) {
00240 sum_of_weights += 1./redchi2;
00241 number_of_hits += 1.;
00242 }
00243 }
00244 else {
00245 sum_of_weights += 1.;
00246 number_of_hits += 1.;
00247 }
00248 }
00249 return sum_of_weights;
00250 }
00251
00252
00253 bool MuonResiduals6DOFFitter::fit(Alignable *ali)
00254 {
00255 initialize_table();
00256 sumofweights();
00257
00258 double resx_std = 0.5;
00259 double resy_std = 3.0;
00260 double resslopex_std = 0.002;
00261 double resslopey_std = 0.005;
00262
00263 int nums[16] = {kAlignX, kAlignY, kAlignZ, kAlignPhiX, kAlignPhiY, kAlignPhiZ, kResidXSigma, kResidYSigma, kResSlopeXSigma, kResSlopeYSigma,
00264 kAlphaX, kAlphaY, kResidXGamma, kResidYGamma, kResSlopeXGamma, kResSlopeYGamma};
00265 std::string names[16] = {"AlignX","AlignY","AlignZ","AlignPhiX","AlignPhiY","AlignPhiZ", "ResidXSigma","ResidYSigma","ResSlopeXSigma","ResSlopeYSigma",
00266 "AlphaX","AlphaY", "ResidXGamma","ResidYGamma","ResSlopeXGamma","ResSlopeYGamma"};
00267 double starts[16] = {0., 0., 0., 0., 0., 0., resx_std, resy_std, resslopex_std, resslopey_std,
00268 0., 0., 0.1*resx_std, 0.1*resy_std, 0.1*resslopex_std, 0.1*resslopey_std};
00269 double steps[16] = {0.1, 0.1, 0.1, 0.001, 0.001, 0.001, 0.001*resx_std, 0.001*resy_std, 0.001*resslopex_std, 0.001*resslopey_std,
00270 0.001, 0.001, 0.01*resx_std, 0.01*resy_std, 0.01*resslopex_std, 0.01*resslopey_std};
00271 double lows[16] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
00272 -1., -1., 0., 0., 0., 0.};
00273 double highs[16] = {0., 0., 0., 0., 0., 0., 10., 10., 0.1, 0.1,
00274 1.,1., 0., 0., 0., 0.};
00275
00276 std::vector<int> num(nums, nums+6);
00277 std::vector<std::string> name(names, names+6);
00278 std::vector<double> start(starts, starts+6);
00279 std::vector<double> step(steps, steps+6);
00280 std::vector<double> low(lows, lows+6);
00281 std::vector<double> high(highs, highs+6);
00282
00283 bool add_alpha = ( residualsModel() == kPureGaussian2D );
00284 bool add_gamma = ( residualsModel() == kROOTVoigt || residualsModel() == kPowerLawTails );
00285
00286 int idx[8], ni = 0;
00287 if (useRes() == k1111) {
00288 for(ni=0; ni<4; ni++) idx[ni] = ni+6;
00289 if (add_alpha) for(; ni<6; ni++) idx[ni] = ni+6;
00290 else if (add_gamma) for(; ni<8; ni++) idx[ni] = ni+8;
00291 if (!add_alpha) fix(kAlphaX);
00292 if (!add_alpha) fix(kAlphaY);
00293 }
00294 else if (useRes() == k1110) {
00295 for(ni=0; ni<3; ni++) idx[ni] = ni+6;
00296 if (add_alpha) idx[ni++] = 10;
00297 else if (add_gamma) for(; ni<6; ni++) idx[ni] = ni+9;
00298 fix(kResSlopeYSigma);
00299 fix(kAlphaY);
00300 if (!add_alpha) fix(kAlphaX);
00301 }
00302 else if (useRes() == k1100) {
00303 for(ni=0; ni<2; ni++) idx[ni] = ni+6;
00304 if (add_gamma) for(; ni<4; ni++) idx[ni] = ni+10;
00305 fix(kResSlopeXSigma);
00306 fix(kResSlopeYSigma);
00307 fix(kAlphaX);
00308 fix(kAlphaY);
00309 }
00310 else if (useRes() == k1010) {
00311 idx[ni++] = 6; idx[ni++] = 8;
00312 if (add_alpha) idx[ni++] = 10;
00313 if (add_gamma) { idx[ni++] = 12; idx[ni++] = 14; }
00314 fix(kResidYSigma);
00315 fix(kResSlopeYSigma);
00316 fix(kAlphaY);
00317 if (!add_alpha) fix(kAlphaX);
00318 }
00319 else if (useRes() == k0010) {
00320 idx[ni++] = 8;
00321 if (add_gamma) idx[ni++] = 14;
00322 fix(kResidXSigma);
00323 fix(kResidYSigma);
00324 fix(kResSlopeYSigma);
00325 fix(kAlphaX);
00326 fix(kAlphaY);
00327 }
00328 for (int i=0; i<ni; i++){
00329 num.push_back(nums[idx[i]]);
00330 name.push_back(names[idx[i]]);
00331 start.push_back(starts[idx[i]]);
00332 step.push_back(steps[idx[i]]);
00333 low.push_back(lows[idx[i]]);
00334 high.push_back(highs[idx[i]]);
00335 }
00336
00337 return dofit(&MuonResiduals6DOFFitter_FCN, num, name, start, step, low, high);
00338 }
00339
00340
00341 double MuonResiduals6DOFFitter::plot(std::string name, TFileDirectory *dir, Alignable *ali)
00342 {
00343 sumofweights();
00344
00345 double mean_residualx = 0., mean_residualy = 0., mean_resslopex = 0., mean_resslopey = 0.;
00346 double mean_trackx = 0., mean_tracky = 0., mean_trackdxdz = 0., mean_trackdydz = 0.;
00347 double sum_w = 0.;
00348
00349 for (std::vector<double*>::const_iterator rit = residuals_begin(); rit != residuals_end(); ++rit)
00350 {
00351 const double redchi2 = (*rit)[kRedChi2];
00352 double weight = 1./redchi2;
00353 if (!m_weightAlignment) weight = 1.;
00354
00355 if (!m_weightAlignment || TMath::Prob(redchi2*12, 12) < 0.99)
00356 {
00357 double factor_w = 1./(sum_w + weight);
00358 mean_residualx = factor_w * (sum_w * mean_residualx + weight * (*rit)[kResidX]);
00359 mean_residualy = factor_w * (sum_w * mean_residualy + weight * (*rit)[kResidY]);
00360 mean_resslopex = factor_w * (sum_w * mean_resslopex + weight * (*rit)[kResSlopeX]);
00361 mean_resslopey = factor_w * (sum_w * mean_resslopey + weight * (*rit)[kResSlopeY]);
00362 mean_trackx = factor_w * (sum_w * mean_trackx + weight * (*rit)[kPositionX]);
00363 mean_tracky = factor_w * (sum_w * mean_tracky + weight * (*rit)[kPositionY]);
00364 mean_trackdxdz = factor_w * (sum_w * mean_trackdxdz + weight * (*rit)[kAngleX]);
00365 mean_trackdydz = factor_w * (sum_w * mean_trackdydz + weight * (*rit)[kAngleY]);
00366 sum_w += weight;
00367 }
00368 }
00369
00370 std::string name_x, name_y, name_dxdz, name_dydz, name_x_raw, name_y_raw, name_dxdz_raw, name_dydz_raw, name_x_cut, name_y_cut, name_alphax, name_alphay;
00371 std::string name_x_trackx, name_y_trackx, name_dxdz_trackx, name_dydz_trackx;
00372 std::string name_x_tracky, name_y_tracky, name_dxdz_tracky, name_dydz_tracky;
00373 std::string name_x_trackdxdz, name_y_trackdxdz, name_dxdz_trackdxdz, name_dydz_trackdxdz;
00374 std::string name_x_trackdydz, name_y_trackdydz, name_dxdz_trackdydz, name_dydz_trackdydz;
00375
00376 name_x = name + "_x";
00377 name_y = name + "_y";
00378 name_dxdz = name + "_dxdz";
00379 name_dydz = name + "_dydz";
00380 name_x_raw = name + "_x_raw";
00381 name_y_raw = name + "_y_raw";
00382 name_dxdz_raw = name + "_dxdz_raw";
00383 name_dydz_raw = name + "_dydz_raw";
00384 name_x_cut = name + "_x_cut";
00385 name_y_cut = name + "_y_cut";
00386 name_alphax = name + "_alphax";
00387 name_alphay = name + "_alphay";
00388 name_x_trackx = name + "_x_trackx";
00389 name_y_trackx = name + "_y_trackx";
00390 name_dxdz_trackx = name + "_dxdz_trackx";
00391 name_dydz_trackx = name + "_dydz_trackx";
00392 name_x_tracky = name + "_x_tracky";
00393 name_y_tracky = name + "_y_tracky";
00394 name_dxdz_tracky = name + "_dxdz_tracky";
00395 name_dydz_tracky = name + "_dydz_tracky";
00396 name_x_trackdxdz = name + "_x_trackdxdz";
00397 name_y_trackdxdz = name + "_y_trackdxdz";
00398 name_dxdz_trackdxdz = name + "_dxdz_trackdxdz";
00399 name_dydz_trackdxdz = name + "_dydz_trackdxdz";
00400 name_x_trackdydz = name + "_x_trackdydz";
00401 name_y_trackdydz = name + "_y_trackdydz";
00402 name_dxdz_trackdydz = name + "_dxdz_trackdydz";
00403 name_dydz_trackdydz = name + "_dydz_trackdydz";
00404
00405 double width = ali->surface().width();
00406 double length = ali->surface().length();
00407 int bins = 200;
00408 double min_x = -100., max_x = 100.;
00409 double min_y = -100., max_y = 100.;
00410 double min_dxdz = -75., max_dxdz = 75.;
00411 double min_dydz = -150., max_dydz = 150.;
00412 double min_trackx = -width/2., max_trackx = width/2.;
00413 double min_tracky = -length/2., max_tracky = length/2.;
00414 double min_trackdxdz = -1.5, max_trackdxdz = 1.5;
00415 double min_trackdydz = -1.5, max_trackdydz = 1.5;
00416
00417 TH1F *hist_x = dir->make<TH1F>(name_x.c_str(), "", bins, min_x, max_x);
00418 TH1F *hist_y = dir->make<TH1F>(name_y.c_str(), "", bins, min_y, max_y);
00419 TH1F *hist_dxdz = dir->make<TH1F>(name_dxdz.c_str(), "", bins, min_dxdz, max_dxdz);
00420 TH1F *hist_dydz = dir->make<TH1F>(name_dydz.c_str(), "", bins, min_dydz, max_dydz);
00421 TH1F *hist_x_raw = dir->make<TH1F>(name_x_raw.c_str(), "", bins, min_x, max_x);
00422 TH1F *hist_y_raw = dir->make<TH1F>(name_y_raw.c_str(), "", bins, min_y, max_y);
00423 TH1F *hist_dxdz_raw = dir->make<TH1F>(name_dxdz_raw.c_str(), "", bins, min_dxdz, max_dxdz);
00424 TH1F *hist_dydz_raw = dir->make<TH1F>(name_dydz_raw.c_str(), "", bins, min_dydz, max_dydz);
00425 TH1F *hist_x_cut = dir->make<TH1F>(name_x_cut.c_str(), "", bins, min_x, max_x);
00426 TH1F *hist_y_cut = dir->make<TH1F>(name_y_cut.c_str(), "", bins, min_y, max_y);
00427 TH2F *hist_alphax = dir->make<TH2F>(name_alphax.c_str(), "", 50, 50, 50, 50, -50., 50.);
00428 TH2F *hist_alphay = dir->make<TH2F>(name_alphay.c_str(), "", 75, 100, 100, 75, -75., 75.);
00429
00430 TProfile *hist_x_trackx = dir->make<TProfile>(name_x_trackx.c_str(), "", 50, min_trackx, max_trackx, min_x, max_x);
00431 TProfile *hist_y_trackx = dir->make<TProfile>(name_y_trackx.c_str(), "", 50, min_trackx, max_trackx, min_y, max_y);
00432 TProfile *hist_dxdz_trackx = dir->make<TProfile>(name_dxdz_trackx.c_str(), "", 50, min_trackx, max_trackx, min_dxdz, max_dxdz);
00433 TProfile *hist_dydz_trackx = dir->make<TProfile>(name_dydz_trackx.c_str(), "", 50, min_trackx, max_trackx, min_dydz, max_dydz);
00434 TProfile *hist_x_tracky = dir->make<TProfile>(name_x_tracky.c_str(), "", 50, min_tracky, max_tracky, min_x, max_x);
00435 TProfile *hist_y_tracky = dir->make<TProfile>(name_y_tracky.c_str(), "", 50, min_tracky, max_tracky, min_y, max_y);
00436 TProfile *hist_dxdz_tracky = dir->make<TProfile>(name_dxdz_tracky.c_str(), "", 50, min_tracky, max_tracky, min_dxdz, max_dxdz);
00437 TProfile *hist_dydz_tracky = dir->make<TProfile>(name_dydz_tracky.c_str(), "", 50, min_tracky, max_tracky, min_dydz, max_dydz);
00438 TProfile *hist_x_trackdxdz = dir->make<TProfile>(name_x_trackdxdz.c_str(), "", 250, min_trackdxdz, max_trackdxdz, min_x, max_x);
00439 TProfile *hist_y_trackdxdz = dir->make<TProfile>(name_y_trackdxdz.c_str(), "", 250, min_trackdxdz, max_trackdxdz, min_y, max_y);
00440 TProfile *hist_dxdz_trackdxdz = dir->make<TProfile>(name_dxdz_trackdxdz.c_str(), "", 250, min_trackdxdz, max_trackdxdz, min_dxdz, max_dxdz);
00441 TProfile *hist_dydz_trackdxdz = dir->make<TProfile>(name_dydz_trackdxdz.c_str(), "", 250, min_trackdxdz, max_trackdxdz, min_dydz, max_dydz);
00442 TProfile *hist_x_trackdydz = dir->make<TProfile>(name_x_trackdydz.c_str(), "", 250, min_trackdydz, max_trackdydz, min_x, max_x);
00443 TProfile *hist_y_trackdydz = dir->make<TProfile>(name_y_trackdydz.c_str(), "", 250, min_trackdydz, max_trackdydz, min_y, max_y);
00444 TProfile *hist_dxdz_trackdydz = dir->make<TProfile>(name_dxdz_trackdydz.c_str(), "", 250, min_trackdydz, max_trackdydz, min_dxdz, max_dxdz);
00445 TProfile *hist_dydz_trackdydz = dir->make<TProfile>(name_dydz_trackdydz.c_str(), "", 250, min_trackdydz, max_trackdydz, min_dydz, max_dydz);
00446
00447 hist_x_trackx->SetAxisRange(-10., 10., "Y");
00448 hist_y_trackx->SetAxisRange(-20., 20., "Y");
00449 hist_dxdz_trackx->SetAxisRange(-10., 10., "Y");
00450 hist_dydz_trackx->SetAxisRange(-20., 20., "Y");
00451 hist_x_tracky->SetAxisRange(-10., 10., "Y");
00452 hist_y_tracky->SetAxisRange(-20., 20., "Y");
00453 hist_dxdz_tracky->SetAxisRange(-10., 10., "Y");
00454 hist_dydz_tracky->SetAxisRange(-20., 20., "Y");
00455 hist_x_trackdxdz->SetAxisRange(-10., 10., "Y");
00456 hist_y_trackdxdz->SetAxisRange(-20., 20., "Y");
00457 hist_dxdz_trackdxdz->SetAxisRange(-10., 10., "Y");
00458 hist_dydz_trackdxdz->SetAxisRange(-20., 20., "Y");
00459 hist_x_trackdydz->SetAxisRange(-10., 10., "Y");
00460 hist_y_trackdydz->SetAxisRange(-20., 20., "Y");
00461 hist_dxdz_trackdydz->SetAxisRange(-10., 10., "Y");
00462 hist_dydz_trackdydz->SetAxisRange(-20., 20., "Y");
00463
00464 name_x += "_fit";
00465 name_y += "_fit";
00466 name_dxdz += "_fit";
00467 name_dydz += "_fit";
00468 name_alphax += "_fit";
00469 name_alphay += "_fit";
00470 name_x_trackx += "_fit";
00471 name_y_trackx += "_fit";
00472 name_dxdz_trackx += "_fit";
00473 name_dydz_trackx += "_fit";
00474 name_x_tracky += "_fit";
00475 name_y_tracky += "_fit";
00476 name_dxdz_tracky += "_fit";
00477 name_dydz_tracky += "_fit";
00478 name_x_trackdxdz += "_fit";
00479 name_y_trackdxdz += "_fit";
00480 name_dxdz_trackdxdz += "_fit";
00481 name_dydz_trackdxdz += "_fit";
00482 name_x_trackdydz += "_fit";
00483 name_y_trackdydz += "_fit";
00484 name_dxdz_trackdydz += "_fit";
00485 name_dydz_trackdydz += "_fit";
00486
00487 TF1 *fit_x = NULL;
00488 TF1 *fit_y = NULL;
00489 TF1 *fit_dxdz = NULL;
00490 TF1 *fit_dydz = NULL;
00491 if (residualsModel() == kPureGaussian || residualsModel() == kPureGaussian2D) {
00492 fit_x = new TF1(name_x.c_str(), MuonResidualsFitter_pureGaussian_TF1, min_x, max_x, 3);
00493 fit_x->SetParameters(sum_of_weights * (max_x - min_x)/bins, 10.*value(kAlignX), 10.*value(kResidXSigma));
00494 const double er_x[3] = {0., 10.*errorerror(kAlignX), 10.*errorerror(kResidXSigma)};
00495 fit_x->SetParErrors(er_x);
00496 fit_y = new TF1(name_y.c_str(), MuonResidualsFitter_pureGaussian_TF1, min_y, max_y, 3);
00497 fit_y->SetParameters(sum_of_weights * (max_y - min_y)/bins, 10.*value(kAlignY), 10.*value(kResidYSigma));
00498 const double er_y[3] = {0., 10.*errorerror(kAlignY), 10.*errorerror(kResidYSigma)};
00499 fit_y->SetParErrors(er_y);
00500 fit_dxdz = new TF1(name_dxdz.c_str(), MuonResidualsFitter_pureGaussian_TF1, min_dxdz, max_dxdz, 3);
00501 fit_dxdz->SetParameters(sum_of_weights * (max_dxdz - min_dxdz)/bins, 1000.*value(kAlignPhiY), 1000.*value(kResSlopeXSigma));
00502 const double er_dxdz[3] = {0., 1000.*errorerror(kAlignPhiY), 1000.*errorerror(kResSlopeXSigma)};
00503 fit_dxdz->SetParErrors(er_dxdz);
00504 fit_dydz = new TF1(name_dydz.c_str(), MuonResidualsFitter_pureGaussian_TF1, min_dydz, max_dydz, 3);
00505 fit_dydz->SetParameters(sum_of_weights * (max_dydz - min_dydz)/bins, -1000.*value(kAlignPhiX), 1000.*value(kResSlopeYSigma));
00506 const double er_dydz[3] = {0., 1000.*errorerror(kAlignPhiX), 1000.*errorerror(kResSlopeYSigma)};
00507 fit_dydz->SetParErrors(er_dydz);
00508 }
00509 else if (residualsModel() == kPowerLawTails) {
00510 fit_x = new TF1(name_x.c_str(), MuonResidualsFitter_powerLawTails_TF1, min_x, max_x, 4);
00511 fit_x->SetParameters(sum_of_weights * (max_x - min_x)/bins, 10.*value(kAlignX), 10.*value(kResidXSigma), 10.*value(kResidXGamma));
00512 fit_y = new TF1(name_y.c_str(), MuonResidualsFitter_powerLawTails_TF1, min_y, max_y, 4);
00513 fit_y->SetParameters(sum_of_weights * (max_y - min_y)/bins, 10.*value(kAlignY), 10.*value(kResidYSigma), 10.*value(kResidYGamma));
00514 fit_dxdz = new TF1(name_dxdz.c_str(), MuonResidualsFitter_powerLawTails_TF1, min_dxdz, max_dxdz, 4);
00515 fit_dxdz->SetParameters(sum_of_weights * (max_dxdz - min_dxdz)/bins, 1000.*value(kAlignPhiY), 1000.*value(kResSlopeXSigma), 1000.*value(kResSlopeXGamma));
00516 fit_dydz = new TF1(name_dydz.c_str(), MuonResidualsFitter_powerLawTails_TF1, min_dydz, max_dydz, 4);
00517 fit_dydz->SetParameters(sum_of_weights * (max_dydz - min_dydz)/bins, -1000.*value(kAlignPhiX), 1000.*value(kResSlopeYSigma), 1000.*value(kResSlopeYGamma));
00518 }
00519 else if (residualsModel() == kROOTVoigt) {
00520 fit_x = new TF1(name_x.c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_x, max_x, 4);
00521 fit_x->SetParameters(sum_of_weights * (max_x - min_x)/bins, 10.*value(kAlignX), 10.*value(kResidXSigma), 10.*value(kResidXGamma));
00522 fit_y = new TF1(name_y.c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_y, max_y, 4);
00523 fit_y->SetParameters(sum_of_weights * (max_y - min_y)/bins, 10.*value(kAlignY), 10.*value(kResidYSigma), 10.*value(kResidYGamma));
00524 fit_dxdz = new TF1(name_dxdz.c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_dxdz, max_dxdz, 4);
00525 fit_dxdz->SetParameters(sum_of_weights * (max_dxdz - min_dxdz)/bins, 1000.*value(kAlignPhiY), 1000.*value(kResSlopeXSigma), 1000.*value(kResSlopeXGamma));
00526 fit_dydz = new TF1(name_dydz.c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_dydz, max_dydz, 4);
00527 fit_dydz->SetParameters(sum_of_weights * (max_dydz - min_dydz)/bins, -1000.*value(kAlignPhiX), 1000.*value(kResSlopeYSigma), 1000.*value(kResSlopeYGamma));
00528 }
00529 else if (residualsModel() == kGaussPowerTails) {
00530 fit_x = new TF1(name_x.c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_x, max_x, 3);
00531 fit_x->SetParameters(sum_of_weights * (max_x - min_x)/bins, 10.*value(kAlignX), 10.*value(kResidXSigma));
00532 fit_y = new TF1(name_y.c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_y, max_y, 3);
00533 fit_y->SetParameters(sum_of_weights * (max_y - min_y)/bins, 10.*value(kAlignY), 10.*value(kResidYSigma));
00534 fit_dxdz = new TF1(name_dxdz.c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_dxdz, max_dxdz, 3);
00535 fit_dxdz->SetParameters(sum_of_weights * (max_dxdz - min_dxdz)/bins, 1000.*value(kAlignPhiY), 1000.*value(kResSlopeXSigma));
00536 fit_dydz = new TF1(name_dydz.c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_dydz, max_dydz, 3);
00537 fit_dydz->SetParameters(sum_of_weights * (max_dydz - min_dydz)/bins, -1000.*value(kAlignPhiX), 1000.*value(kResSlopeYSigma));
00538 }
00539 else { assert(false); }
00540
00541 fit_x->SetLineColor(2); fit_x->SetLineWidth(2); fit_x->Write();
00542 fit_y->SetLineColor(2); fit_y->SetLineWidth(2); fit_y->Write();
00543 fit_dxdz->SetLineColor(2); fit_dxdz->SetLineWidth(2); fit_dxdz->Write();
00544 fit_dydz->SetLineColor(2); fit_dydz->SetLineWidth(2); fit_dydz->Write();
00545
00546 TF1 *fit_alphax = new TF1(name_alphax.c_str(), "[0] + x*[1]", min_dxdz, max_dxdz);
00547 TF1 *fit_alphay = new TF1(name_alphay.c_str(), "[0] + x*[1]", min_dydz, max_dydz);
00548 double aX = 10.*value(kAlignX), bX = 10.*value(kAlphaX)/1000.;
00549 double aY = 10.*value(kAlignY), bY = 10.*value(kAlphaY)/1000.;
00550 if (residualsModel() == kPureGaussian2D)
00551 {
00552 double sx = 10.*value(kResidXSigma), sy = 1000.*value(kResSlopeXSigma), r = value(kAlphaX);
00553 aX = mean_residualx;
00554 bX = 0.;
00555 if ( sx != 0. )
00556 {
00557 bX = 1./(sy/sx*r);
00558 aX = - bX * mean_resslopex;
00559 }
00560 sx = 10.*value(kResidYSigma); sy = 1000.*value(kResSlopeYSigma); r = value(kAlphaY);
00561 aY = mean_residualx;
00562 bY = 0.;
00563 if ( sx != 0. )
00564 {
00565 bY = 1./(sy/sx*r);
00566 aY = - bY * mean_resslopey;
00567 }
00568 }
00569 fit_alphax->SetParameters(aX, bX);
00570 fit_alphay->SetParameters(aY, bY);
00571
00572 fit_alphax->SetLineColor(2); fit_alphax->SetLineWidth(2); fit_alphax->Write();
00573 fit_alphay->SetLineColor(2); fit_alphay->SetLineWidth(2); fit_alphay->Write();
00574
00575 TProfile *fit_x_trackx = dir->make<TProfile>(name_x_trackx.c_str(), "", 100, min_trackx, max_trackx);
00576 TProfile *fit_y_trackx = dir->make<TProfile>(name_y_trackx.c_str(), "", 100, min_trackx, max_trackx);
00577 TProfile *fit_dxdz_trackx = dir->make<TProfile>(name_dxdz_trackx.c_str(), "", 100, min_trackx, max_trackx);
00578 TProfile *fit_dydz_trackx = dir->make<TProfile>(name_dydz_trackx.c_str(), "", 100, min_trackx, max_trackx);
00579 TProfile *fit_x_tracky = dir->make<TProfile>(name_x_tracky.c_str(), "", 100, min_tracky, max_tracky);
00580 TProfile *fit_y_tracky = dir->make<TProfile>(name_y_tracky.c_str(), "", 100, min_tracky, max_tracky);
00581 TProfile *fit_dxdz_tracky = dir->make<TProfile>(name_dxdz_tracky.c_str(), "", 100, min_tracky, max_tracky);
00582 TProfile *fit_dydz_tracky = dir->make<TProfile>(name_dydz_tracky.c_str(), "", 100, min_tracky, max_tracky);
00583 TProfile *fit_x_trackdxdz = dir->make<TProfile>(name_x_trackdxdz.c_str(), "", 500, min_trackdxdz, max_trackdxdz);
00584 TProfile *fit_y_trackdxdz = dir->make<TProfile>(name_y_trackdxdz.c_str(), "", 500, min_trackdxdz, max_trackdxdz);
00585 TProfile *fit_dxdz_trackdxdz = dir->make<TProfile>(name_dxdz_trackdxdz.c_str(), "", 500, min_trackdxdz, max_trackdxdz);
00586 TProfile *fit_dydz_trackdxdz = dir->make<TProfile>(name_dydz_trackdxdz.c_str(), "", 500, min_trackdxdz, max_trackdxdz);
00587 TProfile *fit_x_trackdydz = dir->make<TProfile>(name_x_trackdydz.c_str(), "", 500, min_trackdydz, max_trackdydz);
00588 TProfile *fit_y_trackdydz = dir->make<TProfile>(name_y_trackdydz.c_str(), "", 500, min_trackdydz, max_trackdydz);
00589 TProfile *fit_dxdz_trackdydz = dir->make<TProfile>(name_dxdz_trackdydz.c_str(), "", 500, min_trackdydz, max_trackdydz);
00590 TProfile *fit_dydz_trackdydz = dir->make<TProfile>(name_dydz_trackdydz.c_str(), "", 500, min_trackdydz, max_trackdydz);
00591
00592 fit_x_trackx->SetLineColor(2); fit_x_trackx->SetLineWidth(2);
00593 fit_y_trackx->SetLineColor(2); fit_y_trackx->SetLineWidth(2);
00594 fit_dxdz_trackx->SetLineColor(2); fit_dxdz_trackx->SetLineWidth(2);
00595 fit_dydz_trackx->SetLineColor(2); fit_dydz_trackx->SetLineWidth(2);
00596 fit_x_tracky->SetLineColor(2); fit_x_tracky->SetLineWidth(2);
00597 fit_y_tracky->SetLineColor(2); fit_y_tracky->SetLineWidth(2);
00598 fit_dxdz_tracky->SetLineColor(2); fit_dxdz_tracky->SetLineWidth(2);
00599 fit_dydz_tracky->SetLineColor(2); fit_dydz_tracky->SetLineWidth(2);
00600 fit_x_trackdxdz->SetLineColor(2); fit_x_trackdxdz->SetLineWidth(2);
00601 fit_y_trackdxdz->SetLineColor(2); fit_y_trackdxdz->SetLineWidth(2);
00602 fit_dxdz_trackdxdz->SetLineColor(2); fit_dxdz_trackdxdz->SetLineWidth(2);
00603 fit_dydz_trackdxdz->SetLineColor(2); fit_dydz_trackdxdz->SetLineWidth(2);
00604 fit_x_trackdydz->SetLineColor(2); fit_x_trackdydz->SetLineWidth(2);
00605 fit_y_trackdydz->SetLineColor(2); fit_y_trackdydz->SetLineWidth(2);
00606 fit_dxdz_trackdydz->SetLineColor(2); fit_dxdz_trackdydz->SetLineWidth(2);
00607 fit_dydz_trackdydz->SetLineColor(2); fit_dydz_trackdydz->SetLineWidth(2);
00608
00609 name_x_trackx += "line";
00610 name_y_trackx += "line";
00611 name_dxdz_trackx += "line";
00612 name_dydz_trackx += "line";
00613 name_x_tracky += "line";
00614 name_y_tracky += "line";
00615 name_dxdz_tracky += "line";
00616 name_dydz_tracky += "line";
00617 name_x_trackdxdz += "line";
00618 name_y_trackdxdz += "line";
00619 name_dxdz_trackdxdz += "line";
00620 name_dydz_trackdxdz += "line";
00621 name_x_trackdydz += "line";
00622 name_y_trackdydz += "line";
00623 name_dxdz_trackdydz += "line";
00624 name_dydz_trackdydz += "line";
00625
00626 TF1 *fitline_x_trackx = new TF1(name_x_trackx.c_str(), residual_x_trackx_TF1, min_trackx, max_trackx, 14);
00627 TF1 *fitline_y_trackx = new TF1(name_y_trackx.c_str(), residual_y_trackx_TF1, min_trackx, max_trackx, 14);
00628 TF1 *fitline_dxdz_trackx = new TF1(name_dxdz_trackx.c_str(), residual_dxdz_trackx_TF1, min_trackx, max_trackx, 14);
00629 TF1 *fitline_dydz_trackx = new TF1(name_dydz_trackx.c_str(), residual_dydz_trackx_TF1, min_trackx, max_trackx, 14);
00630 TF1 *fitline_x_tracky = new TF1(name_x_tracky.c_str(), residual_x_tracky_TF1, min_tracky, max_tracky, 14);
00631 TF1 *fitline_y_tracky = new TF1(name_y_tracky.c_str(), residual_y_tracky_TF1, min_tracky, max_tracky, 14);
00632 TF1 *fitline_dxdz_tracky = new TF1(name_dxdz_tracky.c_str(), residual_dxdz_tracky_TF1, min_tracky, max_tracky, 14);
00633 TF1 *fitline_dydz_tracky = new TF1(name_dydz_tracky.c_str(), residual_dydz_tracky_TF1, min_tracky, max_tracky, 14);
00634 TF1 *fitline_x_trackdxdz = new TF1(name_x_trackdxdz.c_str(), residual_x_trackdxdz_TF1, min_trackdxdz, max_trackdxdz, 14);
00635 TF1 *fitline_y_trackdxdz = new TF1(name_y_trackdxdz.c_str(), residual_y_trackdxdz_TF1, min_trackdxdz, max_trackdxdz, 14);
00636 TF1 *fitline_dxdz_trackdxdz = new TF1(name_dxdz_trackdxdz.c_str(), residual_dxdz_trackdxdz_TF1, min_trackdxdz, max_trackdxdz, 14);
00637 TF1 *fitline_dydz_trackdxdz = new TF1(name_dydz_trackdxdz.c_str(), residual_dydz_trackdxdz_TF1, min_trackdxdz, max_trackdxdz, 14);
00638 TF1 *fitline_x_trackdydz = new TF1(name_x_trackdydz.c_str(), residual_x_trackdydz_TF1, min_trackdydz, max_trackdydz, 14);
00639 TF1 *fitline_y_trackdydz = new TF1(name_y_trackdydz.c_str(), residual_y_trackdydz_TF1, min_trackdydz, max_trackdydz, 14);
00640 TF1 *fitline_dxdz_trackdydz = new TF1(name_dxdz_trackdydz.c_str(), residual_dxdz_trackdydz_TF1, min_trackdydz, max_trackdydz, 14);
00641 TF1 *fitline_dydz_trackdydz = new TF1(name_dydz_trackdydz.c_str(), residual_dydz_trackdydz_TF1, min_trackdydz, max_trackdydz, 14);
00642
00643 std::vector<TF1*> fitlines;
00644 fitlines.push_back(fitline_x_trackx);
00645 fitlines.push_back(fitline_y_trackx);
00646 fitlines.push_back(fitline_dxdz_trackx);
00647 fitlines.push_back(fitline_dydz_trackx);
00648 fitlines.push_back(fitline_x_tracky);
00649 fitlines.push_back(fitline_y_tracky);
00650 fitlines.push_back(fitline_dxdz_tracky);
00651 fitlines.push_back(fitline_dydz_tracky);
00652 fitlines.push_back(fitline_x_trackdxdz);
00653 fitlines.push_back(fitline_y_trackdxdz);
00654 fitlines.push_back(fitline_dxdz_trackdxdz);
00655 fitlines.push_back(fitline_dydz_trackdxdz);
00656 fitlines.push_back(fitline_x_trackdydz);
00657 fitlines.push_back(fitline_y_trackdydz);
00658 fitlines.push_back(fitline_dxdz_trackdydz);
00659 fitlines.push_back(fitline_dydz_trackdydz);
00660
00661 double fitparameters[14] = {value(kAlignX), value(kAlignY), value(kAlignZ), value(kAlignPhiX), value(kAlignPhiY), value(kAlignPhiZ),
00662 mean_trackx, mean_tracky, mean_trackdxdz, mean_trackdydz,
00663 value(kAlphaX), mean_resslopex, value(kAlphaY), mean_resslopey};
00664 if (residualsModel() == kPureGaussian2D) fitparameters[10] = fitparameters[12] = 0.;
00665
00666 for(std::vector<TF1*>::const_iterator itr = fitlines.begin(); itr != fitlines.end(); itr++)
00667 {
00668 (*itr)->SetParameters(fitparameters);
00669 (*itr)->SetLineColor(2);
00670 (*itr)->SetLineWidth(2);
00671 (*itr)->Write();
00672 }
00673
00674 for (std::vector<double*>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter)
00675 {
00676 const double residX = (*resiter)[kResidX];
00677 const double residY = (*resiter)[kResidY];
00678 const double resslopeX = (*resiter)[kResSlopeX];
00679 const double resslopeY = (*resiter)[kResSlopeY];
00680 const double positionX = (*resiter)[kPositionX];
00681 const double positionY = (*resiter)[kPositionY];
00682 const double angleX = (*resiter)[kAngleX];
00683 const double angleY = (*resiter)[kAngleY];
00684 const double redchi2 = (*resiter)[kRedChi2];
00685 double weight = 1./redchi2;
00686 if (!m_weightAlignment) weight = 1.;
00687
00688 if (!m_weightAlignment || TMath::Prob(redchi2*12, 12) < 0.99) {
00689
00690 hist_alphax->Fill(1000.*resslopeX, 10.*residX);
00691 hist_alphay->Fill(1000.*resslopeY, 10.*residY);
00692
00693 double coefX = value(kAlphaX), coefY = value(kAlphaY);
00694 if (residualsModel() == kPureGaussian || residualsModel() == kPureGaussian2D) coefX = coefY = 0.;
00695 double geom_residX = residual_x(value(kAlignX), value(kAlignY), value(kAlignZ), value(kAlignPhiX), value(kAlignPhiY), value(kAlignPhiZ), positionX, positionY, angleX, angleY, coefX, resslopeX);
00696 hist_x->Fill(10.*(residX - geom_residX + value(kAlignX)), weight);
00697 hist_x_trackx->Fill(positionX, 10.*residX, weight);
00698 hist_x_tracky->Fill(positionY, 10.*residX, weight);
00699 hist_x_trackdxdz->Fill(angleX, 10.*residX, weight);
00700 hist_x_trackdydz->Fill(angleY, 10.*residX, weight);
00701 fit_x_trackx->Fill(positionX, 10.*geom_residX, weight);
00702 fit_x_tracky->Fill(positionY, 10.*geom_residX, weight);
00703 fit_x_trackdxdz->Fill(angleX, 10.*geom_residX, weight);
00704 fit_x_trackdydz->Fill(angleY, 10.*geom_residX, weight);
00705
00706 double geom_residY = residual_y(value(kAlignX), value(kAlignY), value(kAlignZ), value(kAlignPhiX), value(kAlignPhiY), value(kAlignPhiZ), positionX, positionY, angleX, angleY, coefY, resslopeY);
00707 hist_y->Fill(10.*(residY - geom_residY + value(kAlignY)), weight);
00708 hist_y_trackx->Fill(positionX, 10.*residY, weight);
00709 hist_y_tracky->Fill(positionY, 10.*residY, weight);
00710 hist_y_trackdxdz->Fill(angleX, 10.*residY, weight);
00711 hist_y_trackdydz->Fill(angleY, 10.*residY, weight);
00712 fit_y_trackx->Fill(positionX, 10.*geom_residY, weight);
00713 fit_y_tracky->Fill(positionY, 10.*geom_residY, weight);
00714 fit_y_trackdxdz->Fill(angleX, 10.*geom_residY, weight);
00715 fit_y_trackdydz->Fill(angleY, 10.*geom_residY, weight);
00716
00717 double geom_resslopeX = residual_dxdz(value(kAlignX), value(kAlignY), value(kAlignZ), value(kAlignPhiX), value(kAlignPhiY), value(kAlignPhiZ), positionX, positionY, angleX, angleY);
00718 hist_dxdz->Fill(1000.*(resslopeX - geom_resslopeX + value(kAlignPhiY)), weight);
00719 hist_dxdz_trackx->Fill(positionX, 1000.*resslopeX, weight);
00720 hist_dxdz_tracky->Fill(positionY, 1000.*resslopeX, weight);
00721 hist_dxdz_trackdxdz->Fill(angleX, 1000.*resslopeX, weight);
00722 hist_dxdz_trackdydz->Fill(angleY, 1000.*resslopeX, weight);
00723 fit_dxdz_trackx->Fill(positionX, 1000.*geom_resslopeX, weight);
00724 fit_dxdz_tracky->Fill(positionY, 1000.*geom_resslopeX, weight);
00725 fit_dxdz_trackdxdz->Fill(angleX, 1000.*geom_resslopeX, weight);
00726 fit_dxdz_trackdydz->Fill(angleY, 1000.*geom_resslopeX, weight);
00727
00728 double geom_resslopeY = residual_dydz(value(kAlignX), value(kAlignY), value(kAlignZ), value(kAlignPhiX), value(kAlignPhiY), value(kAlignPhiZ), positionX, positionY, angleX, angleY);
00729 hist_dydz->Fill(1000.*(resslopeY - geom_resslopeY - value(kAlignPhiX)), weight);
00730 hist_dydz_trackx->Fill(positionX, 1000.*resslopeY, weight);
00731 hist_dydz_tracky->Fill(positionY, 1000.*resslopeY, weight);
00732 hist_dydz_trackdxdz->Fill(angleX, 1000.*resslopeY, weight);
00733 hist_dydz_trackdydz->Fill(angleY, 1000.*resslopeY, weight);
00734 fit_dydz_trackx->Fill(positionX, 1000.*geom_resslopeY, weight);
00735 fit_dydz_tracky->Fill(positionY, 1000.*geom_resslopeY, weight);
00736 fit_dydz_trackdxdz->Fill(angleX, 1000.*geom_resslopeY, weight);
00737 fit_dydz_trackdydz->Fill(angleY, 1000.*geom_resslopeY, weight);
00738 }
00739
00740 hist_x_raw->Fill(10.*residX);
00741 hist_y_raw->Fill(10.*residY);
00742 hist_dxdz_raw->Fill(1000.*resslopeX);
00743 hist_dydz_raw->Fill(1000.*resslopeY);
00744 if (fabs(resslopeX) < 0.005) hist_x_cut->Fill(10.*residX);
00745 if (fabs(resslopeY) < 0.030) hist_y_cut->Fill(10.*residY);
00746 }
00747
00748 double chi2 = 0.;
00749 double ndof = 0.;
00750 for (int i = 1; i <= hist_x->GetNbinsX(); i++) {
00751 double xi = hist_x->GetBinCenter(i);
00752 double yi = hist_x->GetBinContent(i);
00753 double yerri = hist_x->GetBinError(i);
00754 double yth = fit_x->Eval(xi);
00755 if (yerri > 0.) {
00756 chi2 += pow((yth - yi)/yerri, 2);
00757 ndof += 1.;
00758 }
00759 }
00760 for (int i = 1; i <= hist_y->GetNbinsX(); i++) {
00761 double xi = hist_y->GetBinCenter(i);
00762 double yi = hist_y->GetBinContent(i);
00763 double yerri = hist_y->GetBinError(i);
00764 double yth = fit_y->Eval(xi);
00765 if (yerri > 0.) {
00766 chi2 += pow((yth - yi)/yerri, 2);
00767 ndof += 1.;
00768 }
00769 }
00770 for (int i = 1; i <= hist_dxdz->GetNbinsX(); i++) {
00771 double xi = hist_dxdz->GetBinCenter(i);
00772 double yi = hist_dxdz->GetBinContent(i);
00773 double yerri = hist_dxdz->GetBinError(i);
00774 double yth = fit_dxdz->Eval(xi);
00775 if (yerri > 0.) {
00776 chi2 += pow((yth - yi)/yerri, 2);
00777 ndof += 1.;
00778 }
00779 }
00780 for (int i = 1; i <= hist_dydz->GetNbinsX(); i++) {
00781 double xi = hist_dydz->GetBinCenter(i);
00782 double yi = hist_dydz->GetBinContent(i);
00783 double yerri = hist_dydz->GetBinError(i);
00784 double yth = fit_dydz->Eval(xi);
00785 if (yerri > 0.) {
00786 chi2 += pow((yth - yi)/yerri, 2);
00787 ndof += 1.;
00788 }
00789 }
00790 ndof -= npar();
00791
00792 return (ndof > 0. ? chi2 / ndof : -1.);
00793 }
00794
00795
00796 TTree * MuonResiduals6DOFFitter::readNtuple(std::string fname, unsigned int wheel, unsigned int station, unsigned int sector, unsigned int preselected)
00797 {
00798 TFile *f = new TFile(fname.c_str());
00799 TTree *t = (TTree*)f->Get("mual_ttree");
00800
00801
00802 TFile *tmpf = new TFile("small_tree_fit.root","recreate");
00803 assert(tmpf!=0);
00804
00805
00806 TTree *tt = t->CopyTree(Form("is_dt && ring_wheel==%d && station==%d && sector==%d && select==%d", wheel, station, sector, (bool)preselected));
00807
00808 MuonAlignmentTreeRow r;
00809 tt->SetBranchAddress("res_x", &r.res_x);
00810 tt->SetBranchAddress("res_slope_x", &r.res_slope_x);
00811 tt->SetBranchAddress("res_y", &r.res_y);
00812 tt->SetBranchAddress("res_slope_y", &r.res_slope_y);
00813 tt->SetBranchAddress("pos_x", &r.pos_x);
00814 tt->SetBranchAddress("pos_y", &r.pos_y);
00815 tt->SetBranchAddress("angle_x", &r.angle_x);
00816 tt->SetBranchAddress("angle_y", &r.angle_y);
00817 tt->SetBranchAddress("pz", &r.pz);
00818 tt->SetBranchAddress("pt", &r.pt);
00819 tt->SetBranchAddress("q", &r.q);
00820
00821 Long64_t nentries = tt->GetEntries();
00822 for (Long64_t i=0;i<nentries; i++)
00823 {
00824 tt->GetEntry(i);
00825 double *rdata = new double[MuonResiduals6DOFFitter::kNData];
00826 rdata[kResidX] = r.res_x;
00827 rdata[kResidY] = r.res_y;
00828 rdata[kResSlopeX] = r.res_slope_x;
00829 rdata[kResSlopeY] = r.res_slope_y;
00830 rdata[kPositionX] = r.pos_x;
00831 rdata[kPositionY] = r.pos_y;
00832 rdata[kAngleX] = r.angle_x;
00833 rdata[kAngleY] = r.angle_y;
00834 rdata[kRedChi2] = 0.1;
00835 rdata[kPz] = r.pz;
00836 rdata[kPt] = r.pt;
00837 rdata[kCharge] = r.q;
00838 fill(rdata);
00839 }
00840 delete f;
00841
00842 return tt;
00843 }
1.7.3