CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_2_7_hltpatch1/src/Alignment/MuonAlignmentAlgorithms/src/MuonResidualsBfieldAngleFitter.cc

Go to the documentation of this file.
00001 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResidualsBfieldAngleFitter.h"
00002 
00003 static TMinuit *MuonResidualsBfieldAngleFitter_TMinuit;
00004 
00005 void MuonResidualsBfieldAngleFitter::inform(TMinuit *tMinuit) {
00006   MuonResidualsBfieldAngleFitter_TMinuit = tMinuit;
00007 }
00008 
00009 void MuonResidualsBfieldAngleFitter_FCN(int &npar, double *gin, double &fval, double *par, int iflag) {
00010   MuonResidualsFitterFitInfo *fitinfo = (MuonResidualsFitterFitInfo*)(MuonResidualsBfieldAngleFitter_TMinuit->GetObjectFit());
00011   MuonResidualsFitter *fitter = fitinfo->fitter();
00012 
00013   fval = 0.;
00014   for (std::vector<double*>::const_iterator resiter = fitter->residuals_begin();  resiter != fitter->residuals_end();  ++resiter) {
00015     const double residual = (*resiter)[MuonResidualsBfieldAngleFitter::kResidual];
00016     const double qoverpt = (*resiter)[MuonResidualsBfieldAngleFitter::kQoverPt];
00017     const double qoverpz = (*resiter)[MuonResidualsBfieldAngleFitter::kQoverPz];
00018     
00019     double center = 0.;
00020     center += par[MuonResidualsBfieldAngleFitter::kAngle];
00021     center += par[MuonResidualsBfieldAngleFitter::kBfrompt] * qoverpt;
00022     center += par[MuonResidualsBfieldAngleFitter::kBfrompz] * qoverpz;
00023     center += par[MuonResidualsBfieldAngleFitter::kdEdx] * (1./qoverpt/qoverpt + 1./qoverpz/qoverpz) * (qoverpt > 0. ? 1. : -1.);
00024     
00025     if (fitter->residualsModel() == MuonResidualsFitter::kPureGaussian) {
00026       fval += -MuonResidualsFitter_logPureGaussian(residual, center, par[MuonResidualsBfieldAngleFitter::kSigma]);
00027     }
00028     else if (fitter->residualsModel() == MuonResidualsFitter::kPowerLawTails) {
00029       fval += -MuonResidualsFitter_logPowerLawTails(residual, center, par[MuonResidualsBfieldAngleFitter::kSigma], par[MuonResidualsBfieldAngleFitter::kGamma]);
00030     }
00031     else if (fitter->residualsModel() == MuonResidualsFitter::kROOTVoigt) {
00032       fval += -MuonResidualsFitter_logROOTVoigt(residual, center, par[MuonResidualsBfieldAngleFitter::kSigma], par[MuonResidualsBfieldAngleFitter::kGamma]);
00033     }
00034     else if (fitter->residualsModel() == MuonResidualsFitter::kGaussPowerTails) {
00035       fval += -MuonResidualsFitter_logGaussPowerTails(residual, center, par[MuonResidualsBfieldAngleFitter::kSigma]);
00036     }
00037     else { assert(false); }
00038   }
00039 }
00040 
00041 bool MuonResidualsBfieldAngleFitter::fit(Alignable *ali) {
00042   initialize_table();  // if not already initialized
00043 
00044   double sum_x = 0.;
00045   double sum_xx = 0.;
00046   int N = 0;
00047 
00048   for (std::vector<double*>::const_iterator resiter = residuals_begin();  resiter != residuals_end();  ++resiter) {
00049     const double residual = (*resiter)[kResidual];
00050 //     const double qoverpt = (*resiter)[kQoverPt];
00051 
00052     if (fabs(residual) < 0.1) {  // truncate at 100 mrad
00053       sum_x += residual;
00054       sum_xx += residual*residual;
00055       N++;
00056     }
00057   }
00058 
00059   if (N < m_minHits) return false;
00060 
00061   // truncated mean and stdev to seed the fit
00062   double mean = sum_x/double(N);
00063   double stdev = sqrt(sum_xx/double(N) - pow(sum_x/double(N), 2));
00064 
00065   // refine the standard deviation calculation
00066   sum_x = 0.;
00067   sum_xx = 0.;
00068   N = 0;
00069   for (std::vector<double*>::const_iterator resiter = residuals_begin();  resiter != residuals_end();  ++resiter) {
00070     const double residual = (*resiter)[kResidual];
00071     if (mean - 1.5*stdev < residual  &&  residual < mean + 1.5*stdev) {
00072       sum_x += residual;
00073       sum_xx += residual*residual;
00074       N++;
00075     }
00076   }
00077   mean = sum_x/double(N);
00078   stdev = sqrt(sum_xx/double(N) - pow(sum_x/double(N), 2));
00079 
00080   sum_x = 0.;
00081   sum_xx = 0.;
00082   N = 0;
00083   for (std::vector<double*>::const_iterator resiter = residuals_begin();  resiter != residuals_end();  ++resiter) {
00084     const double residual = (*resiter)[kResidual];
00085     if (mean - 1.5*stdev < residual  &&  residual < mean + 1.5*stdev) {
00086       sum_x += residual;
00087       sum_xx += residual*residual;
00088       N++;
00089     }
00090   }
00091   mean = sum_x/double(N);
00092   stdev = sqrt(sum_xx/double(N) - pow(sum_x/double(N), 2));
00093 
00094   std::vector<int> parNum;
00095   std::vector<std::string> parName;
00096   std::vector<double> start;
00097   std::vector<double> step;
00098   std::vector<double> low;
00099   std::vector<double> high;
00100 
00101   parNum.push_back(kAngle);     parName.push_back(std::string("angle"));     start.push_back(mean);   step.push_back(0.1);               low.push_back(0.);    high.push_back(0.);
00102   parNum.push_back(kBfrompt);   parName.push_back(std::string("bfrompt"));   start.push_back(0.);       step.push_back(0.1*stdev/0.05);  low.push_back(0.);    high.push_back(0.);
00103   parNum.push_back(kBfrompz);   parName.push_back(std::string("bfrompz"));   start.push_back(0.);       step.push_back(0.1*stdev/0.05);  low.push_back(0.);    high.push_back(0.);
00104   parNum.push_back(kdEdx);      parName.push_back(std::string("dEdx"));      start.push_back(0.);       step.push_back(0.1*stdev/0.05);  low.push_back(0.);    high.push_back(0.);
00105   parNum.push_back(kSigma);     parName.push_back(std::string("sigma"));     start.push_back(stdev);  step.push_back(0.1*stdev);       low.push_back(0.);    high.push_back(0.);
00106   if (residualsModel() != kPureGaussian && residualsModel() != kGaussPowerTails) {
00107   parNum.push_back(kGamma);     parName.push_back(std::string("gamma"));     start.push_back(stdev);  step.push_back(0.1*stdev);       low.push_back(0.);    high.push_back(0.);
00108   }
00109 
00110   return dofit(&MuonResidualsBfieldAngleFitter_FCN, parNum, parName, start, step, low, high);
00111 }
00112 
00113 double MuonResidualsBfieldAngleFitter::plot(std::string name, TFileDirectory *dir, Alignable *ali) {
00114   std::stringstream raw_name, narrowed_name, qoverpt_name, qoverpz_name, psquared_name;
00115   raw_name << name << "_raw";
00116   narrowed_name << name << "_narrowed";
00117   qoverpt_name << name << "_qoverpt";
00118   qoverpz_name << name << "_qoverpz";
00119   psquared_name << name << "_psquared";
00120 
00121   TH1F *raw_hist = dir->make<TH1F>(raw_name.str().c_str(), (raw_name.str() + std::string(" (mrad)")).c_str(), 100, -100., 100.);
00122   TH1F *narrowed_hist = dir->make<TH1F>(narrowed_name.str().c_str(), (narrowed_name.str() + std::string(" (mrad)")).c_str(), 100, -100., 100.);
00123   TProfile *qoverpt_hist = dir->make<TProfile>(qoverpt_name.str().c_str(), (qoverpt_name.str() + std::string(" (mrad)")).c_str(), 100, -0.05, 0.05);
00124   TProfile *qoverpz_hist = dir->make<TProfile>(qoverpz_name.str().c_str(), (qoverpz_name.str() + std::string(" (mrad)")).c_str(), 100, -0.05, 0.05);
00125   TProfile *psquared_hist = dir->make<TProfile>(psquared_name.str().c_str(), (psquared_name.str() + std::string(" (mrad)")).c_str(), 100, -0.05, 0.05);
00126 
00127   narrowed_name << "fit";
00128   qoverpt_name << "fit";
00129   qoverpz_name << "fit";
00130   psquared_name << "fit";
00131 
00132   double scale_factor = double(numResiduals()) * (100. - -100.)/100;   // (max - min)/nbins
00133 
00134   TF1 *narrowed_fit = NULL;
00135   if (residualsModel() == kPureGaussian) {
00136     narrowed_fit = new TF1(narrowed_name.str().c_str(), MuonResidualsFitter_pureGaussian_TF1, -100., 100., 3);
00137     narrowed_fit->SetParameters(scale_factor, value(kAngle) * 1000., value(kSigma) * 1000.);
00138     narrowed_fit->Write();
00139   }
00140   else if (residualsModel() == kPowerLawTails) {
00141     narrowed_fit = new TF1(narrowed_name.str().c_str(), MuonResidualsFitter_powerLawTails_TF1, -100., 100., 4);
00142     narrowed_fit->SetParameters(scale_factor, value(kAngle) * 1000., value(kSigma) * 1000., value(kGamma) * 1000.);
00143     narrowed_fit->Write();
00144   }
00145   else if (residualsModel() == kROOTVoigt) {
00146     narrowed_fit = new TF1(narrowed_name.str().c_str(), MuonResidualsFitter_ROOTVoigt_TF1, -100., 100., 4);
00147     narrowed_fit->SetParameters(scale_factor, value(kAngle) * 1000., value(kSigma) * 1000., value(kGamma) * 1000.);
00148     narrowed_fit->Write();
00149   }
00150   else if (residualsModel() == kGaussPowerTails) {
00151     narrowed_fit = new TF1(narrowed_name.str().c_str(), MuonResidualsFitter_GaussPowerTails_TF1, -100., 100., 3);
00152     narrowed_fit->SetParameters(scale_factor, value(kAngle) * 1000., value(kSigma) * 1000.);
00153     narrowed_fit->Write();
00154   }
00155 
00156   TF1 *qoverpt_fit = new TF1(qoverpt_name.str().c_str(), "[0]+x*[1]", -0.05, 0.05);
00157   qoverpt_fit->SetParameters(value(kAngle) * 1000., value(kBfrompt) * 1000.);
00158   qoverpt_fit->Write();
00159 
00160   TF1 *qoverpz_fit = new TF1(qoverpz_name.str().c_str(), "[0]+x*[1]", -0.05, 0.05);
00161   qoverpz_fit->SetParameters(value(kAngle) * 1000., value(kBfrompz) * 1000.);
00162   qoverpz_fit->Write();
00163 
00164   TF1 *psquared_fit = new TF1(psquared_name.str().c_str(), "[0]+[1]*x**2", -0.05, 0.05);
00165   psquared_fit->SetParameters(value(kAngle) * 1000., value(kdEdx) * 1000.);
00166   psquared_fit->Write();
00167 
00168   for (std::vector<double*>::const_iterator resiter = residuals_begin();  resiter != residuals_end();  ++resiter) {
00169     const double raw_residual = (*resiter)[kResidual];
00170     const double qoverpt = (*resiter)[kQoverPt];
00171     const double qoverpz = (*resiter)[kQoverPz];
00172     const double psquared = (1./qoverpt/qoverpt + 1./qoverpz/qoverpz) * (qoverpt > 0. ? 1. : -1.);
00173 
00174     double qoverpt_correction = value(kBfrompt) * qoverpt;
00175     double qoverpz_correction = value(kBfrompz) * qoverpz;
00176     double dEdx_correction = value(kdEdx) * psquared;
00177     double corrected_residual = raw_residual - qoverpt_correction - qoverpz_correction - dEdx_correction;
00178 
00179     raw_hist->Fill(raw_residual * 1000.);
00180     narrowed_hist->Fill(corrected_residual * 1000.);
00181 
00182     qoverpt_hist->Fill(qoverpt, (raw_residual - qoverpz_correction - dEdx_correction) * 1000.);
00183     qoverpz_hist->Fill(qoverpz, (raw_residual - qoverpt_correction - dEdx_correction) * 1000.);
00184     psquared_hist->Fill(psquared, (raw_residual - qoverpt_correction - qoverpz_correction) * 1000.);
00185   }
00186 
00187   return 0.;
00188 }