MuonAlignmentAlgorithms/interface/MuonResidualsTwoBin.h

0001 #ifndef Alignment_MuonAlignmentAlgorithms_MuonResidualsTwoBin_H
0002 #define Alignment_MuonAlignmentAlgorithms_MuonResidualsTwoBin_H
0003
0004 /** \class MuonResidualsTwoBin
0005  *  $Date: 2011/04/15 21:51:13 $
0006  *  $Revision: 1.10 $
0007  *  \author J. Pivarski - Texas A&M University <pivarski@physics.tamu.edu>
0008  */
0009
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResidualsFitter.h"
0012 #include "Alignment/CommonAlignment/interface/Alignable.h"
0013 #include "TMath.h"
0014
0015 class MuonResidualsTwoBin {
0016 public:
0017   MuonResidualsTwoBin(bool twoBin, MuonResidualsFitter *pos, MuonResidualsFitter *neg)
0018       : m_twoBin(twoBin), m_pos(pos), m_neg(neg){};
0019   ~MuonResidualsTwoBin() {
0020     if (m_pos != nullptr)
0021       delete m_pos;
0022     if (m_neg != nullptr)
0023       delete m_neg;
0024   };
0025
0026   int residualsModel() const {
0027     assert(m_pos->residualsModel() == m_neg->residualsModel());
0028     return m_pos->residualsModel();
0029   };
0030   long numResidualsPos() const { return m_pos->numResiduals(); };
0031   long numResidualsNeg() const { return m_neg->numResiduals(); };
0032   int npar() {
0033     assert(m_pos->npar() == m_neg->npar());
0034     return m_pos->npar();
0035   };
0036   int ndata() {
0037     assert(m_pos->ndata() == m_neg->ndata());
0038     return m_pos->ndata();
0039   };
0040   int type() const {
0041     assert(m_pos->type() == m_neg->type());
0042     return m_pos->type();
0043   };
0044   int useRes() const { return m_pos->useRes(); };
0045
0046   void fix(int parNum, bool value = true) {
0047     m_pos->fix(parNum, value);
0048     m_neg->fix(parNum, value);
0049   };
0050
0051   bool fixed(int parNum) { return m_pos->fixed(parNum) && m_neg->fixed(parNum); };
0052
0053   void setPrintLevel(int printLevel) const {
0054     m_pos->setPrintLevel(printLevel);
0055     m_neg->setPrintLevel(printLevel);
0056   }
0057
0058   void setStrategy(int strategy) const {
0059     m_pos->setStrategy(strategy);
0060     m_neg->setStrategy(strategy);
0061   }
0062
0063   void fill(char charge, double *residual) {
0064     if (!m_twoBin || charge > 0)
0065       m_pos->fill(residual);
0066     else
0067       m_neg->fill(residual);
0068   };
0069
0070   bool fit(Alignable *ali) { return (m_twoBin ? (m_pos->fit(ali) && m_neg->fit(ali)) : m_pos->fit(ali)); };
0071   double value(int parNum) {
0072     return (m_twoBin ? ((m_pos->value(parNum) + m_neg->value(parNum)) / 2.) : m_pos->value(parNum));
0073   };
0074   double errorerror(int parNum) {
0075     return (m_twoBin ? (sqrt(pow(m_pos->errorerror(parNum), 2.) + pow(m_neg->errorerror(parNum), 2.)) / 2.)
0076                      : m_pos->errorerror(parNum));
0077   };
0078   double antisym(int parNum) { return (m_twoBin ? ((m_pos->value(parNum) - m_neg->value(parNum)) / 2.) : 0.); };
0079   double loglikelihood() {
0080     return (m_twoBin ? (m_pos->loglikelihood() + m_neg->loglikelihood()) : m_pos->loglikelihood());
0081   };
0082   double numsegments() { return (m_twoBin ? (m_pos->numsegments() + m_neg->numsegments()) : m_pos->numsegments()); };
0083   double sumofweights() {
0084     return (m_twoBin ? (m_pos->sumofweights() + m_neg->sumofweights()) : m_pos->sumofweights());
0085   };
0086
0087   // demonstration plots
0088   double plot(std::string name, TFileDirectory *dir, Alignable *ali) {
0089     if (m_twoBin) {
0090       std::string namePos = name + std::string("Pos");
0091       std::string nameNeg = name + std::string("Neg");
0092       double output = 0.;
0093       output += m_pos->plot(namePos, dir, ali);
0094       output += m_neg->plot(nameNeg, dir, ali);
0095       return output;
0096     } else {
0097       return m_pos->plot(name, dir, ali);
0098     }
0099   };
0100
0101   // I/O of temporary files for collect mode
0102   void write(FILE *file, int which = 0) {
0103     if (m_twoBin) {
0104       m_pos->write(file, 2 * which);
0105       m_neg->write(file, 2 * which + 1);
0106     } else {
0107       m_pos->write(file, which);
0108     }
0109   };
0110   void read(FILE *file, int which = 0) {
0111     if (m_twoBin) {
0112       m_pos->read(file, 2 * which);
0113       m_neg->read(file, 2 * which + 1);
0114     } else {
0115       m_pos->read(file, which);
0116     }
0117   };
0118
0119   double median(int which) {
0120     std::vector<double> residuals;
0121     for (std::vector<double *>::const_iterator r = residualsPos_begin(); r != residualsPos_end(); ++r) {
0122       residuals.push_back((*r)[which]);
0123     }
0124     if (m_twoBin) {
0125       for (std::vector<double *>::const_iterator r = residualsNeg_begin(); r != residualsNeg_end(); ++r) {
0126         residuals.push_back((*r)[which]);
0127       }
0128     }
0129     std::sort(residuals.begin(), residuals.end());
0130     int length = residuals.size();
0131     return residuals[length / 2];
0132   };
0133
0134   double mean(int which, double truncate) {
0135     double sum = 0.;
0136     double n = 0.;
0137     for (std::vector<double *>::const_iterator r = residualsPos_begin(); r != residualsPos_end(); ++r) {
0138       double value = (*r)[which];
0139       if (fabs(value) < truncate) {
0140         sum += value;
0141         n += 1.;
0142       }
0143     }
0144     if (m_twoBin) {
0145       for (std::vector<double *>::const_iterator r = residualsNeg_begin(); r != residualsNeg_end(); ++r) {
0146         double value = (*r)[which];
0147         if (fabs(value) < truncate) {
0148           sum += value;
0149           n += 1.;
0150         }
0151       }
0152     }
0153     return sum / n;
0154   };
0155
0156   double wmean(int which, int whichredchi2, double truncate) {
0157     double sum = 0.;
0158     double n = 0.;
0159     for (std::vector<double *>::const_iterator r = residualsPos_begin(); r != residualsPos_end(); ++r) {
0160       double value = (*r)[which];
0161       if (fabs(value) < truncate) {
0162         double weight = 1. / (*r)[whichredchi2];
0163         if (TMath::Prob(1. / weight * 12, 12) < 0.99) {
0164           sum += weight * value;
0165           n += weight;
0166         }
0167       }
0168     }
0169     if (m_twoBin) {
0170       for (std::vector<double *>::const_iterator r = residualsNeg_begin(); r != residualsNeg_end(); ++r) {
0171         double value = (*r)[which];
0172         if (fabs(value) < truncate) {
0173           double weight = 1. / (*r)[whichredchi2];
0174           if (TMath::Prob(1. / weight * 12, 12) < 0.99) {
0175             sum += weight * value;
0176             n += weight;
0177           }
0178         }
0179       }
0180     }
0181     return sum / n;
0182   };
0183
0184   double stdev(int which, double truncate) {
0185     double sum2 = 0.;
0186     double sum = 0.;
0187     double n = 0.;
0188     for (std::vector<double *>::const_iterator r = residualsPos_begin(); r != residualsPos_end(); ++r) {
0189       double value = (*r)[which];
0190       if (fabs(value) < truncate) {
0191         sum2 += value * value;
0192         sum += value;
0193         n += 1.;
0194       }
0195     }
0196     if (m_twoBin) {
0197       for (std::vector<double *>::const_iterator r = residualsNeg_begin(); r != residualsNeg_end(); ++r) {
0198         double value = (*r)[which];
0199         if (fabs(value) < truncate) {
0200           sum2 += value * value;
0201           sum += value;
0202           n += 1.;
0203         }
0204       }
0205     }
0206     return sqrt(sum2 / n - pow(sum / n, 2));
0207   };
0208
0209   void plotsimple(std::string name, TFileDirectory *dir, int which, double multiplier) {
0210     if (m_twoBin) {
0211       std::string namePos = name + std::string("Pos");
0212       std::string nameNeg = name + std::string("Neg");
0213       m_pos->plotsimple(namePos, dir, which, multiplier);
0214       m_neg->plotsimple(nameNeg, dir, which, multiplier);
0215     } else {
0216       m_pos->plotsimple(name, dir, which, multiplier);
0217     }
0218   };
0219
0220   void plotweighted(std::string name, TFileDirectory *dir, int which, int whichredchi2, double multiplier) {
0221     if (m_twoBin) {
0222       std::string namePos = name + std::string("Pos");
0223       std::string nameNeg = name + std::string("Neg");
0224       m_pos->plotweighted(namePos, dir, which, whichredchi2, multiplier);
0225       m_neg->plotweighted(nameNeg, dir, which, whichredchi2, multiplier);
0226     } else {
0227       m_pos->plotweighted(name, dir, which, whichredchi2, multiplier);
0228     }
0229   };
0230
0231   void selectPeakResiduals(double nsigma, int nvar, int *vars) {
0232     if (m_twoBin) {
0233       m_pos->selectPeakResiduals(nsigma, nvar, vars);
0234       m_neg->selectPeakResiduals(nsigma, nvar, vars);
0235     } else {
0236       m_pos->selectPeakResiduals(nsigma, nvar, vars);
0237     }
0238   }
0239
0240   void correctBField() {
0241     m_pos->correctBField();
0242     //if (m_twoBin) m_neg->correctBField();
0243   };
0244
0245   void fiducialCuts() { m_pos->fiducialCuts(); };
0246
0247   void eraseNotSelectedResiduals() {
0248     if (m_twoBin) {
0249       m_pos->eraseNotSelectedResiduals();
0250       m_neg->eraseNotSelectedResiduals();
0251     } else {
0252       m_pos->eraseNotSelectedResiduals();
0253     }
0254   }
0255
0256   std::vector<double *>::const_iterator residualsPos_begin() const { return m_pos->residuals_begin(); };
0257   std::vector<double *>::const_iterator residualsPos_end() const { return m_pos->residuals_end(); };
0258   std::vector<double *>::const_iterator residualsNeg_begin() const { return m_neg->residuals_begin(); };
0259   std::vector<double *>::const_iterator residualsNeg_end() const { return m_neg->residuals_end(); };
0260
0261   std::vector<bool>::const_iterator residualsPos_ok_begin() const { return m_pos->selectedResidualsFlags().begin(); };
0262   std::vector<bool>::const_iterator residualsPos_ok_end() const { return m_pos->selectedResidualsFlags().end(); };
0263   std::vector<bool>::const_iterator residualsNeg_ok_begin() const { return m_neg->selectedResidualsFlags().begin(); };
0264   std::vector<bool>::const_iterator residualsNeg_ok_end() const { return m_neg->selectedResidualsFlags().end(); };
0265
0266 protected:
0267   bool m_twoBin;
0268   MuonResidualsFitter *m_pos, *m_neg;
0269 };
0270
0271 #endif  // Alignment_MuonAlignmentAlgorithms_MuonResidualsTwoBin_H