CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/Alignment/MuonAlignmentAlgorithms/interface/MuonResidualsTwoBin.h

Go to the documentation of this file.
00001 #ifndef Alignment_MuonAlignmentAlgorithms_MuonResidualsTwoBin_H
00002 #define Alignment_MuonAlignmentAlgorithms_MuonResidualsTwoBin_H
00003 
00010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00011 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResidualsFitter.h"
00012 #include "Alignment/CommonAlignment/interface/Alignable.h"
00013 #include "TMath.h"
00014 
00015 class MuonResidualsTwoBin {
00016 public:
00017   MuonResidualsTwoBin(bool twoBin, MuonResidualsFitter *pos, MuonResidualsFitter *neg): m_twoBin(twoBin), m_pos(pos), m_neg(neg) {};
00018   ~MuonResidualsTwoBin() {
00019     if (m_pos != NULL) delete m_pos;
00020     if (m_neg != NULL) delete m_neg;
00021   };
00022 
00023   int residualsModel() const { assert(m_pos->residualsModel() == m_neg->residualsModel());  return m_pos->residualsModel(); };
00024   long numResidualsPos() const { return m_pos->numResiduals(); };
00025   long numResidualsNeg() const { return m_neg->numResiduals(); };
00026   int npar() { assert(m_pos->npar() == m_neg->npar());  return m_pos->npar(); };
00027   int ndata() { assert(m_pos->ndata() == m_neg->ndata());  return m_pos->ndata(); };
00028   int type() const { assert(m_pos->type() == m_neg->type());  return m_pos->type(); };
00029   int useRes() const { return m_pos->useRes(); };
00030 
00031   void fix(int parNum, bool value=true) {
00032     m_pos->fix(parNum, value);
00033     m_neg->fix(parNum, value);
00034   };
00035 
00036   bool fixed(int parNum) {
00037     return m_pos->fixed(parNum)  &&  m_neg->fixed(parNum);
00038   };
00039 
00040   void setPrintLevel(int printLevel) const {
00041     m_pos->setPrintLevel(printLevel);
00042     m_neg->setPrintLevel(printLevel);
00043   }
00044 
00045   void setStrategy(int strategy) const {
00046     m_pos->setStrategy(strategy);
00047     m_neg->setStrategy(strategy);
00048   }
00049 
00050   void fill(char charge, double *residual) {
00051     if (!m_twoBin  ||  charge > 0) m_pos->fill(residual);
00052     else m_neg->fill(residual);
00053   };
00054 
00055   bool fit(Alignable *ali) {
00056     return (m_twoBin ? (m_pos->fit(ali)  &&  m_neg->fit(ali)) : m_pos->fit(ali));
00057   };
00058   double value(int parNum) {
00059     return (m_twoBin ? ((m_pos->value(parNum) + m_neg->value(parNum)) / 2.) : m_pos->value(parNum));
00060   };
00061   double errorerror(int parNum) {
00062     return (m_twoBin ? (sqrt(pow(m_pos->errorerror(parNum), 2.) + pow(m_neg->errorerror(parNum), 2.)) / 2.) : m_pos->errorerror(parNum));
00063   };
00064   double antisym(int parNum) {
00065     return (m_twoBin ? ((m_pos->value(parNum) - m_neg->value(parNum)) / 2.) : 0.);
00066   };
00067   double loglikelihood() {
00068     return (m_twoBin ? (m_pos->loglikelihood() + m_neg->loglikelihood()) : m_pos->loglikelihood());
00069   };
00070   double numsegments() {
00071     return (m_twoBin ? (m_pos->numsegments() + m_neg->numsegments()) : m_pos->numsegments());
00072   };
00073   double sumofweights() {
00074     return (m_twoBin ? (m_pos->sumofweights() + m_neg->sumofweights()) : m_pos->sumofweights());
00075   };
00076 
00077   // demonstration plots
00078   double plot(std::string name, TFileDirectory *dir, Alignable *ali) {
00079     if (m_twoBin) {
00080       std::string namePos = name + std::string("Pos");
00081       std::string nameNeg = name + std::string("Neg");
00082       double output = 0.;
00083       output += m_pos->plot(namePos, dir, ali);
00084       output += m_neg->plot(nameNeg, dir, ali);
00085       return output;
00086     }
00087     else {
00088       return m_pos->plot(name, dir, ali);
00089     }
00090   };
00091 
00092   // I/O of temporary files for collect mode
00093   void write(FILE *file, int which=0) {
00094     if (m_twoBin) {
00095       m_pos->write(file, 2*which);
00096       m_neg->write(file, 2*which + 1);
00097     }
00098     else {
00099       m_pos->write(file, which);
00100     }
00101   };
00102   void read(FILE *file, int which=0) {
00103     if (m_twoBin) {
00104       m_pos->read(file, 2*which);
00105       m_neg->read(file, 2*which + 1);
00106     }
00107     else {
00108       m_pos->read(file, which);
00109     }
00110   };
00111 
00112   double median(int which) {
00113     std::vector<double> residuals;
00114     for (std::vector<double*>::const_iterator r = residualsPos_begin();  r != residualsPos_end();  ++r) {
00115       residuals.push_back((*r)[which]);
00116     }
00117     if (m_twoBin) {
00118       for (std::vector<double*>::const_iterator r = residualsNeg_begin();  r != residualsNeg_end();  ++r) {
00119         residuals.push_back((*r)[which]);
00120       }
00121     }
00122     std::sort(residuals.begin(), residuals.end());
00123     int length = residuals.size();
00124     return residuals[length/2];
00125   };
00126 
00127   double mean(int which, double truncate) {
00128     double sum = 0.;
00129     double n = 0.;
00130     for (std::vector<double*>::const_iterator r = residualsPos_begin();  r != residualsPos_end();  ++r) {
00131       double value = (*r)[which];
00132       if (fabs(value) < truncate) {
00133         sum += value;
00134         n += 1.;
00135       }
00136     }
00137     if (m_twoBin) {
00138       for (std::vector<double*>::const_iterator r = residualsNeg_begin();  r != residualsNeg_end();  ++r) {
00139         double value = (*r)[which];
00140         if (fabs(value) < truncate) {
00141           sum += value;
00142           n += 1.;
00143         }
00144       }
00145     }
00146     return sum/n;
00147   };
00148 
00149   double wmean(int which, int whichredchi2, double truncate) {
00150     double sum = 0.;
00151     double n = 0.;
00152     for (std::vector<double*>::const_iterator r = residualsPos_begin();  r != residualsPos_end();  ++r) {
00153       double value = (*r)[which];
00154       if (fabs(value) < truncate) {
00155         double weight = 1./(*r)[whichredchi2];
00156         if (TMath::Prob(1./weight*12, 12) < 0.99) {
00157           sum += weight*value;
00158           n += weight;
00159         }
00160       }
00161     }
00162     if (m_twoBin) {
00163       for (std::vector<double*>::const_iterator r = residualsNeg_begin();  r != residualsNeg_end();  ++r) {
00164         double value = (*r)[which];
00165         if (fabs(value) < truncate) {
00166           double weight = 1./(*r)[whichredchi2];
00167           if (TMath::Prob(1./weight*12, 12) < 0.99) {
00168             sum += weight*value;
00169             n += weight;
00170           }
00171         }
00172       }
00173     }
00174     return sum/n;
00175   };
00176 
00177   double stdev(int which, double truncate) {
00178     double sum2 = 0.;
00179     double sum = 0.;
00180     double n = 0.;
00181     for (std::vector<double*>::const_iterator r = residualsPos_begin();  r != residualsPos_end();  ++r) {
00182       double value = (*r)[which];
00183       if (fabs(value) < truncate) {
00184         sum2 += value*value;
00185         sum += value;
00186         n += 1.;
00187       }
00188     }
00189     if (m_twoBin) {
00190       for (std::vector<double*>::const_iterator r = residualsNeg_begin();  r != residualsNeg_end();  ++r) {
00191         double value = (*r)[which];
00192         if (fabs(value) < truncate) {
00193           sum2 += value*value;
00194           sum += value;
00195           n += 1.;
00196         }
00197       }
00198     }
00199     return sqrt(sum2/n - pow(sum/n, 2));
00200   };
00201 
00202   void plotsimple(std::string name, TFileDirectory *dir, int which, double multiplier) {
00203     if (m_twoBin) {
00204       std::string namePos = name + std::string("Pos");
00205       std::string nameNeg = name + std::string("Neg");
00206       m_pos->plotsimple(namePos, dir, which, multiplier);
00207       m_neg->plotsimple(nameNeg, dir, which, multiplier);
00208     }
00209     else {
00210       m_pos->plotsimple(name, dir, which, multiplier);
00211     }
00212   };
00213 
00214   void plotweighted(std::string name, TFileDirectory *dir, int which, int whichredchi2, double multiplier) {
00215     if (m_twoBin) {
00216       std::string namePos = name + std::string("Pos");
00217       std::string nameNeg = name + std::string("Neg");
00218       m_pos->plotweighted(namePos, dir, which, whichredchi2, multiplier);
00219       m_neg->plotweighted(nameNeg, dir, which, whichredchi2, multiplier);
00220     }
00221     else {
00222       m_pos->plotweighted(name, dir, which, whichredchi2, multiplier);
00223     }
00224   };
00225 
00226   void selectPeakResiduals(double nsigma, int nvar, int *vars)
00227   {
00228     if (m_twoBin) {
00229       m_pos->selectPeakResiduals(nsigma, nvar, vars);
00230       m_neg->selectPeakResiduals(nsigma, nvar, vars);
00231     }
00232     else {
00233       m_pos->selectPeakResiduals(nsigma, nvar, vars);
00234     }
00235   }
00236   
00237   void correctBField()
00238   {
00239     m_pos->correctBField();
00240     //if (m_twoBin) m_neg->correctBField();
00241   };
00242 
00243   void eraseNotSelectedResiduals()
00244   {
00245     if (m_twoBin) {
00246       m_pos->eraseNotSelectedResiduals();
00247       m_neg->eraseNotSelectedResiduals();
00248     }
00249     else {
00250       m_pos->eraseNotSelectedResiduals();
00251     }
00252   }
00253 
00254   std::vector<double*>::const_iterator residualsPos_begin() const { return m_pos->residuals_begin(); };
00255   std::vector<double*>::const_iterator residualsPos_end() const { return m_pos->residuals_end(); };
00256   std::vector<double*>::const_iterator residualsNeg_begin() const { return m_neg->residuals_begin(); };
00257   std::vector<double*>::const_iterator residualsNeg_end() const { return m_neg->residuals_end(); };
00258 
00259   std::vector<bool>::const_iterator residualsPos_ok_begin() const { return m_pos->selectedResidualsFlags().begin(); };
00260   std::vector<bool>::const_iterator residualsPos_ok_end() const { return m_pos->selectedResidualsFlags().end(); };
00261   std::vector<bool>::const_iterator residualsNeg_ok_begin() const { return m_neg->selectedResidualsFlags().begin(); };
00262   std::vector<bool>::const_iterator residualsNeg_ok_end() const { return m_neg->selectedResidualsFlags().end(); };
00263 
00264 protected:
00265   bool m_twoBin;
00266   MuonResidualsFitter *m_pos, *m_neg;
00267 };
00268 
00269 #endif // Alignment_MuonAlignmentAlgorithms_MuonResidualsTwoBin_H