Reference Manual

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