MuonResiduals6DOFrphiFitter Class Reference

#include <MuonResiduals6DOFrphiFitter.h>

Inheritance diagram for MuonResiduals6DOFrphiFitter:

Public Types
enum	{   kAlignX = 0, kAlignY, kAlignZ, kAlignPhiX,   kAlignPhiY, kAlignPhiZ, kResidSigma, kResSlopeSigma,   kAlpha, kResidGamma, kResSlopeGamma, kNPar }
enum	{   kResid = 0, kResSlope, kPositionX, kPositionY,   kAngleX, kAngleY, kRedChi2, kPz,   kPt, kCharge, 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
	MuonResiduals6DOFrphiFitter (int residualsModel, int minHits, int useResiduals, bool weightAlignment=true)
	MuonResiduals6DOFrphiFitter (int residualsModel, int minHits, int useResiduals, const CSCGeometry *cscGeometry, 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 endcap, unsigned int station, unsigned int ring, unsigned int chamber, unsigned int preselected=1)
double	sumofweights () override
int	type () const override
	~MuonResiduals6DOFrphiFitter () 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 21:29:15 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 19 of file MuonResiduals6DOFrphiFitter.h.

Member Enumeration Documentation

anonymous enum

anonymous enum

Enumerator
kAlignX
kAlignY
kAlignZ
kAlignPhiX
kAlignPhiY
kAlignPhiZ
kResidSigma
kResSlopeSigma
kAlpha
kResidGamma
kResSlopeGamma
kNPar

Definition at line 21 of file MuonResiduals6DOFrphiFitter.h.

       {
    kAlignX = 0,
    kAlignY,
    kAlignZ,
    kAlignPhiX,
    kAlignPhiY,
    kAlignPhiZ,
    kResidSigma,
    kResSlopeSigma,
    kAlpha,
    kResidGamma,
    kResSlopeGamma,
    kNPar
  };

anonymous enum

anonymous enum

Enumerator
kResid
kResSlope
kPositionX
kPositionY
kAngleX
kAngleY
kRedChi2
kPz
kPt
kCharge
kNData

Definition at line 36 of file MuonResiduals6DOFrphiFitter.h.

{ kResid = 0, kResSlope, kPositionX, kPositionY, kAngleX, kAngleY, kRedChi2, kPz, kPt, kCharge, kNData };

Constructor & Destructor Documentation

MuonResiduals6DOFrphiFitter() [1/2]

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

inline

Definition at line 38 of file MuonResiduals6DOFrphiFitter.h.

      : MuonResidualsFitter(residualsModel, minHits, useResiduals, weightAlignment) {}

MuonResiduals6DOFrphiFitter() [2/2]

MuonResiduals6DOFrphiFitter::MuonResiduals6DOFrphiFitter ( int  residualsModel, int  minHits, int  useResiduals, const CSCGeometry *  cscGeometry, bool  weightAlignment = true  )

inline

Definition at line 42 of file MuonResiduals6DOFrphiFitter.h.

      : MuonResidualsFitter(residualsModel, minHits, useResiduals, weightAlignment) {}

~MuonResiduals6DOFrphiFitter()

MuonResiduals6DOFrphiFitter::~MuonResiduals6DOFrphiFitter ( )

inlineoverride

Definition at line 47 of file MuonResiduals6DOFrphiFitter.h.

{}

Member Function Documentation

correctBField()

void MuonResiduals6DOFrphiFitter::correctBField ( )

overridevirtual

Implements MuonResidualsFitter.

Definition at line 192 of file MuonResiduals6DOFrphiFitter.cc.

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

{ MuonResidualsFitter::correctBField(kPz, kCharge); }

fit()

bool MuonResiduals6DOFrphiFitter::fit ( Alignable *  ali )

overridevirtual

Implements MuonResidualsFitter.

Definition at line 215 of file MuonResiduals6DOFrphiFitter.cc.

References MuonResidualsFitter::dofit(), MuonResidualsFitter::fix(), Alignable::globalPosition(), LaserClient_cfi::high, mps_fire::i, heavyIonCSV_trainingSettings::idx, MuonResidualsFitter::initialize_table(), MuonResidualsFitter::k0010, MuonResidualsFitter::k1010, MuonResidualsFitter::k1100, MuonResidualsFitter::k1110, MuonResidualsFitter::k1111, kAlignPhiX, kAlignPhiY, kAlignPhiZ, kAlignX, kAlignZ, kAlpha, MuonResidualsFitter::kPowerLawTails, MuonResidualsFitter::kPureGaussian2D, kResidGamma, kResidSigma, kResSlopeGamma, kResSlopeSigma, MuonResidualsFitter::kROOTVoigt, LaserClient_cfi::low, MuonResiduals6DOFrphiFitter_FCN(), Skims_PA_cff::name, names, EgammaValidation_cff::num, funct::pow(), MuonResidualsFitter::residualsModel(), mathSSE::sqrt(), command_line::start, customisers::steps, AlCaHLTBitMon_QueryRunRegistry::string, sumofweights(), MuonResidualsFitter::useRes(), PV3DBase< T, PVType, FrameType >::x(), and PV3DBase< T, PVType, FrameType >::y().

Referenced by trackingPlots.Iteration::modules().

                                                    {
  //cscGeometry = m_cscGeometry;
  //cscDetId = CSCDetId(ali->geomDetId().rawId());

#ifndef STANDALONE_FITTER
  csc_R = sqrt(pow(ali->globalPosition().x(), 2) + pow(ali->globalPosition().y(), 2));
#else
  csc_R = 200;  // not important what number it is, as we usually not use the DOF where it matters
#endif

  initialize_table();  // if not already initialized
  sumofweights();

  double res_std = 0.5;
  double resslope_std = 0.002;

  int nums[10] = {kAlignX,
                  kAlignZ,
                  kAlignPhiX,
                  kAlignPhiY,
                  kAlignPhiZ,
                  kResidSigma,
                  kResSlopeSigma,
                  kAlpha,
                  kResidGamma,
                  kResSlopeGamma};
  std::string names[10] = {"AlignX",
                           "AlignZ",
                           "AlignPhiX",
                           "AlignPhiY",
                           "AlignPhiZ",
                           "ResidSigma",
                           "ResSlopeSigma",
                           "Alpha",
                           "ResidGamma",
                           "ResSlopeGamma"};
  double starts[10] = {0., 0., 0., 0., 0., res_std, resslope_std, 0., 0.1 * res_std, 0.1 * resslope_std};
  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};
  double lows[10] = {0., 0., 0., 0., 0., 0., 0., -1., 0., 0.};
  double highs[10] = {0., 0., 0., 0., 0., 10., 0.1, 1., 0., 0.};

  std::vector<int> num(nums, nums + 5);
  std::vector<std::string> name(names, names + 5);
  std::vector<double> start(starts, starts + 5);
  std::vector<double> step(steps, steps + 5);
  std::vector<double> low(lows, lows + 5);
  std::vector<double> high(highs, highs + 5);

  bool add_alpha = (residualsModel() == kPureGaussian2D);
  bool add_gamma = (residualsModel() == kROOTVoigt || residualsModel() == kPowerLawTails);

  int idx[4], ni = 0;
  if (useRes() == k1111 || useRes() == k1110 || useRes() == k1010) {
    for (ni = 0; ni < 2; ni++)
      idx[ni] = ni + 5;
    if (add_alpha)
      idx[ni++] = 7;
    else if (add_gamma)
      for (; ni < 4; ni++)
        idx[ni] = ni + 6;
    if (!add_alpha)
      fix(kAlpha);
  } else if (useRes() == k1100) {
    idx[ni++] = 5;
    if (add_gamma)
      idx[ni++] = 8;
    fix(kResSlopeSigma);
    fix(kAlpha);
  } else if (useRes() == k0010) {
    idx[ni++] = 6;
    if (add_gamma)
      idx[ni++] = 9;
    fix(kResidSigma);
    fix(kAlpha);
  }
  for (int i = 0; i < ni; i++) {
    num.push_back(nums[idx[i]]);
    name.push_back(names[idx[i]]);
    start.push_back(starts[idx[i]]);
    step.push_back(steps[idx[i]]);
    low.push_back(lows[idx[i]]);
    high.push_back(highs[idx[i]]);
  }

  return dofit(&MuonResiduals6DOFrphiFitter_FCN, num, name, start, step, low, high);
}

inform()

void MuonResiduals6DOFrphiFitter::inform ( TMinuit *  tMinuit )

overrideprotectedvirtual

Implements MuonResidualsFitter.

Definition at line 194 of file MuonResiduals6DOFrphiFitter.cc.

{ minuit = tMinuit; }

ndata()

int MuonResiduals6DOFrphiFitter::ndata ( )

inlineoverridevirtual

Implements MuonResidualsFitter.

Definition at line 62 of file MuonResiduals6DOFrphiFitter.h.

References kNData.

{ return kNData; }

npar()

int MuonResiduals6DOFrphiFitter::npar ( )

inlineoverridevirtual

Implements MuonResidualsFitter.

Definition at line 51 of file MuonResiduals6DOFrphiFitter.h.

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

Referenced by plot().

                      {
    if (residualsModel() == kPureGaussian || residualsModel() == kPureGaussian2D ||
        residualsModel() == kGaussPowerTails)
      return kNPar - 2;
    else if (residualsModel() == kPowerLawTails)
      return kNPar;
    else if (residualsModel() == kROOTVoigt)
      return kNPar;
    else
      assert(false);
  }

plot()

double MuonResiduals6DOFrphiFitter::plot ( std::string  name, TFileDirectory *  dir, Alignable *  ali  )

overridevirtual

Implements MuonResidualsFitter.

Definition at line 303 of file MuonResiduals6DOFrphiFitter.cc.

References a, cms::cuda::assert(), b, hltPixelTracks_cff::chi2, DeadROC_duringRun::dir, MuonResidualsFitter::errorerror(), mps_fire::i, kAlignPhiX, kAlignPhiY, kAlignPhiZ, kAlignX, kAlignY, kAlignZ, kAlpha, kAngleX, kAngleY, MuonResidualsFitter::kGaussPowerTails, kPositionX, kPositionY, MuonResidualsFitter::kPowerLawTails, MuonResidualsFitter::kPureGaussian, MuonResidualsFitter::kPureGaussian2D, kRedChi2, kResid, kResidGamma, kResidSigma, kResSlope, kResSlopeGamma, kResSlopeSigma, MuonResidualsFitter::kROOTVoigt, AlignableSurface::length(), MuonResidualsFitter::m_weightAlignment, MuonResidualsFitter_GaussPowerTails_TF1(), MuonResidualsFitter_powerLawTails_TF1(), MuonResidualsFitter_pureGaussian_TF1(), MuonResidualsFitter_ROOTVoigt_TF1(), Skims_PA_cff::name, ndof, npar(), funct::pow(), alignCSCRings::r, MuonResidualsFitter::residuals_begin(), MuonResidualsFitter::residuals_end(), MuonResidualsFitter::residualsModel(), AlCaHLTBitMon_QueryRunRegistry::string, sumofweights(), Alignable::surface(), fftjetcommon_cfi::sx, fftjetcommon_cfi::sy, MuonResidualsFitter::value(), mps_merge::weight, ApeEstimator_cff::width, AlignableSurface::width(), and protons_cff::xi.

                                                                                            {
  sumofweights();

  double mean_residual = 0., mean_resslope = 0.;
  double mean_trackx = 0., mean_tracky = 0., mean_trackdxdz = 0., mean_trackdydz = 0.;
  double sum_w = 0.;

  for (std::vector<double *>::const_iterator rit = residuals_begin(); rit != residuals_end(); ++rit) {
    const double redchi2 = (*rit)[kRedChi2];
    double weight = 1. / redchi2;
    if (!m_weightAlignment)
      weight = 1.;

    if (!m_weightAlignment || TMath::Prob(redchi2 * 6, 6) < 0.99)  // no spikes allowed
    {
      double factor_w = 1. / (sum_w + weight);
      mean_residual = factor_w * (sum_w * mean_residual + weight * (*rit)[kResid]);
      mean_resslope = factor_w * (sum_w * mean_resslope + weight * (*rit)[kResSlope]);
      mean_trackx = factor_w * (sum_w * mean_trackx + weight * (*rit)[kPositionX]);
      mean_tracky = factor_w * (sum_w * mean_tracky + weight * (*rit)[kPositionY]);
      mean_trackdxdz = factor_w * (sum_w * mean_trackdxdz + weight * (*rit)[kAngleX]);
      mean_trackdydz = factor_w * (sum_w * mean_trackdydz + weight * (*rit)[kAngleY]);
      sum_w += weight;
    }
  }

  std::string name_residual, name_resslope, name_residual_raw, name_resslope_raw, name_residual_cut, name_alpha;
  std::string name_residual_trackx, name_resslope_trackx;
  std::string name_residual_tracky, name_resslope_tracky;
  std::string name_residual_trackdxdz, name_resslope_trackdxdz;
  std::string name_residual_trackdydz, name_resslope_trackdydz;

  name_residual = name + "_residual";
  name_resslope = name + "_resslope";
  name_residual_raw = name + "_residual_raw";
  name_resslope_raw = name + "_resslope_raw";
  name_residual_cut = name + "_residual_cut";
  name_alpha = name + "_alpha";
  name_residual_trackx = name + "_residual_trackx";
  name_resslope_trackx = name + "_resslope_trackx";
  name_residual_tracky = name + "_residual_tracky";
  name_resslope_tracky = name + "_resslope_tracky";
  name_residual_trackdxdz = name + "_residual_trackdxdz";
  name_resslope_trackdxdz = name + "_resslope_trackdxdz";
  name_residual_trackdydz = name + "_residual_trackdydz";
  name_resslope_trackdydz = name + "_resslope_trackdydz";

  double width = ali->surface().width();
  double length = ali->surface().length();
  int bins_residual = 100, bins_resslope = 100;
  double min_residual = -50., max_residual = 50.;
  double min_resslope = -50., max_resslope = 50.;
  double min_trackx = -width / 2., max_trackx = width / 2.;
  double min_tracky = -length / 2., max_tracky = length / 2.;
  double min_trackdxdz = -1.5, max_trackdxdz = 1.5;
  double min_trackdydz = -1.5, max_trackdydz = 1.5;

  TH1F *hist_residual = dir->make<TH1F>(name_residual.c_str(), "", bins_residual, min_residual, max_residual);
  TH1F *hist_resslope = dir->make<TH1F>(name_resslope.c_str(), "", bins_resslope, min_resslope, max_resslope);
  TH1F *hist_residual_raw = dir->make<TH1F>(name_residual_raw.c_str(), "", bins_residual, min_residual, max_residual);
  TH1F *hist_resslope_raw = dir->make<TH1F>(name_resslope_raw.c_str(), "", bins_resslope, min_resslope, max_resslope);
  TH1F *hist_residual_cut = dir->make<TH1F>(name_residual_cut.c_str(), "", bins_residual, min_residual, max_residual);
  TH2F *hist_alpha = dir->make<TH2F>(name_alpha.c_str(), "", 50, min_resslope, max_resslope, 50, -50., 50.);

  TProfile *hist_residual_trackx =
      dir->make<TProfile>(name_residual_trackx.c_str(), "", 50, min_trackx, max_trackx, min_residual, max_residual);
  TProfile *hist_resslope_trackx =
      dir->make<TProfile>(name_resslope_trackx.c_str(), "", 50, min_trackx, max_trackx, min_resslope, max_resslope);
  TProfile *hist_residual_tracky =
      dir->make<TProfile>(name_residual_tracky.c_str(), "", 50, min_tracky, max_tracky, min_residual, max_residual);
  TProfile *hist_resslope_tracky =
      dir->make<TProfile>(name_resslope_tracky.c_str(), "", 50, min_tracky, max_tracky, min_resslope, max_resslope);
  TProfile *hist_residual_trackdxdz = dir->make<TProfile>(
      name_residual_trackdxdz.c_str(), "", 500, min_trackdxdz, max_trackdxdz, min_residual, max_residual);
  TProfile *hist_resslope_trackdxdz = dir->make<TProfile>(
      name_resslope_trackdxdz.c_str(), "", 500, min_trackdxdz, max_trackdxdz, min_resslope, max_resslope);
  TProfile *hist_residual_trackdydz = dir->make<TProfile>(
      name_residual_trackdydz.c_str(), "", 500, min_trackdydz, max_trackdydz, min_residual, max_residual);
  TProfile *hist_resslope_trackdydz = dir->make<TProfile>(
      name_resslope_trackdydz.c_str(), "", 500, min_trackdydz, max_trackdydz, min_resslope, max_resslope);

  hist_residual_trackx->SetAxisRange(-10., 10., "Y");
  hist_resslope_trackx->SetAxisRange(-10., 10., "Y");
  hist_residual_tracky->SetAxisRange(-10., 10., "Y");
  hist_resslope_tracky->SetAxisRange(-10., 10., "Y");
  hist_residual_trackdxdz->SetAxisRange(-10., 10., "Y");
  hist_resslope_trackdxdz->SetAxisRange(-10., 10., "Y");
  hist_residual_trackdydz->SetAxisRange(-10., 10., "Y");
  hist_resslope_trackdydz->SetAxisRange(-10., 10., "Y");

  name_residual += "_fit";
  name_resslope += "_fit";
  name_alpha += "_fit";
  name_residual_trackx += "_fit";
  name_resslope_trackx += "_fit";
  name_residual_tracky += "_fit";
  name_resslope_tracky += "_fit";
  name_residual_trackdxdz += "_fit";
  name_resslope_trackdxdz += "_fit";
  name_residual_trackdydz += "_fit";
  name_resslope_trackdydz += "_fit";

  TF1 *fit_residual = nullptr;
  TF1 *fit_resslope = nullptr;
  if (residualsModel() == kPureGaussian || residualsModel() == kPureGaussian2D) {
    fit_residual = new TF1(name_residual.c_str(), MuonResidualsFitter_pureGaussian_TF1, min_residual, max_residual, 3);
    fit_residual->SetParameters(
        sum_of_weights * (max_residual - min_residual) / bins_residual, 10. * value(kAlignX), 10. * value(kResidSigma));
    const double er_res[3] = {0., 10. * errorerror(kAlignX), 10. * errorerror(kResidSigma)};
    fit_residual->SetParErrors(er_res);
    fit_resslope = new TF1(name_resslope.c_str(), MuonResidualsFitter_pureGaussian_TF1, min_resslope, max_resslope, 3);
    fit_resslope->SetParameters(sum_of_weights * (max_resslope - min_resslope) / bins_resslope,
                                1000. * value(kAlignPhiY),
                                1000. * value(kResSlopeSigma));
    const double er_resslope[3] = {0., 1000. * errorerror(kAlignPhiY), 1000. * errorerror(kResSlopeSigma)};
    fit_resslope->SetParErrors(er_resslope);
  } else if (residualsModel() == kPowerLawTails) {
    fit_residual = new TF1(name_residual.c_str(), MuonResidualsFitter_powerLawTails_TF1, min_residual, max_residual, 4);
    fit_residual->SetParameters(sum_of_weights * (max_residual - min_residual) / bins_residual,
                                10. * value(kAlignX),
                                10. * value(kResidSigma),
                                10. * value(kResidGamma));
    fit_resslope = new TF1(name_resslope.c_str(), MuonResidualsFitter_powerLawTails_TF1, min_resslope, max_resslope, 4);
    fit_resslope->SetParameters(sum_of_weights * (max_resslope - min_resslope) / bins_resslope,
                                1000. * value(kAlignPhiY),
                                1000. * value(kResSlopeSigma),
                                1000. * value(kResSlopeGamma));
  } else if (residualsModel() == kROOTVoigt) {
    fit_residual = new TF1(name_residual.c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_residual, max_residual, 4);
    fit_residual->SetParameters(sum_of_weights * (max_residual - min_residual) / bins_residual,
                                10. * value(kAlignX),
                                10. * value(kResidSigma),
                                10. * value(kResidGamma));
    fit_resslope = new TF1(name_resslope.c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_resslope, max_resslope, 4);
    fit_resslope->SetParameters(sum_of_weights * (max_resslope - min_resslope) / bins_resslope,
                                1000. * value(kAlignPhiY),
                                1000. * value(kResSlopeSigma),
                                1000. * value(kResSlopeGamma));
  } else if (residualsModel() == kGaussPowerTails) {
    fit_residual =
        new TF1(name_residual.c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_residual, max_residual, 3);
    fit_residual->SetParameters(
        sum_of_weights * (max_residual - min_residual) / bins_residual, 10. * value(kAlignX), 10. * value(kResidSigma));
    fit_resslope =
        new TF1(name_resslope.c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_resslope, max_resslope, 3);
    fit_resslope->SetParameters(sum_of_weights * (max_resslope - min_resslope) / bins_resslope,
                                1000. * value(kAlignPhiY),
                                1000. * value(kResSlopeSigma));
  } else {
    assert(false);
  }

  fit_residual->SetLineColor(2);
  fit_residual->SetLineWidth(2);
  fit_residual->Write();
  fit_resslope->SetLineColor(2);
  fit_resslope->SetLineWidth(2);
  fit_resslope->Write();

  TF1 *fit_alpha = new TF1(name_alpha.c_str(), "[0] + x*[1]", min_resslope, max_resslope);
  double a = 10. * value(kAlignX), b = 10. * value(kAlpha) / 1000.;
  if (residualsModel() == kPureGaussian2D) {
    double sx = 10. * value(kResidSigma), sy = 1000. * value(kResSlopeSigma), r = value(kAlpha);
    a = mean_residual;
    b = 0.;
    if (sx != 0.) {
      b = 1. / (sy / sx * r);
      a = -b * mean_resslope;
    }
  }
  fit_alpha->SetParameters(a, b);
  fit_alpha->SetLineColor(2);
  fit_alpha->SetLineWidth(2);
  fit_alpha->Write();

  TProfile *fit_residual_trackx = dir->make<TProfile>(name_residual_trackx.c_str(), "", 100, min_trackx, max_trackx);
  TProfile *fit_resslope_trackx = dir->make<TProfile>(name_resslope_trackx.c_str(), "", 100, min_trackx, max_trackx);
  TProfile *fit_residual_tracky = dir->make<TProfile>(name_residual_tracky.c_str(), "", 100, min_tracky, max_tracky);
  TProfile *fit_resslope_tracky = dir->make<TProfile>(name_resslope_tracky.c_str(), "", 100, min_tracky, max_tracky);
  TProfile *fit_residual_trackdxdz =
      dir->make<TProfile>(name_residual_trackdxdz.c_str(), "", 500, min_trackdxdz, max_trackdxdz);
  TProfile *fit_resslope_trackdxdz =
      dir->make<TProfile>(name_resslope_trackdxdz.c_str(), "", 500, min_trackdxdz, max_trackdxdz);
  TProfile *fit_residual_trackdydz =
      dir->make<TProfile>(name_residual_trackdydz.c_str(), "", 500, min_trackdydz, max_trackdydz);
  TProfile *fit_resslope_trackdydz =
      dir->make<TProfile>(name_resslope_trackdydz.c_str(), "", 500, min_trackdydz, max_trackdydz);

  fit_residual_trackx->SetLineColor(2);
  fit_residual_trackx->SetLineWidth(2);
  fit_resslope_trackx->SetLineColor(2);
  fit_resslope_trackx->SetLineWidth(2);
  fit_residual_tracky->SetLineColor(2);
  fit_residual_tracky->SetLineWidth(2);
  fit_resslope_tracky->SetLineColor(2);
  fit_resslope_tracky->SetLineWidth(2);
  fit_residual_trackdxdz->SetLineColor(2);
  fit_residual_trackdxdz->SetLineWidth(2);
  fit_resslope_trackdxdz->SetLineColor(2);
  fit_resslope_trackdxdz->SetLineWidth(2);
  fit_residual_trackdydz->SetLineColor(2);
  fit_residual_trackdydz->SetLineWidth(2);
  fit_resslope_trackdydz->SetLineColor(2);
  fit_resslope_trackdydz->SetLineWidth(2);

  name_residual_trackx += "line";
  name_resslope_trackx += "line";
  name_residual_tracky += "line";
  name_resslope_tracky += "line";
  name_residual_trackdxdz += "line";
  name_resslope_trackdxdz += "line";
  name_residual_trackdydz += "line";
  name_resslope_trackdydz += "line";

  TF1 *fitline_residual_trackx = new TF1(name_residual_trackx.c_str(), residual_trackx_TF1, min_trackx, max_trackx, 12);
  TF1 *fitline_resslope_trackx = new TF1(name_resslope_trackx.c_str(), resslope_trackx_TF1, min_trackx, max_trackx, 12);
  TF1 *fitline_residual_tracky = new TF1(name_residual_tracky.c_str(), residual_tracky_TF1, min_tracky, max_tracky, 12);
  TF1 *fitline_resslope_tracky = new TF1(name_resslope_tracky.c_str(), resslope_tracky_TF1, min_tracky, max_tracky, 12);
  TF1 *fitline_residual_trackdxdz =
      new TF1(name_residual_trackdxdz.c_str(), residual_trackdxdz_TF1, min_trackdxdz, max_trackdxdz, 12);
  TF1 *fitline_resslope_trackdxdz =
      new TF1(name_resslope_trackdxdz.c_str(), resslope_trackdxdz_TF1, min_trackdxdz, max_trackdxdz, 12);
  TF1 *fitline_residual_trackdydz =
      new TF1(name_residual_trackdydz.c_str(), residual_trackdydz_TF1, min_trackdydz, max_trackdydz, 12);
  TF1 *fitline_resslope_trackdydz =
      new TF1(name_resslope_trackdydz.c_str(), resslope_trackdydz_TF1, min_trackdydz, max_trackdydz, 12);

  std::vector<TF1 *> fitlines;
  fitlines.push_back(fitline_residual_trackx);
  fitlines.push_back(fitline_resslope_trackx);
  fitlines.push_back(fitline_residual_tracky);
  fitlines.push_back(fitline_resslope_tracky);
  fitlines.push_back(fitline_residual_trackdxdz);
  fitlines.push_back(fitline_resslope_trackdxdz);
  fitlines.push_back(fitline_residual_trackdydz);
  fitlines.push_back(fitline_resslope_trackdydz);

  double fitparameters[12] = {value(kAlignX),
                              value(kAlignY),
                              value(kAlignZ),
                              value(kAlignPhiX),
                              value(kAlignPhiY),
                              value(kAlignPhiZ),
                              mean_trackx,
                              mean_tracky,
                              mean_trackdxdz,
                              mean_trackdydz,
                              value(kAlpha),
                              mean_resslope};
  if (residualsModel() == kPureGaussian2D)
    fitparameters[10] = 0.;

  for (std::vector<TF1 *>::const_iterator itr = fitlines.begin(); itr != fitlines.end(); itr++) {
    (*itr)->SetParameters(fitparameters);
    (*itr)->SetLineColor(2);
    (*itr)->SetLineWidth(2);
    (*itr)->Write();
  }

  for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
    const double resid = (*resiter)[kResid];
    const double resslope = (*resiter)[kResSlope];
    const double positionX = (*resiter)[kPositionX];
    const double positionY = (*resiter)[kPositionY];
    const double angleX = (*resiter)[kAngleX];
    const double angleY = (*resiter)[kAngleY];
    const double redchi2 = (*resiter)[kRedChi2];
    double weight = 1. / redchi2;
    if (!m_weightAlignment)
      weight = 1.;

    if (!m_weightAlignment || TMath::Prob(redchi2 * 6, 6) < 0.99) {  // no spikes allowed
      hist_alpha->Fill(1000. * resslope, 10. * resid);

      double coeff = value(kAlpha);
      if (residualsModel() == kPureGaussian || residualsModel() == kPureGaussian2D)
        coeff = 0.;
      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);
      hist_residual->Fill(10. * (resid - geom_resid + value(kAlignX)), weight);
      hist_residual_trackx->Fill(positionX, 10. * resid, weight);
      hist_residual_tracky->Fill(positionY, 10. * resid, weight);
      hist_residual_trackdxdz->Fill(angleX, 10. * resid, weight);
      hist_residual_trackdydz->Fill(angleY, 10. * resid, weight);
      fit_residual_trackx->Fill(positionX, 10. * geom_resid, weight);
      fit_residual_tracky->Fill(positionY, 10. * geom_resid, weight);
      fit_residual_trackdxdz->Fill(angleX, 10. * geom_resid, weight);
      fit_residual_trackdydz->Fill(angleY, 10. * geom_resid, weight);

      double geom_resslope = getResSlope(value(kAlignX),
                                         value(kAlignY),
                                         value(kAlignZ),
                                         value(kAlignPhiX),
                                         value(kAlignPhiY),
                                         value(kAlignPhiZ),
                                         positionX,
                                         positionY,
                                         angleX,
                                         angleY,
                                         effectiveR(positionX, positionY));
      hist_resslope->Fill(1000. * (resslope - geom_resslope + value(kAlignPhiY)), weight);
      hist_resslope_trackx->Fill(positionX, 1000. * resslope, weight);
      hist_resslope_tracky->Fill(positionY, 1000. * resslope, weight);
      hist_resslope_trackdxdz->Fill(angleX, 1000. * resslope, weight);
      hist_resslope_trackdydz->Fill(angleY, 1000. * resslope, weight);
      fit_resslope_trackx->Fill(positionX, 1000. * geom_resslope, weight);
      fit_resslope_tracky->Fill(positionY, 1000. * geom_resslope, weight);
      fit_resslope_trackdxdz->Fill(angleX, 1000. * geom_resslope, weight);
      fit_resslope_trackdydz->Fill(angleY, 1000. * geom_resslope, weight);
    }

    hist_residual_raw->Fill(10. * resid);
    hist_resslope_raw->Fill(1000. * resslope);
    if (fabs(resslope) < 0.005)
      hist_residual_cut->Fill(10. * resid);
  }

  double chi2 = 0.;
  double ndof = 0.;
  for (int i = 1; i <= hist_residual->GetNbinsX(); i++) {
    double xi = hist_residual->GetBinCenter(i);
    double yi = hist_residual->GetBinContent(i);
    double yerri = hist_residual->GetBinError(i);
    double yth = fit_residual->Eval(xi);
    if (yerri > 0.) {
      chi2 += pow((yth - yi) / yerri, 2);
      ndof += 1.;
    }
  }
  for (int i = 1; i <= hist_resslope->GetNbinsX(); i++) {
    double xi = hist_resslope->GetBinCenter(i);
    double yi = hist_resslope->GetBinContent(i);
    double yerri = hist_resslope->GetBinError(i);
    double yth = fit_resslope->Eval(xi);
    if (yerri > 0.) {
      chi2 += pow((yth - yi) / yerri, 2);
      ndof += 1.;
    }
  }
  ndof -= npar();

  return (ndof > 0. ? chi2 / ndof : -1.);
}

readNtuple()

TTree * MuonResiduals6DOFrphiFitter::readNtuple	(	std::string	fname,
		unsigned int	endcap,
		unsigned int	station,
		unsigned int	ring,
		unsigned int	chamber,
		unsigned int	preselected = 1
	)

Definition at line 658 of file MuonResiduals6DOFrphiFitter.cc.

References cms::cuda::assert(), relativeConstraints::chamber, makeMuonMisalignmentScenario::endcap, f, MuonResidualsFitter::fill(), alignmentValidation::fname, mps_fire::i, kAngleX, kAngleY, kCharge, kNData, kPositionX, kPositionY, kPt, kPz, kRedChi2, kResid, kResSlope, alignCSCRings::r, relativeConstraints::ring, relativeConstraints::station, submitPVValidationJobs::t, and groupFilesInBlocks::tt.

                                                                         {
  TFile *f = new TFile(fname.c_str());
  TTree *t = (TTree *)f->Get("mual_ttree");

  // Create  new temporary file
  TFile *tmpf = new TFile("small_tree_fit.root", "recreate");
  assert(tmpf != nullptr);

  // filter the tree (temporarily lives in the new file)
  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));

  MuonAlignmentTreeRow r;
  tt->SetBranchAddress("res_x", &r.res_x);
  tt->SetBranchAddress("res_slope_x", &r.res_slope_x);
  tt->SetBranchAddress("pos_x", &r.pos_x);
  tt->SetBranchAddress("pos_y", &r.pos_y);
  tt->SetBranchAddress("angle_x", &r.angle_x);
  tt->SetBranchAddress("angle_y", &r.angle_y);
  tt->SetBranchAddress("pz", &r.pz);
  tt->SetBranchAddress("pt", &r.pt);
  tt->SetBranchAddress("q", &r.q);

  Long64_t nentries = tt->GetEntries();
  for (Long64_t i = 0; i < nentries; i++) {
    tt->GetEntry(i);
    double *rdata = new double[MuonResiduals6DOFrphiFitter::kNData];
    rdata[kResid] = r.res_x;
    rdata[kResSlope] = r.res_slope_x;
    rdata[kPositionX] = r.pos_x;
    rdata[kPositionY] = r.pos_y;
    rdata[kAngleX] = r.angle_x;
    rdata[kAngleY] = r.angle_y;
    rdata[kRedChi2] = 0.1;
    rdata[kPz] = r.pz;
    rdata[kPt] = r.pt;
    rdata[kCharge] = r.q;
    fill(rdata);
  }
  delete f;
  //delete tmpf;
  return tt;
}

sumofweights()

double MuonResiduals6DOFrphiFitter::sumofweights ( )

overridevirtual

Implements MuonResidualsFitter.

Definition at line 196 of file MuonResiduals6DOFrphiFitter.cc.

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

Referenced by fit(), and plot().

                                                 {
  sum_of_weights = 0.;
  number_of_hits = 0.;
  weight_alignment = m_weightAlignment;
  for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
    if (m_weightAlignment) {
      double redchi2 = (*resiter)[kRedChi2];
      if (TMath::Prob(redchi2 * 6, 6) < 0.99) {  // no spikes allowed
        sum_of_weights += 1. / redchi2;
        number_of_hits += 1.;
      }
    } else {
      sum_of_weights += 1.;
      number_of_hits += 1.;
    }
  }
  return sum_of_weights;
}

type()

int MuonResiduals6DOFrphiFitter::type ( ) const

inlineoverridevirtual

Implements MuonResidualsFitter.

Definition at line 49 of file MuonResiduals6DOFrphiFitter.h.

References MuonResidualsFitter::k6DOFrphi.

{ return MuonResidualsFitter::k6DOFrphi; }