CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
List of all members | Public Types | Public Member Functions | Protected Member Functions
MuonResidualsBfieldAngleFitter Class Reference

#include <MuonResidualsBfieldAngleFitter.h>

Inheritance diagram for MuonResidualsBfieldAngleFitter:
MuonResidualsFitter

Public Types

enum  {
  kAngle = 0, kBfrompt, kBfrompz, kdEdx,
  kSigma, kGamma, kNPar
}
 
enum  { kResidual = 0, kQoverPt, kQoverPz, kNData }
 
- Public Types inherited from MuonResidualsFitter
enum  {
  kPureGaussian, kPowerLawTails, kROOTVoigt, kGaussPowerTails,
  kPureGaussian2D
}
 
enum  {
  k1DOF, k5DOF, k6DOF, k6DOFrphi,
  kPositionFitter, kAngleFitter, kAngleBfieldFitter
}
 
enum  {
  k1111, k1110, k1100, k1010,
  k0010
}
 

Public Member Functions

bool fit (Alignable *ali)
 
 MuonResidualsBfieldAngleFitter (int residualsModel, int minHitsPerRegion, int useResiduals, bool weightAlignment=true)
 
int ndata ()
 
int npar ()
 
double plot (std::string name, TFileDirectory *dir, Alignable *ali)
 
double sumofweights ()
 
int type () const
 
- Public Member Functions inherited from MuonResidualsFitter
void computeHistogramRangeAndBinning (int which, int &nbins, double &a, double &b)
 
virtual void correctBField ()=0
 
virtual void correctBField (int idx_momentum, int idx_q)
 
TMatrixDSym correlationMatrix ()
 
double covarianceElement (int parNum1, int parNum2)
 
TMatrixDSym covarianceMatrix ()
 
void eraseNotSelectedResiduals ()
 
double errorerror (int parNum)
 
void fill (double *residual)
 
void fix (int parNum, bool val=true)
 
bool fixed (int parNum)
 
void histogramChi2GaussianFit (int which, double &fit_mean, double &fit_sigma)
 
double loglikelihood ()
 
 MuonResidualsFitter (int residualsModel, int minHits, int useResiduals, bool weightAlignment=true)
 
int nfixed ()
 
long numResiduals () const
 
long numsegments ()
 
int parNum2parIdx (int parNum)
 
void plotsimple (std::string name, TFileDirectory *dir, int which, double multiplier)
 
void plotweighted (std::string name, TFileDirectory *dir, int which, int whichredchi2, double multiplier)
 
void read (FILE *file, int which=0)
 
std::vector< double * >
::const_iterator 
residuals_begin () const
 
std::vector< double * >
::const_iterator 
residuals_end () const
 
int residualsModel () const
 
std::vector< bool > & selectedResidualsFlags ()
 
void selectPeakResiduals (double nsigma, int nvar, int *vars)
 
void selectPeakResiduals_simple (double nsigma, int nvar, int *vars)
 
void setPrintLevel (int printLevel)
 
void setStrategy (int strategy)
 
int useRes () const
 
double value (int parNum)
 
void write (FILE *file, int which=0)
 
virtual ~MuonResidualsFitter ()
 

Protected Member Functions

void inform (TMinuit *tMinuit)
 
- Protected Member Functions inherited from MuonResidualsFitter
bool dofit (void(*fcn)(int &, double *, double &, double *, int), std::vector< int > &parNum, std::vector< std::string > &parName, std::vector< double > &start, std::vector< double > &step, std::vector< double > &low, std::vector< double > &high)
 
void initialize_table ()
 

Additional Inherited Members

- Protected Attributes inherited from MuonResidualsFitter
double m_center [20]
 
TMatrixDSym m_cov
 
std::vector< double > m_error
 
std::vector< bool > m_fixed
 
double m_loglikelihood
 
int m_minHits
 
std::map< int, int > m_parNum2parIdx
 
int m_printLevel
 
double m_radii [20]
 
std::vector< double * > m_residuals
 
std::vector< bool > m_residuals_ok
 
int m_residualsModel
 
int m_strategy
 
int m_useResiduals
 
std::vector< double > m_value
 
bool m_weightAlignment
 

Detailed Description

Date:
2010/03/12 22:23:26
Revision:
1.5
Author
J. Pivarski - Texas A&M University pivar.nosp@m.ski@.nosp@m.physi.nosp@m.cs.t.nosp@m.amu.e.nosp@m.du

Definition at line 12 of file MuonResidualsBfieldAngleFitter.h.

Member Enumeration Documentation

anonymous enum
anonymous enum

Constructor & Destructor Documentation

MuonResidualsBfieldAngleFitter::MuonResidualsBfieldAngleFitter ( int  residualsModel,
int  minHitsPerRegion,
int  useResiduals,
bool  weightAlignment = true 
)
inline

Definition at line 31 of file MuonResidualsBfieldAngleFitter.h.

tuple weightAlignment
Definition: align_cfg.py:30
MuonResidualsFitter(int residualsModel, int minHits, int useResiduals, bool weightAlignment=true)
list useResiduals
Definition: align_cfg.py:36

Member Function Documentation

bool MuonResidualsBfieldAngleFitter::fit ( Alignable ali)
virtual

Implements MuonResidualsFitter.

Definition at line 41 of file MuonResidualsBfieldAngleFitter.cc.

References MuonResidualsFitter::dofit(), MuonResidualsFitter::initialize_table(), kAngle, kBfrompt, kBfrompz, kdEdx, kGamma, MuonResidualsFitter::kGaussPowerTails, MuonResidualsFitter::kPureGaussian, kResidual, kSigma, MuonResidualsFitter::m_minHits, timingPdfMaker::mean, MuonResidualsBfieldAngleFitter_FCN(), N, funct::pow(), MuonResidualsFitter::residuals_begin(), MuonResidualsFitter::residuals_end(), MuonResidualsFitter::residualsModel(), mathSSE::sqrt(), dqm_diff::start, plotscripts::stdev(), relval_parameters_module::step, AlCaHLTBitMon_QueryRunRegistry::string, makeMuonMisalignmentScenario::sum_x, and makeMuonMisalignmentScenario::sum_xx.

41  {
42  initialize_table(); // if not already initialized
43 
44  double sum_x = 0.;
45  double sum_xx = 0.;
46  int N = 0;
47 
48  for (std::vector<double*>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
49  const double residual = (*resiter)[kResidual];
50 // const double qoverpt = (*resiter)[kQoverPt];
51 
52  if (fabs(residual) < 0.1) { // truncate at 100 mrad
53  sum_x += residual;
54  sum_xx += residual*residual;
55  N++;
56  }
57  }
58 
59  if (N < m_minHits) return false;
60 
61  // truncated mean and stdev to seed the fit
62  double mean = sum_x/double(N);
63  double stdev = sqrt(sum_xx/double(N) - pow(sum_x/double(N), 2));
64 
65  // refine the standard deviation calculation
66  sum_x = 0.;
67  sum_xx = 0.;
68  N = 0;
69  for (std::vector<double*>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
70  const double residual = (*resiter)[kResidual];
71  if (mean - 1.5*stdev < residual && residual < mean + 1.5*stdev) {
72  sum_x += residual;
73  sum_xx += residual*residual;
74  N++;
75  }
76  }
77  mean = sum_x/double(N);
78  stdev = sqrt(sum_xx/double(N) - pow(sum_x/double(N), 2));
79 
80  sum_x = 0.;
81  sum_xx = 0.;
82  N = 0;
83  for (std::vector<double*>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
84  const double residual = (*resiter)[kResidual];
85  if (mean - 1.5*stdev < residual && residual < mean + 1.5*stdev) {
86  sum_x += residual;
87  sum_xx += residual*residual;
88  N++;
89  }
90  }
91  mean = sum_x/double(N);
92  stdev = sqrt(sum_xx/double(N) - pow(sum_x/double(N), 2));
93 
94  std::vector<int> parNum;
95  std::vector<std::string> parName;
96  std::vector<double> start;
97  std::vector<double> step;
98  std::vector<double> low;
99  std::vector<double> high;
100 
101  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.);
102  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.);
103  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.);
104  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.);
105  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.);
107  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.);
108  }
109 
110  return dofit(&MuonResidualsBfieldAngleFitter_FCN, parNum, parName, start, step, low, high);
111 }
tuple start
Check for commandline option errors.
Definition: dqm_diff.py:58
bool dofit(void(*fcn)(int &, double *, double &, double *, int), std::vector< int > &parNum, std::vector< std::string > &parName, std::vector< double > &start, std::vector< double > &step, std::vector< double > &low, std::vector< double > &high)
T sqrt(T t)
Definition: SSEVec.h:48
#define N
Definition: blowfish.cc:9
std::vector< double * >::const_iterator residuals_end() const
std::vector< double * >::const_iterator residuals_begin() const
void MuonResidualsBfieldAngleFitter_FCN(int &npar, double *gin, double &fval, double *par, int iflag)
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
void MuonResidualsBfieldAngleFitter::inform ( TMinuit *  tMinuit)
protectedvirtual

Implements MuonResidualsFitter.

Definition at line 5 of file MuonResidualsBfieldAngleFitter.cc.

References MuonResidualsBfieldAngleFitter_TMinuit.

5  {
7 }
static TMinuit * MuonResidualsBfieldAngleFitter_TMinuit
int MuonResidualsBfieldAngleFitter::ndata ( )
inlinevirtual
int MuonResidualsBfieldAngleFitter::npar ( )
inlinevirtual
double MuonResidualsBfieldAngleFitter::plot ( std::string  name,
TFileDirectory dir,
Alignable ali 
)
virtual

Implements MuonResidualsFitter.

Definition at line 113 of file MuonResidualsBfieldAngleFitter.cc.

References kAngle, kBfrompt, kBfrompz, kdEdx, kGamma, MuonResidualsFitter::kGaussPowerTails, MuonResidualsFitter::kPowerLawTails, MuonResidualsFitter::kPureGaussian, kQoverPt, kQoverPz, kResidual, MuonResidualsFitter::kROOTVoigt, kSigma, TFileDirectory::make(), MuonResidualsFitter_GaussPowerTails_TF1(), MuonResidualsFitter_powerLawTails_TF1(), MuonResidualsFitter_pureGaussian_TF1(), MuonResidualsFitter_ROOTVoigt_TF1(), NULL, MuonResidualsFitter::numResiduals(), MuonResidualsFitter::residuals_begin(), MuonResidualsFitter::residuals_end(), MuonResidualsFitter::residualsModel(), AlCaHLTBitMon_QueryRunRegistry::string, and MuonResidualsFitter::value().

Referenced by cuy.FindIssue::__init__().

113  {
114  std::stringstream raw_name, narrowed_name, qoverpt_name, qoverpz_name, psquared_name;
115  raw_name << name << "_raw";
116  narrowed_name << name << "_narrowed";
117  qoverpt_name << name << "_qoverpt";
118  qoverpz_name << name << "_qoverpz";
119  psquared_name << name << "_psquared";
120 
121  TH1F *raw_hist = dir->make<TH1F>(raw_name.str().c_str(), (raw_name.str() + std::string(" (mrad)")).c_str(), 100, -100., 100.);
122  TH1F *narrowed_hist = dir->make<TH1F>(narrowed_name.str().c_str(), (narrowed_name.str() + std::string(" (mrad)")).c_str(), 100, -100., 100.);
123  TProfile *qoverpt_hist = dir->make<TProfile>(qoverpt_name.str().c_str(), (qoverpt_name.str() + std::string(" (mrad)")).c_str(), 100, -0.05, 0.05);
124  TProfile *qoverpz_hist = dir->make<TProfile>(qoverpz_name.str().c_str(), (qoverpz_name.str() + std::string(" (mrad)")).c_str(), 100, -0.05, 0.05);
125  TProfile *psquared_hist = dir->make<TProfile>(psquared_name.str().c_str(), (psquared_name.str() + std::string(" (mrad)")).c_str(), 100, -0.05, 0.05);
126 
127  narrowed_name << "fit";
128  qoverpt_name << "fit";
129  qoverpz_name << "fit";
130  psquared_name << "fit";
131 
132  double scale_factor = double(numResiduals()) * (100. - -100.)/100; // (max - min)/nbins
133 
134  TF1 *narrowed_fit = NULL;
135  if (residualsModel() == kPureGaussian) {
136  narrowed_fit = new TF1(narrowed_name.str().c_str(), MuonResidualsFitter_pureGaussian_TF1, -100., 100., 3);
137  narrowed_fit->SetParameters(scale_factor, value(kAngle) * 1000., value(kSigma) * 1000.);
138  narrowed_fit->Write();
139  }
140  else if (residualsModel() == kPowerLawTails) {
141  narrowed_fit = new TF1(narrowed_name.str().c_str(), MuonResidualsFitter_powerLawTails_TF1, -100., 100., 4);
142  narrowed_fit->SetParameters(scale_factor, value(kAngle) * 1000., value(kSigma) * 1000., value(kGamma) * 1000.);
143  narrowed_fit->Write();
144  }
145  else if (residualsModel() == kROOTVoigt) {
146  narrowed_fit = new TF1(narrowed_name.str().c_str(), MuonResidualsFitter_ROOTVoigt_TF1, -100., 100., 4);
147  narrowed_fit->SetParameters(scale_factor, value(kAngle) * 1000., value(kSigma) * 1000., value(kGamma) * 1000.);
148  narrowed_fit->Write();
149  }
150  else if (residualsModel() == kGaussPowerTails) {
151  narrowed_fit = new TF1(narrowed_name.str().c_str(), MuonResidualsFitter_GaussPowerTails_TF1, -100., 100., 3);
152  narrowed_fit->SetParameters(scale_factor, value(kAngle) * 1000., value(kSigma) * 1000.);
153  narrowed_fit->Write();
154  }
155 
156  TF1 *qoverpt_fit = new TF1(qoverpt_name.str().c_str(), "[0]+x*[1]", -0.05, 0.05);
157  qoverpt_fit->SetParameters(value(kAngle) * 1000., value(kBfrompt) * 1000.);
158  qoverpt_fit->Write();
159 
160  TF1 *qoverpz_fit = new TF1(qoverpz_name.str().c_str(), "[0]+x*[1]", -0.05, 0.05);
161  qoverpz_fit->SetParameters(value(kAngle) * 1000., value(kBfrompz) * 1000.);
162  qoverpz_fit->Write();
163 
164  TF1 *psquared_fit = new TF1(psquared_name.str().c_str(), "[0]+[1]*x**2", -0.05, 0.05);
165  psquared_fit->SetParameters(value(kAngle) * 1000., value(kdEdx) * 1000.);
166  psquared_fit->Write();
167 
168  for (std::vector<double*>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
169  const double raw_residual = (*resiter)[kResidual];
170  const double qoverpt = (*resiter)[kQoverPt];
171  const double qoverpz = (*resiter)[kQoverPz];
172  const double psquared = (1./qoverpt/qoverpt + 1./qoverpz/qoverpz) * (qoverpt > 0. ? 1. : -1.);
173 
174  double qoverpt_correction = value(kBfrompt) * qoverpt;
175  double qoverpz_correction = value(kBfrompz) * qoverpz;
176  double dEdx_correction = value(kdEdx) * psquared;
177  double corrected_residual = raw_residual - qoverpt_correction - qoverpz_correction - dEdx_correction;
178 
179  raw_hist->Fill(raw_residual * 1000.);
180  narrowed_hist->Fill(corrected_residual * 1000.);
181 
182  qoverpt_hist->Fill(qoverpt, (raw_residual - qoverpz_correction - dEdx_correction) * 1000.);
183  qoverpz_hist->Fill(qoverpz, (raw_residual - qoverpt_correction - dEdx_correction) * 1000.);
184  psquared_hist->Fill(psquared, (raw_residual - qoverpt_correction - qoverpz_correction) * 1000.);
185  }
186 
187  return 0.;
188 }
Double_t MuonResidualsFitter_powerLawTails_TF1(Double_t *xvec, Double_t *par)
#define NULL
Definition: scimark2.h:8
double value(int parNum)
Double_t MuonResidualsFitter_GaussPowerTails_TF1(Double_t *xvec, Double_t *par)
T * make(const Args &...args) const
make new ROOT object
Double_t MuonResidualsFitter_pureGaussian_TF1(Double_t *xvec, Double_t *par)
std::vector< double * >::const_iterator residuals_end() const
Double_t MuonResidualsFitter_ROOTVoigt_TF1(Double_t *xvec, Double_t *par)
std::vector< double * >::const_iterator residuals_begin() const
double MuonResidualsBfieldAngleFitter::sumofweights ( )
inlinevirtual

Implements MuonResidualsFitter.

Definition at line 44 of file MuonResidualsBfieldAngleFitter.h.

References MuonResidualsFitter::numResiduals().

44 { return numResiduals(); }
int MuonResidualsBfieldAngleFitter::type ( ) const
inlinevirtual