CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
List of all members | Public Types | Public Member Functions | Protected Member Functions
MuonResiduals6DOFFitter Class Reference

#include <MuonResiduals6DOFFitter.h>

Inheritance diagram for MuonResiduals6DOFFitter:
MuonResidualsFitter

Public Types

enum  {
  kAlignX = 0, kAlignY, kAlignZ, kAlignPhiX,
  kAlignPhiY, kAlignPhiZ, kResidXSigma, kResidYSigma,
  kResSlopeXSigma, kResSlopeYSigma, kAlphaX, kAlphaY,
  kResidXGamma, kResidYGamma, kResSlopeXGamma, kResSlopeYGamma,
  kNPar
}
 
enum  {
  kResidX = 0, kResidY, kResSlopeX, kResSlopeY,
  kPositionX, kPositionY, kAngleX, kAngleY,
  kRedChi2, kPz, kPt, kCharge,
  kStation, kWheel, kSector, kChambW,
  kChambl, 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, k1000, k0100
}
 

Public Member Functions

void correctBField () override
 
bool fit (Alignable *ali) override
 
 MuonResiduals6DOFFitter (int residualsModel, int minHits, int useResiduals, bool weightAlignment=true)
 
int ndata () override
 
int npar () override
 
double plot (std::string name, TFileDirectory *dir, Alignable *ali) override
 
TTree * readNtuple (std::string fname, unsigned int wheel, unsigned int station, unsigned int sector, unsigned int preselected=1)
 
double sumofweights () override
 
int type () const override
 
 ~MuonResiduals6DOFFitter () override
 
- Public Member Functions inherited from MuonResidualsFitter
void computeHistogramRangeAndBinning (int which, int &nbins, double &a, double &b)
 
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 fiducialCuts (double xMin=-80.0, double xMax=80.0, double yMin=-80.0, double yMax=80.0, bool fidcut1=false)
 
void fill (double *residual)
 
void fix (int parNum, bool dofix=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 setInitialValue (int parNum, double value)
 
void setPrintLevel (int printLevel)
 
void setStrategy (int strategy)
 
int useRes (int pattern=-1)
 
double value (int parNum)
 
void write (FILE *file, int which=0)
 
virtual ~MuonResidualsFitter ()
 

Protected Member Functions

void inform (TMinuit *tMinuit) override
 
- 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, double > m_parNum2InitValue
 
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: Thu Apr 16 14:20:58 CDT 2009

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 18 of file MuonResiduals6DOFFitter.h.

Member Enumeration Documentation

anonymous enum
Enumerator
kAlignX 
kAlignY 
kAlignZ 
kAlignPhiX 
kAlignPhiY 
kAlignPhiZ 
kResidXSigma 
kResidYSigma 
kResSlopeXSigma 
kResSlopeYSigma 
kAlphaX 
kAlphaY 
kResidXGamma 
kResidYGamma 
kResSlopeXGamma 
kResSlopeYGamma 
kNPar 

Definition at line 20 of file MuonResiduals6DOFFitter.h.

20  {
21  kAlignX = 0,
22  kAlignY,
23  kAlignZ,
24  kAlignPhiX,
25  kAlignPhiY,
26  kAlignPhiZ,
27  kResidXSigma,
28  kResidYSigma,
29  kResSlopeXSigma,
30  kResSlopeYSigma,
31  kAlphaX,
32  kAlphaY,
33  kResidXGamma,
34  kResidYGamma,
35  kResSlopeXGamma,
36  kResSlopeYGamma,
37  kNPar
38  };
anonymous enum
Enumerator
kResidX 
kResidY 
kResSlopeX 
kResSlopeY 
kPositionX 
kPositionY 
kAngleX 
kAngleY 
kRedChi2 
kPz 
kPt 
kCharge 
kStation 
kWheel 
kSector 
kChambW 
kChambl 
kNData 

Definition at line 40 of file MuonResiduals6DOFFitter.h.

40  {
41  kResidX = 0,
42  kResidY,
43  kResSlopeX,
44  kResSlopeY,
45  kPositionX,
46  kPositionY,
47  kAngleX,
48  kAngleY,
49  kRedChi2,
50  kPz,
51  kPt,
52  kCharge,
53  kStation,
54  kWheel,
55  kSector,
56  kChambW,
57  kChambl,
58  kNData
59  };

Constructor & Destructor Documentation

MuonResiduals6DOFFitter::MuonResiduals6DOFFitter ( int  residualsModel,
int  minHits,
int  useResiduals,
bool  weightAlignment = true 
)
inline

Definition at line 61 of file MuonResiduals6DOFFitter.h.

62  : MuonResidualsFitter(residualsModel, minHits, useResiduals, weightAlignment) {}
align_cfg.weightAlignment
tuple weightAlignment
Definition: align_cfg.py:30
OccupancyClient_cfi.minHits
int minHits
Definition: OccupancyClient_cfi.py:5
MuonResidualsFitter::MuonResidualsFitter
MuonResidualsFitter(int residualsModel, int minHits, int useResiduals, bool weightAlignment=true)
Definition: MuonResidualsFitter.cc:151
align_cfg.useResiduals
list useResiduals
Definition: align_cfg.py:36
MuonResiduals6DOFFitter::~MuonResiduals6DOFFitter ( )
inlineoverride

Definition at line 63 of file MuonResiduals6DOFFitter.h.

63 {}

Member Function Documentation

void MuonResiduals6DOFFitter::correctBField ( )
overridevirtual

Implements MuonResidualsFitter.

Definition at line 264 of file MuonResiduals6DOFFitter.cc.

References MuonResidualsFitter::correctBField(), kCharge, and kPt.

264 { MuonResidualsFitter::correctBField(kPt, kCharge); }
MuonResidualsFitter::correctBField
virtual void correctBField()=0
bool MuonResiduals6DOFFitter::fit ( Alignable *  ali)
overridevirtual

Implements MuonResidualsFitter.

Definition at line 285 of file MuonResiduals6DOFFitter.cc.

References MuonResidualsFitter::dofit(), MuonResidualsFitter::fix(), mps_fire::i, MuonResidualsFitter::initialize_table(), MuonResidualsFitter::k0010, MuonResidualsFitter::k1010, MuonResidualsFitter::k1100, MuonResidualsFitter::k1110, MuonResidualsFitter::k1111, kAlignPhiX, kAlignPhiY, kAlignPhiZ, kAlignX, kAlignY, kAlignZ, kAlphaX, kAlphaY, MuonResidualsFitter::kPowerLawTails, MuonResidualsFitter::kPureGaussian2D, kResidXGamma, kResidXSigma, kResidYGamma, kResidYSigma, kResSlopeXGamma, kResSlopeXSigma, kResSlopeYGamma, kResSlopeYSigma, MuonResidualsFitter::kROOTVoigt, MuonResiduals6DOFFitter_FCN(), mergeVDriftHistosByStation::name, names, pileupDistInMC::num, MuonResidualsFitter::residualsModel(), command_line::start, relval_machine::steps, AlCaHLTBitMon_QueryRunRegistry::string, sumofweights(), and MuonResidualsFitter::useRes().

Referenced by trackingPlots.Iteration::modules().

285  {
286  initialize_table(); // if not already initialized
287  sumofweights();
288 
289  double resx_std = 0.5;
290  double resy_std = 3.0;
291  double resslopex_std = 0.002;
292  double resslopey_std = 0.005;
293 
294  int nums[16] = {kAlignX,
295  kAlignY,
296  kAlignZ,
297  kAlignPhiX,
298  kAlignPhiY,
299  kAlignPhiZ,
300  kResidXSigma,
301  kResidYSigma,
302  kResSlopeXSigma,
303  kResSlopeYSigma,
304  kAlphaX,
305  kAlphaY,
306  kResidXGamma,
307  kResidYGamma,
308  kResSlopeXGamma,
309  kResSlopeYGamma};
310  std::string names[16] = {"AlignX",
311  "AlignY",
312  "AlignZ",
313  "AlignPhiX",
314  "AlignPhiY",
315  "AlignPhiZ",
316  "ResidXSigma",
317  "ResidYSigma",
318  "ResSlopeXSigma",
319  "ResSlopeYSigma",
320  "AlphaX",
321  "AlphaY",
322  "ResidXGamma",
323  "ResidYGamma",
324  "ResSlopeXGamma",
325  "ResSlopeYGamma"};
326  double starts[16] = {0.,
327  0.,
328  0.,
329  0.,
330  0.,
331  0.,
332  resx_std,
333  resy_std,
334  resslopex_std,
335  resslopey_std,
336  0.,
337  0.,
338  0.1 * resx_std,
339  0.1 * resy_std,
340  0.1 * resslopex_std,
341  0.1 * resslopey_std};
342  double steps[16] = {0.1,
343  0.1,
344  0.1,
345  0.001,
346  0.001,
347  0.001,
348  0.001 * resx_std,
349  0.001 * resy_std,
350  0.001 * resslopex_std,
351  0.001 * resslopey_std,
352  0.001,
353  0.001,
354  0.01 * resx_std,
355  0.01 * resy_std,
356  0.01 * resslopex_std,
357  0.01 * resslopey_std};
358  double lows[16] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -1., -1., 0., 0., 0., 0.};
359  double highs[16] = {0., 0., 0., 0., 0., 0., 10., 10., 0.1, 0.1, 1., 1., 0., 0., 0., 0.};
360 
361  std::vector<int> num(nums, nums + 6);
362  std::vector<std::string> name(names, names + 6);
363  std::vector<double> start(starts, starts + 6);
364  std::vector<double> step(steps, steps + 6);
365  std::vector<double> low(lows, lows + 6);
366  std::vector<double> high(highs, highs + 6);
367 
368  bool add_alpha = (residualsModel() == kPureGaussian2D);
369  bool add_gamma = (residualsModel() == kROOTVoigt || residualsModel() == kPowerLawTails);
370 
371  int idx[8], ni = 0;
372  if (useRes() == k1111) {
373  for (ni = 0; ni < 4; ni++)
374  idx[ni] = ni + 6;
375  if (add_alpha)
376  for (; ni < 6; ni++)
377  idx[ni] = ni + 6;
378  else if (add_gamma)
379  for (; ni < 8; ni++)
380  idx[ni] = ni + 8;
381  if (!add_alpha)
382  fix(kAlphaX);
383  if (!add_alpha)
384  fix(kAlphaY);
385  } else if (useRes() == k1110) {
386  for (ni = 0; ni < 3; ni++)
387  idx[ni] = ni + 6;
388  if (add_alpha)
389  idx[ni++] = 10;
390  else if (add_gamma)
391  for (; ni < 6; ni++)
392  idx[ni] = ni + 9;
393  fix(kResSlopeYSigma);
394  fix(kAlphaY);
395  if (!add_alpha)
396  fix(kAlphaX);
397  } else if (useRes() == k1100) {
398  for (ni = 0; ni < 2; ni++)
399  idx[ni] = ni + 6;
400  if (add_gamma)
401  for (; ni < 4; ni++)
402  idx[ni] = ni + 10;
403  fix(kResSlopeXSigma);
404  fix(kResSlopeYSigma);
405  fix(kAlphaX);
406  fix(kAlphaY);
407  } else if (useRes() == k1010) {
408  idx[ni++] = 6;
409  idx[ni++] = 8;
410  if (add_alpha)
411  idx[ni++] = 10;
412  if (add_gamma) {
413  idx[ni++] = 12;
414  idx[ni++] = 14;
415  }
416  fix(kResidYSigma);
417  fix(kResSlopeYSigma);
418  fix(kAlphaY);
419  if (!add_alpha)
420  fix(kAlphaX);
421  } else if (useRes() == k0010) {
422  idx[ni++] = 8;
423  if (add_gamma)
424  idx[ni++] = 14;
425  fix(kResidXSigma);
426  fix(kResidYSigma);
427  fix(kResSlopeYSigma);
428  fix(kAlphaX);
429  fix(kAlphaY);
430  }
431  for (int i = 0; i < ni; i++) {
432  num.push_back(nums[idx[i]]);
433  name.push_back(names[idx[i]]);
434  start.push_back(starts[idx[i]]);
435  step.push_back(steps[idx[i]]);
436  low.push_back(lows[idx[i]]);
437  high.push_back(highs[idx[i]]);
438  }
439 
440  return dofit(&MuonResiduals6DOFFitter_FCN, num, name, start, step, low, high);
441 }
mps_fire.i
i
Definition: mps_fire.py:428
relval_machine.steps
tuple steps
Definition: relval_machine.py:9
names
const std::string names[nVars_]
Definition: PhotonIDValueMapProducer.cc:124
MuonResidualsFitter::dofit
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)
Definition: MuonResidualsFitter.cc:272
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
command_line.start
start
Definition: command_line.py:165
MuonResidualsFitter::fix
void fix(int parNum, bool dofix=true)
Definition: MuonResidualsFitter.cc:171
MuonResidualsFitter::useRes
int useRes(int pattern=-1)
Definition: MuonResidualsFitter.h:119
step
step
Definition: StallMonitor.cc:98
MuonResiduals6DOFFitter_FCN
void MuonResiduals6DOFFitter_FCN(int &npar, double *gin, double &fval, double *par, int iflag)
Definition: MuonResiduals6DOFFitter.cc:134
void MuonResiduals6DOFFitter::inform ( TMinuit *  tMinuit)
overrideprotectedvirtual

Implements MuonResidualsFitter.

Definition at line 262 of file MuonResiduals6DOFFitter.cc.

262 { minuit = tMinuit; }
int MuonResiduals6DOFFitter::ndata ( )
inlineoverridevirtual

Implements MuonResidualsFitter.

Definition at line 78 of file MuonResiduals6DOFFitter.h.

References kNData.

78 { return kNData; }
int MuonResiduals6DOFFitter::npar ( )
inlineoverridevirtual

Implements MuonResidualsFitter.

Definition at line 67 of file MuonResiduals6DOFFitter.h.

References cms::cuda::assert(), MuonResidualsFitter::kGaussPowerTails, kNPar, MuonResidualsFitter::kPowerLawTails, MuonResidualsFitter::kPureGaussian, MuonResidualsFitter::kPureGaussian2D, MuonResidualsFitter::kROOTVoigt, and MuonResidualsFitter::residualsModel().

Referenced by plot().

67  {
68  if (residualsModel() == kPureGaussian || residualsModel() == kPureGaussian2D ||
69  residualsModel() == kGaussPowerTails)
70  return kNPar - 4;
71  else if (residualsModel() == kPowerLawTails)
72  return kNPar;
73  else if (residualsModel() == kROOTVoigt)
74  return kNPar;
75  else
76  assert(false);
77  }
cms::cuda::assert
assert(be >=bs)
double MuonResiduals6DOFFitter::plot ( std::string  name,
TFileDirectory *  dir,
Alignable *  ali 
)
overridevirtual

Implements MuonResidualsFitter.

Definition at line 443 of file MuonResiduals6DOFFitter.cc.

References cms::cuda::assert(), HLT_FULL_cff::chi2, MuonResidualsFitter::errorerror(), mps_fire::i, kAlignPhiX, kAlignPhiY, kAlignPhiZ, kAlignX, kAlignY, kAlignZ, kAlphaX, kAlphaY, kAngleX, kAngleY, MuonResidualsFitter::kGaussPowerTails, kPositionX, kPositionY, MuonResidualsFitter::kPowerLawTails, MuonResidualsFitter::kPureGaussian, MuonResidualsFitter::kPureGaussian2D, kRedChi2, kResidX, kResidXGamma, kResidXSigma, kResidY, kResidYGamma, kResidYSigma, kResSlopeX, kResSlopeXGamma, kResSlopeXSigma, kResSlopeY, kResSlopeYGamma, kResSlopeYSigma, MuonResidualsFitter::kROOTVoigt, AlignableSurface::length(), MuonResidualsFitter::m_weightAlignment, TFileDirectory::make(), MuonResidualsFitter_GaussPowerTails_TF1(), MuonResidualsFitter_powerLawTails_TF1(), MuonResidualsFitter_pureGaussian_TF1(), MuonResidualsFitter_ROOTVoigt_TF1(), ndof, npar(), funct::pow(), alignCSCRings::r, MuonResidualsFitter::residuals_begin(), MuonResidualsFitter::residuals_end(), MuonResidualsFitter::residualsModel(), AlCaHLTBitMon_QueryRunRegistry::string, sumofweights(), Alignable::surface(), MuonResidualsFitter::value(), histoStyle::weight, and AlignableSurface::width().

443  {
444  sumofweights();
445 
446  double mean_residualx = 0., mean_residualy = 0., mean_resslopex = 0., mean_resslopey = 0.;
447  double mean_trackx = 0., mean_tracky = 0., mean_trackdxdz = 0., mean_trackdydz = 0.;
448  double sum_w = 0.;
449 
450  for (std::vector<double *>::const_iterator rit = residuals_begin(); rit != residuals_end(); ++rit) {
451  const double redchi2 = (*rit)[kRedChi2];
452  double weight = 1. / redchi2;
453  if (!m_weightAlignment)
454  weight = 1.;
455 
456  if (!m_weightAlignment || TMath::Prob(redchi2 * 12, 12) < 0.99) // no spikes allowed
457  {
458  double factor_w = 1. / (sum_w + weight);
459  mean_residualx = factor_w * (sum_w * mean_residualx + weight * (*rit)[kResidX]);
460  mean_residualy = factor_w * (sum_w * mean_residualy + weight * (*rit)[kResidY]);
461  mean_resslopex = factor_w * (sum_w * mean_resslopex + weight * (*rit)[kResSlopeX]);
462  mean_resslopey = factor_w * (sum_w * mean_resslopey + weight * (*rit)[kResSlopeY]);
463  mean_trackx = factor_w * (sum_w * mean_trackx + weight * (*rit)[kPositionX]);
464  mean_tracky = factor_w * (sum_w * mean_tracky + weight * (*rit)[kPositionY]);
465  mean_trackdxdz = factor_w * (sum_w * mean_trackdxdz + weight * (*rit)[kAngleX]);
466  mean_trackdydz = factor_w * (sum_w * mean_trackdydz + weight * (*rit)[kAngleY]);
467  sum_w += weight;
468  }
469  }
470 
471  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,
472  name_y_cut, name_alphax, name_alphay;
473  std::string name_x_trackx, name_y_trackx, name_dxdz_trackx, name_dydz_trackx;
474  std::string name_x_tracky, name_y_tracky, name_dxdz_tracky, name_dydz_tracky;
475  std::string name_x_trackdxdz, name_y_trackdxdz, name_dxdz_trackdxdz, name_dydz_trackdxdz;
476  std::string name_x_trackdydz, name_y_trackdydz, name_dxdz_trackdydz, name_dydz_trackdydz;
477 
478  name_x = name + "_x";
479  name_y = name + "_y";
480  name_dxdz = name + "_dxdz";
481  name_dydz = name + "_dydz";
482  name_x_raw = name + "_x_raw";
483  name_y_raw = name + "_y_raw";
484  name_dxdz_raw = name + "_dxdz_raw";
485  name_dydz_raw = name + "_dydz_raw";
486  name_x_cut = name + "_x_cut";
487  name_y_cut = name + "_y_cut";
488  name_alphax = name + "_alphax";
489  name_alphay = name + "_alphay";
490  name_x_trackx = name + "_x_trackx";
491  name_y_trackx = name + "_y_trackx";
492  name_dxdz_trackx = name + "_dxdz_trackx";
493  name_dydz_trackx = name + "_dydz_trackx";
494  name_x_tracky = name + "_x_tracky";
495  name_y_tracky = name + "_y_tracky";
496  name_dxdz_tracky = name + "_dxdz_tracky";
497  name_dydz_tracky = name + "_dydz_tracky";
498  name_x_trackdxdz = name + "_x_trackdxdz";
499  name_y_trackdxdz = name + "_y_trackdxdz";
500  name_dxdz_trackdxdz = name + "_dxdz_trackdxdz";
501  name_dydz_trackdxdz = name + "_dydz_trackdxdz";
502  name_x_trackdydz = name + "_x_trackdydz";
503  name_y_trackdydz = name + "_y_trackdydz";
504  name_dxdz_trackdydz = name + "_dxdz_trackdydz";
505  name_dydz_trackdydz = name + "_dydz_trackdydz";
506 
507  double width = ali->surface().width();
508  double length = ali->surface().length();
509  int bins = 200;
510  double min_x = -100., max_x = 100.;
511  double min_y = -100., max_y = 100.;
512  double min_dxdz = -75., max_dxdz = 75.;
513  double min_dydz = -150., max_dydz = 150.;
514  double min_trackx = -width / 2., max_trackx = width / 2.;
515  double min_tracky = -length / 2., max_tracky = length / 2.;
516  double min_trackdxdz = -1.5, max_trackdxdz = 1.5;
517  double min_trackdydz = -1.5, max_trackdydz = 1.5;
518 
519  TH1F *hist_x = dir->make<TH1F>(name_x.c_str(), "", bins, min_x, max_x);
520  TH1F *hist_y = dir->make<TH1F>(name_y.c_str(), "", bins, min_y, max_y);
521  TH1F *hist_dxdz = dir->make<TH1F>(name_dxdz.c_str(), "", bins, min_dxdz, max_dxdz);
522  TH1F *hist_dydz = dir->make<TH1F>(name_dydz.c_str(), "", bins, min_dydz, max_dydz);
523  TH1F *hist_x_raw = dir->make<TH1F>(name_x_raw.c_str(), "", bins, min_x, max_x);
524  TH1F *hist_y_raw = dir->make<TH1F>(name_y_raw.c_str(), "", bins, min_y, max_y);
525  TH1F *hist_dxdz_raw = dir->make<TH1F>(name_dxdz_raw.c_str(), "", bins, min_dxdz, max_dxdz);
526  TH1F *hist_dydz_raw = dir->make<TH1F>(name_dydz_raw.c_str(), "", bins, min_dydz, max_dydz);
527  TH1F *hist_x_cut = dir->make<TH1F>(name_x_cut.c_str(), "", bins, min_x, max_x);
528  TH1F *hist_y_cut = dir->make<TH1F>(name_y_cut.c_str(), "", bins, min_y, max_y);
529  TH2F *hist_alphax = dir->make<TH2F>(name_alphax.c_str(), "", 50, 50, 50, 50, -50., 50.);
530  TH2F *hist_alphay = dir->make<TH2F>(name_alphay.c_str(), "", 75, 100, 100, 75, -75., 75.);
531 
532  TProfile *hist_x_trackx = dir->make<TProfile>(name_x_trackx.c_str(), "", 50, min_trackx, max_trackx, min_x, max_x);
533  TProfile *hist_y_trackx = dir->make<TProfile>(name_y_trackx.c_str(), "", 50, min_trackx, max_trackx, min_y, max_y);
534  TProfile *hist_dxdz_trackx =
535  dir->make<TProfile>(name_dxdz_trackx.c_str(), "", 50, min_trackx, max_trackx, min_dxdz, max_dxdz);
536  TProfile *hist_dydz_trackx =
537  dir->make<TProfile>(name_dydz_trackx.c_str(), "", 50, min_trackx, max_trackx, min_dydz, max_dydz);
538  TProfile *hist_x_tracky = dir->make<TProfile>(name_x_tracky.c_str(), "", 50, min_tracky, max_tracky, min_x, max_x);
539  TProfile *hist_y_tracky = dir->make<TProfile>(name_y_tracky.c_str(), "", 50, min_tracky, max_tracky, min_y, max_y);
540  TProfile *hist_dxdz_tracky =
541  dir->make<TProfile>(name_dxdz_tracky.c_str(), "", 50, min_tracky, max_tracky, min_dxdz, max_dxdz);
542  TProfile *hist_dydz_tracky =
543  dir->make<TProfile>(name_dydz_tracky.c_str(), "", 50, min_tracky, max_tracky, min_dydz, max_dydz);
544  TProfile *hist_x_trackdxdz =
545  dir->make<TProfile>(name_x_trackdxdz.c_str(), "", 250, min_trackdxdz, max_trackdxdz, min_x, max_x);
546  TProfile *hist_y_trackdxdz =
547  dir->make<TProfile>(name_y_trackdxdz.c_str(), "", 250, min_trackdxdz, max_trackdxdz, min_y, max_y);
548  TProfile *hist_dxdz_trackdxdz =
549  dir->make<TProfile>(name_dxdz_trackdxdz.c_str(), "", 250, min_trackdxdz, max_trackdxdz, min_dxdz, max_dxdz);
550  TProfile *hist_dydz_trackdxdz =
551  dir->make<TProfile>(name_dydz_trackdxdz.c_str(), "", 250, min_trackdxdz, max_trackdxdz, min_dydz, max_dydz);
552  TProfile *hist_x_trackdydz =
553  dir->make<TProfile>(name_x_trackdydz.c_str(), "", 250, min_trackdydz, max_trackdydz, min_x, max_x);
554  TProfile *hist_y_trackdydz =
555  dir->make<TProfile>(name_y_trackdydz.c_str(), "", 250, min_trackdydz, max_trackdydz, min_y, max_y);
556  TProfile *hist_dxdz_trackdydz =
557  dir->make<TProfile>(name_dxdz_trackdydz.c_str(), "", 250, min_trackdydz, max_trackdydz, min_dxdz, max_dxdz);
558  TProfile *hist_dydz_trackdydz =
559  dir->make<TProfile>(name_dydz_trackdydz.c_str(), "", 250, min_trackdydz, max_trackdydz, min_dydz, max_dydz);
560 
561  hist_x_trackx->SetAxisRange(-10., 10., "Y");
562  hist_y_trackx->SetAxisRange(-20., 20., "Y");
563  hist_dxdz_trackx->SetAxisRange(-10., 10., "Y");
564  hist_dydz_trackx->SetAxisRange(-20., 20., "Y");
565  hist_x_tracky->SetAxisRange(-10., 10., "Y");
566  hist_y_tracky->SetAxisRange(-20., 20., "Y");
567  hist_dxdz_tracky->SetAxisRange(-10., 10., "Y");
568  hist_dydz_tracky->SetAxisRange(-20., 20., "Y");
569  hist_x_trackdxdz->SetAxisRange(-10., 10., "Y");
570  hist_y_trackdxdz->SetAxisRange(-20., 20., "Y");
571  hist_dxdz_trackdxdz->SetAxisRange(-10., 10., "Y");
572  hist_dydz_trackdxdz->SetAxisRange(-20., 20., "Y");
573  hist_x_trackdydz->SetAxisRange(-10., 10., "Y");
574  hist_y_trackdydz->SetAxisRange(-20., 20., "Y");
575  hist_dxdz_trackdydz->SetAxisRange(-10., 10., "Y");
576  hist_dydz_trackdydz->SetAxisRange(-20., 20., "Y");
577 
578  name_x += "_fit";
579  name_y += "_fit";
580  name_dxdz += "_fit";
581  name_dydz += "_fit";
582  name_alphax += "_fit";
583  name_alphay += "_fit";
584  name_x_trackx += "_fit";
585  name_y_trackx += "_fit";
586  name_dxdz_trackx += "_fit";
587  name_dydz_trackx += "_fit";
588  name_x_tracky += "_fit";
589  name_y_tracky += "_fit";
590  name_dxdz_tracky += "_fit";
591  name_dydz_tracky += "_fit";
592  name_x_trackdxdz += "_fit";
593  name_y_trackdxdz += "_fit";
594  name_dxdz_trackdxdz += "_fit";
595  name_dydz_trackdxdz += "_fit";
596  name_x_trackdydz += "_fit";
597  name_y_trackdydz += "_fit";
598  name_dxdz_trackdydz += "_fit";
599  name_dydz_trackdydz += "_fit";
600 
601  TF1 *fit_x = nullptr;
602  TF1 *fit_y = nullptr;
603  TF1 *fit_dxdz = nullptr;
604  TF1 *fit_dydz = nullptr;
605  if (residualsModel() == kPureGaussian || residualsModel() == kPureGaussian2D) {
606  fit_x = new TF1(name_x.c_str(), MuonResidualsFitter_pureGaussian_TF1, min_x, max_x, 3);
607  fit_x->SetParameters(sum_of_weights * (max_x - min_x) / bins, 10. * value(kAlignX), 10. * value(kResidXSigma));
608  const double er_x[3] = {0., 10. * errorerror(kAlignX), 10. * errorerror(kResidXSigma)};
609  fit_x->SetParErrors(er_x);
610  fit_y = new TF1(name_y.c_str(), MuonResidualsFitter_pureGaussian_TF1, min_y, max_y, 3);
611  fit_y->SetParameters(sum_of_weights * (max_y - min_y) / bins, 10. * value(kAlignY), 10. * value(kResidYSigma));
612  const double er_y[3] = {0., 10. * errorerror(kAlignY), 10. * errorerror(kResidYSigma)};
613  fit_y->SetParErrors(er_y);
614  fit_dxdz = new TF1(name_dxdz.c_str(), MuonResidualsFitter_pureGaussian_TF1, min_dxdz, max_dxdz, 3);
615  fit_dxdz->SetParameters(
616  sum_of_weights * (max_dxdz - min_dxdz) / bins, 1000. * value(kAlignPhiY), 1000. * value(kResSlopeXSigma));
617  const double er_dxdz[3] = {0., 1000. * errorerror(kAlignPhiY), 1000. * errorerror(kResSlopeXSigma)};
618  fit_dxdz->SetParErrors(er_dxdz);
619  fit_dydz = new TF1(name_dydz.c_str(), MuonResidualsFitter_pureGaussian_TF1, min_dydz, max_dydz, 3);
620  fit_dydz->SetParameters(
621  sum_of_weights * (max_dydz - min_dydz) / bins, -1000. * value(kAlignPhiX), 1000. * value(kResSlopeYSigma));
622  const double er_dydz[3] = {0., 1000. * errorerror(kAlignPhiX), 1000. * errorerror(kResSlopeYSigma)};
623  fit_dydz->SetParErrors(er_dydz);
624  } else if (residualsModel() == kPowerLawTails) {
625  fit_x = new TF1(name_x.c_str(), MuonResidualsFitter_powerLawTails_TF1, min_x, max_x, 4);
626  fit_x->SetParameters(sum_of_weights * (max_x - min_x) / bins,
627  10. * value(kAlignX),
628  10. * value(kResidXSigma),
629  10. * value(kResidXGamma));
630  fit_y = new TF1(name_y.c_str(), MuonResidualsFitter_powerLawTails_TF1, min_y, max_y, 4);
631  fit_y->SetParameters(sum_of_weights * (max_y - min_y) / bins,
632  10. * value(kAlignY),
633  10. * value(kResidYSigma),
634  10. * value(kResidYGamma));
635  fit_dxdz = new TF1(name_dxdz.c_str(), MuonResidualsFitter_powerLawTails_TF1, min_dxdz, max_dxdz, 4);
636  fit_dxdz->SetParameters(sum_of_weights * (max_dxdz - min_dxdz) / bins,
637  1000. * value(kAlignPhiY),
638  1000. * value(kResSlopeXSigma),
639  1000. * value(kResSlopeXGamma));
640  fit_dydz = new TF1(name_dydz.c_str(), MuonResidualsFitter_powerLawTails_TF1, min_dydz, max_dydz, 4);
641  fit_dydz->SetParameters(sum_of_weights * (max_dydz - min_dydz) / bins,
642  -1000. * value(kAlignPhiX),
643  1000. * value(kResSlopeYSigma),
644  1000. * value(kResSlopeYGamma));
645  } else if (residualsModel() == kROOTVoigt) {
646  fit_x = new TF1(name_x.c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_x, max_x, 4);
647  fit_x->SetParameters(sum_of_weights * (max_x - min_x) / bins,
648  10. * value(kAlignX),
649  10. * value(kResidXSigma),
650  10. * value(kResidXGamma));
651  fit_y = new TF1(name_y.c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_y, max_y, 4);
652  fit_y->SetParameters(sum_of_weights * (max_y - min_y) / bins,
653  10. * value(kAlignY),
654  10. * value(kResidYSigma),
655  10. * value(kResidYGamma));
656  fit_dxdz = new TF1(name_dxdz.c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_dxdz, max_dxdz, 4);
657  fit_dxdz->SetParameters(sum_of_weights * (max_dxdz - min_dxdz) / bins,
658  1000. * value(kAlignPhiY),
659  1000. * value(kResSlopeXSigma),
660  1000. * value(kResSlopeXGamma));
661  fit_dydz = new TF1(name_dydz.c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_dydz, max_dydz, 4);
662  fit_dydz->SetParameters(sum_of_weights * (max_dydz - min_dydz) / bins,
663  -1000. * value(kAlignPhiX),
664  1000. * value(kResSlopeYSigma),
665  1000. * value(kResSlopeYGamma));
666  } else if (residualsModel() == kGaussPowerTails) {
667  fit_x = new TF1(name_x.c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_x, max_x, 3);
668  fit_x->SetParameters(sum_of_weights * (max_x - min_x) / bins, 10. * value(kAlignX), 10. * value(kResidXSigma));
669  fit_y = new TF1(name_y.c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_y, max_y, 3);
670  fit_y->SetParameters(sum_of_weights * (max_y - min_y) / bins, 10. * value(kAlignY), 10. * value(kResidYSigma));
671  fit_dxdz = new TF1(name_dxdz.c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_dxdz, max_dxdz, 3);
672  fit_dxdz->SetParameters(
673  sum_of_weights * (max_dxdz - min_dxdz) / bins, 1000. * value(kAlignPhiY), 1000. * value(kResSlopeXSigma));
674  fit_dydz = new TF1(name_dydz.c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_dydz, max_dydz, 3);
675  fit_dydz->SetParameters(
676  sum_of_weights * (max_dydz - min_dydz) / bins, -1000. * value(kAlignPhiX), 1000. * value(kResSlopeYSigma));
677  } else {
678  assert(false);
679  }
680 
681  fit_x->SetLineColor(2);
682  fit_x->SetLineWidth(2);
683  fit_x->Write();
684  fit_y->SetLineColor(2);
685  fit_y->SetLineWidth(2);
686  fit_y->Write();
687  fit_dxdz->SetLineColor(2);
688  fit_dxdz->SetLineWidth(2);
689  fit_dxdz->Write();
690  fit_dydz->SetLineColor(2);
691  fit_dydz->SetLineWidth(2);
692  fit_dydz->Write();
693 
694  TF1 *fit_alphax = new TF1(name_alphax.c_str(), "[0] + x*[1]", min_dxdz, max_dxdz);
695  TF1 *fit_alphay = new TF1(name_alphay.c_str(), "[0] + x*[1]", min_dydz, max_dydz);
696  double aX = 10. * value(kAlignX), bX = 10. * value(kAlphaX) / 1000.;
697  double aY = 10. * value(kAlignY), bY = 10. * value(kAlphaY) / 1000.;
698  if (residualsModel() == kPureGaussian2D) {
699  double sx = 10. * value(kResidXSigma), sy = 1000. * value(kResSlopeXSigma), r = value(kAlphaX);
700  aX = mean_residualx;
701  bX = 0.;
702  if (sx != 0.) {
703  bX = 1. / (sy / sx * r);
704  aX = -bX * mean_resslopex;
705  }
706  sx = 10. * value(kResidYSigma);
707  sy = 1000. * value(kResSlopeYSigma);
708  r = value(kAlphaY);
709  aY = mean_residualx;
710  bY = 0.;
711  if (sx != 0.) {
712  bY = 1. / (sy / sx * r);
713  aY = -bY * mean_resslopey;
714  }
715  }
716  fit_alphax->SetParameters(aX, bX);
717  fit_alphay->SetParameters(aY, bY);
718 
719  fit_alphax->SetLineColor(2);
720  fit_alphax->SetLineWidth(2);
721  fit_alphax->Write();
722  fit_alphay->SetLineColor(2);
723  fit_alphay->SetLineWidth(2);
724  fit_alphay->Write();
725 
726  TProfile *fit_x_trackx = dir->make<TProfile>(name_x_trackx.c_str(), "", 100, min_trackx, max_trackx);
727  TProfile *fit_y_trackx = dir->make<TProfile>(name_y_trackx.c_str(), "", 100, min_trackx, max_trackx);
728  TProfile *fit_dxdz_trackx = dir->make<TProfile>(name_dxdz_trackx.c_str(), "", 100, min_trackx, max_trackx);
729  TProfile *fit_dydz_trackx = dir->make<TProfile>(name_dydz_trackx.c_str(), "", 100, min_trackx, max_trackx);
730  TProfile *fit_x_tracky = dir->make<TProfile>(name_x_tracky.c_str(), "", 100, min_tracky, max_tracky);
731  TProfile *fit_y_tracky = dir->make<TProfile>(name_y_tracky.c_str(), "", 100, min_tracky, max_tracky);
732  TProfile *fit_dxdz_tracky = dir->make<TProfile>(name_dxdz_tracky.c_str(), "", 100, min_tracky, max_tracky);
733  TProfile *fit_dydz_tracky = dir->make<TProfile>(name_dydz_tracky.c_str(), "", 100, min_tracky, max_tracky);
734  TProfile *fit_x_trackdxdz = dir->make<TProfile>(name_x_trackdxdz.c_str(), "", 500, min_trackdxdz, max_trackdxdz);
735  TProfile *fit_y_trackdxdz = dir->make<TProfile>(name_y_trackdxdz.c_str(), "", 500, min_trackdxdz, max_trackdxdz);
736  TProfile *fit_dxdz_trackdxdz =
737  dir->make<TProfile>(name_dxdz_trackdxdz.c_str(), "", 500, min_trackdxdz, max_trackdxdz);
738  TProfile *fit_dydz_trackdxdz =
739  dir->make<TProfile>(name_dydz_trackdxdz.c_str(), "", 500, min_trackdxdz, max_trackdxdz);
740  TProfile *fit_x_trackdydz = dir->make<TProfile>(name_x_trackdydz.c_str(), "", 500, min_trackdydz, max_trackdydz);
741  TProfile *fit_y_trackdydz = dir->make<TProfile>(name_y_trackdydz.c_str(), "", 500, min_trackdydz, max_trackdydz);
742  TProfile *fit_dxdz_trackdydz =
743  dir->make<TProfile>(name_dxdz_trackdydz.c_str(), "", 500, min_trackdydz, max_trackdydz);
744  TProfile *fit_dydz_trackdydz =
745  dir->make<TProfile>(name_dydz_trackdydz.c_str(), "", 500, min_trackdydz, max_trackdydz);
746 
747  fit_x_trackx->SetLineColor(2);
748  fit_x_trackx->SetLineWidth(2);
749  fit_y_trackx->SetLineColor(2);
750  fit_y_trackx->SetLineWidth(2);
751  fit_dxdz_trackx->SetLineColor(2);
752  fit_dxdz_trackx->SetLineWidth(2);
753  fit_dydz_trackx->SetLineColor(2);
754  fit_dydz_trackx->SetLineWidth(2);
755  fit_x_tracky->SetLineColor(2);
756  fit_x_tracky->SetLineWidth(2);
757  fit_y_tracky->SetLineColor(2);
758  fit_y_tracky->SetLineWidth(2);
759  fit_dxdz_tracky->SetLineColor(2);
760  fit_dxdz_tracky->SetLineWidth(2);
761  fit_dydz_tracky->SetLineColor(2);
762  fit_dydz_tracky->SetLineWidth(2);
763  fit_x_trackdxdz->SetLineColor(2);
764  fit_x_trackdxdz->SetLineWidth(2);
765  fit_y_trackdxdz->SetLineColor(2);
766  fit_y_trackdxdz->SetLineWidth(2);
767  fit_dxdz_trackdxdz->SetLineColor(2);
768  fit_dxdz_trackdxdz->SetLineWidth(2);
769  fit_dydz_trackdxdz->SetLineColor(2);
770  fit_dydz_trackdxdz->SetLineWidth(2);
771  fit_x_trackdydz->SetLineColor(2);
772  fit_x_trackdydz->SetLineWidth(2);
773  fit_y_trackdydz->SetLineColor(2);
774  fit_y_trackdydz->SetLineWidth(2);
775  fit_dxdz_trackdydz->SetLineColor(2);
776  fit_dxdz_trackdydz->SetLineWidth(2);
777  fit_dydz_trackdydz->SetLineColor(2);
778  fit_dydz_trackdydz->SetLineWidth(2);
779 
780  name_x_trackx += "line";
781  name_y_trackx += "line";
782  name_dxdz_trackx += "line";
783  name_dydz_trackx += "line";
784  name_x_tracky += "line";
785  name_y_tracky += "line";
786  name_dxdz_tracky += "line";
787  name_dydz_tracky += "line";
788  name_x_trackdxdz += "line";
789  name_y_trackdxdz += "line";
790  name_dxdz_trackdxdz += "line";
791  name_dydz_trackdxdz += "line";
792  name_x_trackdydz += "line";
793  name_y_trackdydz += "line";
794  name_dxdz_trackdydz += "line";
795  name_dydz_trackdydz += "line";
796 
797  TF1 *fitline_x_trackx = new TF1(name_x_trackx.c_str(), residual_x_trackx_TF1, min_trackx, max_trackx, 14);
798  TF1 *fitline_y_trackx = new TF1(name_y_trackx.c_str(), residual_y_trackx_TF1, min_trackx, max_trackx, 14);
799  TF1 *fitline_dxdz_trackx = new TF1(name_dxdz_trackx.c_str(), residual_dxdz_trackx_TF1, min_trackx, max_trackx, 14);
800  TF1 *fitline_dydz_trackx = new TF1(name_dydz_trackx.c_str(), residual_dydz_trackx_TF1, min_trackx, max_trackx, 14);
801  TF1 *fitline_x_tracky = new TF1(name_x_tracky.c_str(), residual_x_tracky_TF1, min_tracky, max_tracky, 14);
802  TF1 *fitline_y_tracky = new TF1(name_y_tracky.c_str(), residual_y_tracky_TF1, min_tracky, max_tracky, 14);
803  TF1 *fitline_dxdz_tracky = new TF1(name_dxdz_tracky.c_str(), residual_dxdz_tracky_TF1, min_tracky, max_tracky, 14);
804  TF1 *fitline_dydz_tracky = new TF1(name_dydz_tracky.c_str(), residual_dydz_tracky_TF1, min_tracky, max_tracky, 14);
805  TF1 *fitline_x_trackdxdz =
806  new TF1(name_x_trackdxdz.c_str(), residual_x_trackdxdz_TF1, min_trackdxdz, max_trackdxdz, 14);
807  TF1 *fitline_y_trackdxdz =
808  new TF1(name_y_trackdxdz.c_str(), residual_y_trackdxdz_TF1, min_trackdxdz, max_trackdxdz, 14);
809  TF1 *fitline_dxdz_trackdxdz =
810  new TF1(name_dxdz_trackdxdz.c_str(), residual_dxdz_trackdxdz_TF1, min_trackdxdz, max_trackdxdz, 14);
811  TF1 *fitline_dydz_trackdxdz =
812  new TF1(name_dydz_trackdxdz.c_str(), residual_dydz_trackdxdz_TF1, min_trackdxdz, max_trackdxdz, 14);
813  TF1 *fitline_x_trackdydz =
814  new TF1(name_x_trackdydz.c_str(), residual_x_trackdydz_TF1, min_trackdydz, max_trackdydz, 14);
815  TF1 *fitline_y_trackdydz =
816  new TF1(name_y_trackdydz.c_str(), residual_y_trackdydz_TF1, min_trackdydz, max_trackdydz, 14);
817  TF1 *fitline_dxdz_trackdydz =
818  new TF1(name_dxdz_trackdydz.c_str(), residual_dxdz_trackdydz_TF1, min_trackdydz, max_trackdydz, 14);
819  TF1 *fitline_dydz_trackdydz =
820  new TF1(name_dydz_trackdydz.c_str(), residual_dydz_trackdydz_TF1, min_trackdydz, max_trackdydz, 14);
821 
822  std::vector<TF1 *> fitlines;
823  fitlines.push_back(fitline_x_trackx);
824  fitlines.push_back(fitline_y_trackx);
825  fitlines.push_back(fitline_dxdz_trackx);
826  fitlines.push_back(fitline_dydz_trackx);
827  fitlines.push_back(fitline_x_tracky);
828  fitlines.push_back(fitline_y_tracky);
829  fitlines.push_back(fitline_dxdz_tracky);
830  fitlines.push_back(fitline_dydz_tracky);
831  fitlines.push_back(fitline_x_trackdxdz);
832  fitlines.push_back(fitline_y_trackdxdz);
833  fitlines.push_back(fitline_dxdz_trackdxdz);
834  fitlines.push_back(fitline_dydz_trackdxdz);
835  fitlines.push_back(fitline_x_trackdydz);
836  fitlines.push_back(fitline_y_trackdydz);
837  fitlines.push_back(fitline_dxdz_trackdydz);
838  fitlines.push_back(fitline_dydz_trackdydz);
839 
840  double fitparameters[14] = {value(kAlignX),
841  value(kAlignY),
842  value(kAlignZ),
843  value(kAlignPhiX),
844  value(kAlignPhiY),
845  value(kAlignPhiZ),
846  mean_trackx,
847  mean_tracky,
848  mean_trackdxdz,
849  mean_trackdydz,
850  value(kAlphaX),
851  mean_resslopex,
852  value(kAlphaY),
853  mean_resslopey};
854  if (residualsModel() == kPureGaussian2D)
855  fitparameters[10] = fitparameters[12] = 0.;
856 
857  for (std::vector<TF1 *>::const_iterator itr = fitlines.begin(); itr != fitlines.end(); itr++) {
858  (*itr)->SetParameters(fitparameters);
859  (*itr)->SetLineColor(2);
860  (*itr)->SetLineWidth(2);
861  (*itr)->Write();
862  }
863 
864  for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
865  const double residX = (*resiter)[kResidX];
866  const double residY = (*resiter)[kResidY];
867  const double resslopeX = (*resiter)[kResSlopeX];
868  const double resslopeY = (*resiter)[kResSlopeY];
869  const double positionX = (*resiter)[kPositionX];
870  const double positionY = (*resiter)[kPositionY];
871  const double angleX = (*resiter)[kAngleX];
872  const double angleY = (*resiter)[kAngleY];
873  const double redchi2 = (*resiter)[kRedChi2];
874  double weight = 1. / redchi2;
875  if (!m_weightAlignment)
876  weight = 1.;
877 
878  if (!m_weightAlignment || TMath::Prob(redchi2 * 12, 12) < 0.99) { // no spikes allowed
879 
880  hist_alphax->Fill(1000. * resslopeX, 10. * residX);
881  hist_alphay->Fill(1000. * resslopeY, 10. * residY);
882 
883  double coefX = value(kAlphaX), coefY = value(kAlphaY);
884  if (residualsModel() == kPureGaussian || residualsModel() == kPureGaussian2D)
885  coefX = coefY = 0.;
886  double geom_residX = residual_x(value(kAlignX),
887  value(kAlignY),
888  value(kAlignZ),
889  value(kAlignPhiX),
890  value(kAlignPhiY),
891  value(kAlignPhiZ),
892  positionX,
893  positionY,
894  angleX,
895  angleY,
896  coefX,
897  resslopeX);
898  hist_x->Fill(10. * (residX - geom_residX + value(kAlignX)), weight);
899  hist_x_trackx->Fill(positionX, 10. * residX, weight);
900  hist_x_tracky->Fill(positionY, 10. * residX, weight);
901  hist_x_trackdxdz->Fill(angleX, 10. * residX, weight);
902  hist_x_trackdydz->Fill(angleY, 10. * residX, weight);
903  fit_x_trackx->Fill(positionX, 10. * geom_residX, weight);
904  fit_x_tracky->Fill(positionY, 10. * geom_residX, weight);
905  fit_x_trackdxdz->Fill(angleX, 10. * geom_residX, weight);
906  fit_x_trackdydz->Fill(angleY, 10. * geom_residX, weight);
907 
908  double geom_residY = residual_y(value(kAlignX),
909  value(kAlignY),
910  value(kAlignZ),
911  value(kAlignPhiX),
912  value(kAlignPhiY),
913  value(kAlignPhiZ),
914  positionX,
915  positionY,
916  angleX,
917  angleY,
918  coefY,
919  resslopeY);
920  hist_y->Fill(10. * (residY - geom_residY + value(kAlignY)), weight);
921  hist_y_trackx->Fill(positionX, 10. * residY, weight);
922  hist_y_tracky->Fill(positionY, 10. * residY, weight);
923  hist_y_trackdxdz->Fill(angleX, 10. * residY, weight);
924  hist_y_trackdydz->Fill(angleY, 10. * residY, weight);
925  fit_y_trackx->Fill(positionX, 10. * geom_residY, weight);
926  fit_y_tracky->Fill(positionY, 10. * geom_residY, weight);
927  fit_y_trackdxdz->Fill(angleX, 10. * geom_residY, weight);
928  fit_y_trackdydz->Fill(angleY, 10. * geom_residY, weight);
929 
930  double geom_resslopeX = residual_dxdz(value(kAlignX),
931  value(kAlignY),
932  value(kAlignZ),
933  value(kAlignPhiX),
934  value(kAlignPhiY),
935  value(kAlignPhiZ),
936  positionX,
937  positionY,
938  angleX,
939  angleY);
940  hist_dxdz->Fill(1000. * (resslopeX - geom_resslopeX + value(kAlignPhiY)), weight);
941  hist_dxdz_trackx->Fill(positionX, 1000. * resslopeX, weight);
942  hist_dxdz_tracky->Fill(positionY, 1000. * resslopeX, weight);
943  hist_dxdz_trackdxdz->Fill(angleX, 1000. * resslopeX, weight);
944  hist_dxdz_trackdydz->Fill(angleY, 1000. * resslopeX, weight);
945  fit_dxdz_trackx->Fill(positionX, 1000. * geom_resslopeX, weight);
946  fit_dxdz_tracky->Fill(positionY, 1000. * geom_resslopeX, weight);
947  fit_dxdz_trackdxdz->Fill(angleX, 1000. * geom_resslopeX, weight);
948  fit_dxdz_trackdydz->Fill(angleY, 1000. * geom_resslopeX, weight);
949 
950  double geom_resslopeY = residual_dydz(value(kAlignX),
951  value(kAlignY),
952  value(kAlignZ),
953  value(kAlignPhiX),
954  value(kAlignPhiY),
955  value(kAlignPhiZ),
956  positionX,
957  positionY,
958  angleX,
959  angleY);
960  hist_dydz->Fill(1000. * (resslopeY - geom_resslopeY - value(kAlignPhiX)), weight);
961  hist_dydz_trackx->Fill(positionX, 1000. * resslopeY, weight);
962  hist_dydz_tracky->Fill(positionY, 1000. * resslopeY, weight);
963  hist_dydz_trackdxdz->Fill(angleX, 1000. * resslopeY, weight);
964  hist_dydz_trackdydz->Fill(angleY, 1000. * resslopeY, weight);
965  fit_dydz_trackx->Fill(positionX, 1000. * geom_resslopeY, weight);
966  fit_dydz_tracky->Fill(positionY, 1000. * geom_resslopeY, weight);
967  fit_dydz_trackdxdz->Fill(angleX, 1000. * geom_resslopeY, weight);
968  fit_dydz_trackdydz->Fill(angleY, 1000. * geom_resslopeY, weight);
969  }
970 
971  hist_x_raw->Fill(10. * residX);
972  hist_y_raw->Fill(10. * residY);
973  hist_dxdz_raw->Fill(1000. * resslopeX);
974  hist_dydz_raw->Fill(1000. * resslopeY);
975  if (fabs(resslopeX) < 0.005)
976  hist_x_cut->Fill(10. * residX);
977  if (fabs(resslopeY) < 0.030)
978  hist_y_cut->Fill(10. * residY);
979  }
980 
981  double chi2 = 0.;
982  double ndof = 0.;
983  for (int i = 1; i <= hist_x->GetNbinsX(); i++) {
984  double xi = hist_x->GetBinCenter(i);
985  double yi = hist_x->GetBinContent(i);
986  double yerri = hist_x->GetBinError(i);
987  double yth = fit_x->Eval(xi);
988  if (yerri > 0.) {
989  chi2 += pow((yth - yi) / yerri, 2);
990  ndof += 1.;
991  }
992  }
993  for (int i = 1; i <= hist_y->GetNbinsX(); i++) {
994  double xi = hist_y->GetBinCenter(i);
995  double yi = hist_y->GetBinContent(i);
996  double yerri = hist_y->GetBinError(i);
997  double yth = fit_y->Eval(xi);
998  if (yerri > 0.) {
999  chi2 += pow((yth - yi) / yerri, 2);
1000  ndof += 1.;
1001  }
1002  }
1003  for (int i = 1; i <= hist_dxdz->GetNbinsX(); i++) {
1004  double xi = hist_dxdz->GetBinCenter(i);
1005  double yi = hist_dxdz->GetBinContent(i);
1006  double yerri = hist_dxdz->GetBinError(i);
1007  double yth = fit_dxdz->Eval(xi);
1008  if (yerri > 0.) {
1009  chi2 += pow((yth - yi) / yerri, 2);
1010  ndof += 1.;
1011  }
1012  }
1013  for (int i = 1; i <= hist_dydz->GetNbinsX(); i++) {
1014  double xi = hist_dydz->GetBinCenter(i);
1015  double yi = hist_dydz->GetBinContent(i);
1016  double yerri = hist_dydz->GetBinError(i);
1017  double yth = fit_dydz->Eval(xi);
1018  if (yerri > 0.) {
1019  chi2 += pow((yth - yi) / yerri, 2);
1020  ndof += 1.;
1021  }
1022  }
1023  ndof -= npar();
1024 
1025  return (ndof > 0. ? chi2 / ndof : -1.);
1026 }
AlignableSurface::width
align::Scalar width() const
Definition: AlignableSurface.h:30
MuonResidualsFitter::residuals_begin
std::vector< double * >::const_iterator residuals_begin() const
Definition: MuonResidualsFitter.h:190
mps_fire.i
i
Definition: mps_fire.py:428
MuonResidualsFitter_powerLawTails_TF1
Double_t MuonResidualsFitter_powerLawTails_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResidualsFitter.cc:116
MuonResidualsFitter::value
double value(int parNum)
Definition: MuonResidualsFitter.h:145
cms::cuda::assert
assert(be >=bs)
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
HLT_FULL_cff.chi2
tuple chi2
Definition: HLT_FULL_cff.py:9425
MuonResidualsFitter_GaussPowerTails_TF1
Double_t MuonResidualsFitter_GaussPowerTails_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResidualsFitter.cc:142
TFileDirectory::make
T * make(const Args &...args) const
make new ROOT object
Definition: TFileDirectory.h:53
Alignable::surface
const AlignableSurface & surface() const
Return the Surface (global position and orientation) of the object.
Definition: Alignable.h:132
MuonResidualsFitter_pureGaussian_TF1
Double_t MuonResidualsFitter_pureGaussian_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResidualsFitter.cc:38
AlignableSurface::length
align::Scalar length() const
Definition: AlignableSurface.h:32
MuonResidualsFitter::errorerror
double errorerror(int parNum)
Definition: MuonResidualsFitter.h:149
MuonResidualsFitter::residuals_end
std::vector< double * >::const_iterator residuals_end() const
Definition: MuonResidualsFitter.h:191
histoStyle.weight
int weight
Definition: histoStyle.py:51
MuonResidualsFitter_ROOTVoigt_TF1
Double_t MuonResidualsFitter_ROOTVoigt_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResidualsFitter.cc:124
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
TTree * MuonResiduals6DOFFitter::readNtuple ( std::string  fname,
unsigned int  wheel,
unsigned int  station,
unsigned int  sector,
unsigned int  preselected = 1 
)

Definition at line 1028 of file MuonResiduals6DOFFitter.cc.

References MuonResidualsFitter::MuonAlignmentTreeRow::angle_x, MuonResidualsFitter::MuonAlignmentTreeRow::angle_y, cms::cuda::assert(), validate-o2o-wbm::f, MuonResidualsFitter::fill(), mps_fire::i, kAngleX, kAngleY, kCharge, kNData, kPositionX, kPositionY, kPt, kPz, kRedChi2, kResidX, kResidY, kResSlopeX, kResSlopeY, MuonResidualsFitter::MuonAlignmentTreeRow::pos_x, MuonResidualsFitter::MuonAlignmentTreeRow::pos_y, MuonResidualsFitter::MuonAlignmentTreeRow::pt, MuonResidualsFitter::MuonAlignmentTreeRow::pz, MuonResidualsFitter::MuonAlignmentTreeRow::q, alignCSCRings::r, MuonResidualsFitter::MuonAlignmentTreeRow::res_slope_x, MuonResidualsFitter::MuonAlignmentTreeRow::res_slope_y, MuonResidualsFitter::MuonAlignmentTreeRow::res_x, MuonResidualsFitter::MuonAlignmentTreeRow::res_y, submitPVValidationJobs::t, and groupFilesInBlocks::tt.

1029  {
1030  TFile *f = new TFile(fname.c_str());
1031  TTree *t = (TTree *)f->Get("mual_ttree");
1032 
1033  // Create new temporary file
1034  TFile *tmpf = new TFile("small_tree_fit.root", "recreate");
1035  assert(tmpf != nullptr);
1036 
1037  // filter the tree (temporarily lives in the new file)
1038  TTree *tt = t->CopyTree(Form(
1039  "is_dt && ring_wheel==%d && station==%d && sector==%d && select==%d", wheel, station, sector, (bool)preselected));
1040 
1041  MuonAlignmentTreeRow r;
1042  tt->SetBranchAddress("res_x", &r.res_x);
1043  tt->SetBranchAddress("res_slope_x", &r.res_slope_x);
1044  tt->SetBranchAddress("res_y", &r.res_y);
1045  tt->SetBranchAddress("res_slope_y", &r.res_slope_y);
1046  tt->SetBranchAddress("pos_x", &r.pos_x);
1047  tt->SetBranchAddress("pos_y", &r.pos_y);
1048  tt->SetBranchAddress("angle_x", &r.angle_x);
1049  tt->SetBranchAddress("angle_y", &r.angle_y);
1050  tt->SetBranchAddress("pz", &r.pz);
1051  tt->SetBranchAddress("pt", &r.pt);
1052  tt->SetBranchAddress("q", &r.q);
1053 
1054  Long64_t nentries = tt->GetEntries();
1055  for (Long64_t i = 0; i < nentries; i++) {
1056  tt->GetEntry(i);
1057  double *rdata = new double[MuonResiduals6DOFFitter::kNData];
1058  rdata[kResidX] = r.res_x;
1059  rdata[kResidY] = r.res_y;
1060  rdata[kResSlopeX] = r.res_slope_x;
1061  rdata[kResSlopeY] = r.res_slope_y;
1062  rdata[kPositionX] = r.pos_x;
1063  rdata[kPositionY] = r.pos_y;
1064  rdata[kAngleX] = r.angle_x;
1065  rdata[kAngleY] = r.angle_y;
1066  rdata[kRedChi2] = 0.1;
1067  rdata[kPz] = r.pz;
1068  rdata[kPt] = r.pt;
1069  rdata[kCharge] = r.q;
1070  fill(rdata);
1071  }
1072  delete f;
1073  //delete tmpf;
1074  return tt;
1075 }
mps_fire.i
i
Definition: mps_fire.py:428
MuonResidualsFitter::fill
void fill(double *residual)
Definition: MuonResidualsFitter.cc:186
cms::cuda::assert
assert(be >=bs)
alignmentValidation.fname
string fname
main script
Definition: alignmentValidation.py:959
double MuonResiduals6DOFFitter::sumofweights ( )
overridevirtual

Implements MuonResidualsFitter.

Definition at line 266 of file MuonResiduals6DOFFitter.cc.

References kRedChi2, MuonResidualsFitter::m_weightAlignment, MuonResidualsFitter::residuals_begin(), and MuonResidualsFitter::residuals_end().

Referenced by fit(), and plot().

266  {
267  sum_of_weights = 0.;
268  number_of_hits = 0.;
269  weight_alignment = m_weightAlignment;
270  for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
271  if (m_weightAlignment) {
272  double redchi2 = (*resiter)[kRedChi2];
273  if (TMath::Prob(redchi2 * 12, 12) < 0.99) { // no spikes allowed
274  sum_of_weights += 1. / redchi2;
275  number_of_hits += 1.;
276  }
277  } else {
278  sum_of_weights += 1.;
279  number_of_hits += 1.;
280  }
281  }
282  return sum_of_weights;
283 }
MuonResidualsFitter::residuals_begin
std::vector< double * >::const_iterator residuals_begin() const
Definition: MuonResidualsFitter.h:190
MuonResidualsFitter::residuals_end
std::vector< double * >::const_iterator residuals_end() const
Definition: MuonResidualsFitter.h:191
int MuonResiduals6DOFFitter::type ( ) const
inlineoverridevirtual

Implements MuonResidualsFitter.

Definition at line 65 of file MuonResiduals6DOFFitter.h.

References MuonResidualsFitter::k6DOF.

65 { return MuonResidualsFitter::k6DOF; }
Generated for CMSSW Reference Manual by doxygen 1.8.5