File indexing completed on 2024-09-07 04:34:33
0001 #ifndef Alignment_MuonAlignmentAlgorithms_MuonResidualsTwoBin_H
0002 #define Alignment_MuonAlignmentAlgorithms_MuonResidualsTwoBin_H
0003
0004
0005
0006
0007
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
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
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
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
