MuonResidualsTwoBin.h

Go to the documentation of this file.

#ifndef Alignment_MuonAlignmentAlgorithms_MuonResidualsTwoBin_H
#define Alignment_MuonAlignmentAlgorithms_MuonResidualsTwoBin_H

#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "Alignment/MuonAlignmentAlgorithms/interface/MuonResidualsFitter.h"
#include "Alignment/CommonAlignment/interface/Alignable.h"
#include "TMath.h"

class MuonResidualsTwoBin {
public:
  MuonResidualsTwoBin(bool twoBin, MuonResidualsFitter *pos, MuonResidualsFitter *neg)
      : m_twoBin(twoBin), m_pos(pos), m_neg(neg){};
  ~MuonResidualsTwoBin() {
    if (m_pos != nullptr)
      delete m_pos;
    if (m_neg != nullptr)
      delete m_neg;
  };

  int residualsModel() const {
    assert(m_pos->residualsModel() == m_neg->residualsModel());
    return m_pos->residualsModel();
  };
  long numResidualsPos() const { return m_pos->numResiduals(); };
  long numResidualsNeg() const { return m_neg->numResiduals(); };
  int npar() {
    assert(m_pos->npar() == m_neg->npar());
    return m_pos->npar();
  };
  int ndata() {
    assert(m_pos->ndata() == m_neg->ndata());
    return m_pos->ndata();
  };
  int type() const {
    assert(m_pos->type() == m_neg->type());
    return m_pos->type();
  };
  int useRes() const { return m_pos->useRes(); };

  void fix(int parNum, bool value = true) {
    m_pos->fix(parNum, value);
    m_neg->fix(parNum, value);
  };

  bool fixed(int parNum) { return m_pos->fixed(parNum) && m_neg->fixed(parNum); };

  void setPrintLevel(int printLevel) const {
    m_pos->setPrintLevel(printLevel);
    m_neg->setPrintLevel(printLevel);
  }

  void setStrategy(int strategy) const {
    m_pos->setStrategy(strategy);
    m_neg->setStrategy(strategy);
  }

  void fill(char charge, double *residual) {
    if (!m_twoBin || charge > 0)
      m_pos->fill(residual);
    else
      m_neg->fill(residual);
  };

  bool fit(Alignable *ali) { return (m_twoBin ? (m_pos->fit(ali) && m_neg->fit(ali)) : m_pos->fit(ali)); };
  double value(int parNum) {
    return (m_twoBin ? ((m_pos->value(parNum) + m_neg->value(parNum)) / 2.) : m_pos->value(parNum));
  };
  double errorerror(int parNum) {
    return (m_twoBin ? (sqrt(pow(m_pos->errorerror(parNum), 2.) + pow(m_neg->errorerror(parNum), 2.)) / 2.)
                     : m_pos->errorerror(parNum));
  };
  double antisym(int parNum) { return (m_twoBin ? ((m_pos->value(parNum) - m_neg->value(parNum)) / 2.) : 0.); };
  double loglikelihood() {
    return (m_twoBin ? (m_pos->loglikelihood() + m_neg->loglikelihood()) : m_pos->loglikelihood());
  };
  double numsegments() { return (m_twoBin ? (m_pos->numsegments() + m_neg->numsegments()) : m_pos->numsegments()); };
  double sumofweights() {
    return (m_twoBin ? (m_pos->sumofweights() + m_neg->sumofweights()) : m_pos->sumofweights());
  };

  // demonstration plots
  double plot(std::string name, TFileDirectory *dir, Alignable *ali) {
    if (m_twoBin) {
      std::string namePos = name + std::string("Pos");
      std::string nameNeg = name + std::string("Neg");
      double output = 0.;
      output += m_pos->plot(namePos, dir, ali);
      output += m_neg->plot(nameNeg, dir, ali);
      return output;
    } else {
      return m_pos->plot(name, dir, ali);
    }
  };

  // I/O of temporary files for collect mode
  void write(FILE *file, int which = 0) {
    if (m_twoBin) {
      m_pos->write(file, 2 * which);
      m_neg->write(file, 2 * which + 1);
    } else {
      m_pos->write(file, which);
    }
  };
  void read(FILE *file, int which = 0) {
    if (m_twoBin) {
      m_pos->read(file, 2 * which);
      m_neg->read(file, 2 * which + 1);
    } else {
      m_pos->read(file, which);
    }
  };

  double median(int which) {
    std::vector<double> residuals;
    for (std::vector<double *>::const_iterator r = residualsPos_begin(); r != residualsPos_end(); ++r) {
      residuals.push_back((*r)[which]);
    }
    if (m_twoBin) {
      for (std::vector<double *>::const_iterator r = residualsNeg_begin(); r != residualsNeg_end(); ++r) {
        residuals.push_back((*r)[which]);
      }
    }
    std::sort(residuals.begin(), residuals.end());
    int length = residuals.size();
    return residuals[length / 2];
  };

  double mean(int which, double truncate) {
    double sum = 0.;
    double n = 0.;
    for (std::vector<double *>::const_iterator r = residualsPos_begin(); r != residualsPos_end(); ++r) {
      double value = (*r)[which];
      if (fabs(value) < truncate) {
        sum += value;
        n += 1.;
      }
    }
    if (m_twoBin) {
      for (std::vector<double *>::const_iterator r = residualsNeg_begin(); r != residualsNeg_end(); ++r) {
        double value = (*r)[which];
        if (fabs(value) < truncate) {
          sum += value;
          n += 1.;
        }
      }
    }
    return sum / n;
  };

  double wmean(int which, int whichredchi2, double truncate) {
    double sum = 0.;
    double n = 0.;
    for (std::vector<double *>::const_iterator r = residualsPos_begin(); r != residualsPos_end(); ++r) {
      double value = (*r)[which];
      if (fabs(value) < truncate) {
        double weight = 1. / (*r)[whichredchi2];
        if (TMath::Prob(1. / weight * 12, 12) < 0.99) {
          sum += weight * value;
          n += weight;
        }
      }
    }
    if (m_twoBin) {
      for (std::vector<double *>::const_iterator r = residualsNeg_begin(); r != residualsNeg_end(); ++r) {
        double value = (*r)[which];
        if (fabs(value) < truncate) {
          double weight = 1. / (*r)[whichredchi2];
          if (TMath::Prob(1. / weight * 12, 12) < 0.99) {
            sum += weight * value;
            n += weight;
          }
        }
      }
    }
    return sum / n;
  };

  double stdev(int which, double truncate) {
    double sum2 = 0.;
    double sum = 0.;
    double n = 0.;
    for (std::vector<double *>::const_iterator r = residualsPos_begin(); r != residualsPos_end(); ++r) {
      double value = (*r)[which];
      if (fabs(value) < truncate) {
        sum2 += value * value;
        sum += value;
        n += 1.;
      }
    }
    if (m_twoBin) {
      for (std::vector<double *>::const_iterator r = residualsNeg_begin(); r != residualsNeg_end(); ++r) {
        double value = (*r)[which];
        if (fabs(value) < truncate) {
          sum2 += value * value;
          sum += value;
          n += 1.;
        }
      }
    }
    return sqrt(sum2 / n - pow(sum / n, 2));
  };

  void plotsimple(std::string name, TFileDirectory *dir, int which, double multiplier) {
    if (m_twoBin) {
      std::string namePos = name + std::string("Pos");
      std::string nameNeg = name + std::string("Neg");
      m_pos->plotsimple(namePos, dir, which, multiplier);
      m_neg->plotsimple(nameNeg, dir, which, multiplier);
    } else {
      m_pos->plotsimple(name, dir, which, multiplier);
    }
  };

  void plotweighted(std::string name, TFileDirectory *dir, int which, int whichredchi2, double multiplier) {
    if (m_twoBin) {
      std::string namePos = name + std::string("Pos");
      std::string nameNeg = name + std::string("Neg");
      m_pos->plotweighted(namePos, dir, which, whichredchi2, multiplier);
      m_neg->plotweighted(nameNeg, dir, which, whichredchi2, multiplier);
    } else {
      m_pos->plotweighted(name, dir, which, whichredchi2, multiplier);
    }
  };

  void selectPeakResiduals(double nsigma, int nvar, int *vars) {
    if (m_twoBin) {
      m_pos->selectPeakResiduals(nsigma, nvar, vars);
      m_neg->selectPeakResiduals(nsigma, nvar, vars);
    } else {
      m_pos->selectPeakResiduals(nsigma, nvar, vars);
    }
  }

  void correctBField() {
    m_pos->correctBField();
    //if (m_twoBin) m_neg->correctBField();
  };

  void fiducialCuts() { m_pos->fiducialCuts(); };

  void eraseNotSelectedResiduals() {
    if (m_twoBin) {
      m_pos->eraseNotSelectedResiduals();
      m_neg->eraseNotSelectedResiduals();
    } else {
      m_pos->eraseNotSelectedResiduals();
    }
  }

  std::vector<double *>::const_iterator residualsPos_begin() const { return m_pos->residuals_begin(); };
  std::vector<double *>::const_iterator residualsPos_end() const { return m_pos->residuals_end(); };
  std::vector<double *>::const_iterator residualsNeg_begin() const { return m_neg->residuals_begin(); };
  std::vector<double *>::const_iterator residualsNeg_end() const { return m_neg->residuals_end(); };

  std::vector<bool>::const_iterator residualsPos_ok_begin() const { return m_pos->selectedResidualsFlags().begin(); };
  std::vector<bool>::const_iterator residualsPos_ok_end() const { return m_pos->selectedResidualsFlags().end(); };
  std::vector<bool>::const_iterator residualsNeg_ok_begin() const { return m_neg->selectedResidualsFlags().begin(); };
  std::vector<bool>::const_iterator residualsNeg_ok_end() const { return m_neg->selectedResidualsFlags().end(); };

protected:
  bool m_twoBin;
  MuonResidualsFitter *m_pos, *m_neg;
};

#endif  // Alignment_MuonAlignmentAlgorithms_MuonResidualsTwoBin_H