Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-26 05:07:18 

0001 #include "TopQuarkAnalysis/TopKinFitter/interface/TtFullLepKinSolver.h"
0002 #include "TF2.h"
0003 
0004 TtFullLepKinSolver::TtFullLepKinSolver()
0005     : topmass_begin(0), topmass_end(0), topmass_step(0), mw(80.4), mb(4.8), pxmiss_(0), pymiss_(0) {
0006   // That crude parametrisation has been obtained from a fit of O(1000) pythia events.
0007   // It is normalized to 1.
0008   EventShape_ = new TF2("landau2D", "[0]*TMath::Landau(x,[1],[2],0)*TMath::Landau(y,[3],[4],0)", 0, 500, 0, 500);
0009   EventShape_->SetParameters(30.7137, 56.2880, 23.0744, 59.1015, 24.9145);
0010 }
0011 
0012 TtFullLepKinSolver::TtFullLepKinSolver(
0013     const double b, const double e, const double s, const std::vector<double>& nupars, const double mW, const double mB)
0014     : topmass_begin(b), topmass_end(e), topmass_step(s), mw(mW), mb(mB), pxmiss_(0), pymiss_(0) {
0015   EventShape_ = new TF2("landau2D", "[0]*TMath::Landau(x,[1],[2],0)*TMath::Landau(y,[3],[4],0)", 0, 500, 0, 500);
0016   EventShape_->SetParameters(nupars[0], nupars[1], nupars[2], nupars[3], nupars[4]);
0017 }
0018 
0019 //
0020 // destructor
0021 //
0022 TtFullLepKinSolver::~TtFullLepKinSolver() { delete EventShape_; }
0023 
0024 TtDilepEvtSolution TtFullLepKinSolver::addKinSolInfo(TtDilepEvtSolution* asol) {
0025   TtDilepEvtSolution fitsol(*asol);
0026 
0027   //antilepton and lepton
0028   TLorentzVector LV_e, LV_e_;
0029   //b and bbar quark
0030   TLorentzVector LV_b, LV_b_;
0031 
0032   bool hasMCinfo = true;
0033   if (fitsol.getGenN()) {  // protect against non-dilept genevents
0034     genLV_n = TLorentzVector(
0035         fitsol.getGenN()->px(), fitsol.getGenN()->py(), fitsol.getGenN()->pz(), fitsol.getGenN()->energy());
0036   } else
0037     hasMCinfo = false;
0038 
0039   if (fitsol.getGenNbar()) {  // protect against non-dilept genevents
0040     genLV_n_ = TLorentzVector(
0041         fitsol.getGenNbar()->px(), fitsol.getGenNbar()->py(), fitsol.getGenNbar()->pz(), fitsol.getGenNbar()->energy());
0042   } else
0043     hasMCinfo = false;
0044   // if MC is to be used to select the best top mass and is not available,
0045   // then nothing can be done. Stop here.
0046   if (useMCforBest_ && !hasMCinfo)
0047     return fitsol;
0048 
0049   // first lepton
0050   if (fitsol.getWpDecay() == "muon") {
0051     LV_e = TLorentzVector(
0052         fitsol.getMuonp().px(), fitsol.getMuonp().py(), fitsol.getMuonp().pz(), fitsol.getMuonp().energy());
0053   } else if (fitsol.getWpDecay() == "electron") {
0054     LV_e = TLorentzVector(fitsol.getElectronp().px(),
0055                           fitsol.getElectronp().py(),
0056                           fitsol.getElectronp().pz(),
0057                           fitsol.getElectronp().energy());
0058   } else if (fitsol.getWpDecay() == "tau") {
0059     LV_e =
0060         TLorentzVector(fitsol.getTaup().px(), fitsol.getTaup().py(), fitsol.getTaup().pz(), fitsol.getTaup().energy());
0061   }
0062 
0063   // second lepton
0064   if (fitsol.getWmDecay() == "muon") {
0065     LV_e_ = TLorentzVector(
0066         fitsol.getMuonm().px(), fitsol.getMuonm().py(), fitsol.getMuonm().pz(), fitsol.getMuonm().energy());
0067   } else if (fitsol.getWmDecay() == "electron") {
0068     LV_e_ = TLorentzVector(fitsol.getElectronm().px(),
0069                            fitsol.getElectronm().py(),
0070                            fitsol.getElectronm().pz(),
0071                            fitsol.getElectronm().energy());
0072   } else if (fitsol.getWmDecay() == "tau") {
0073     LV_e_ =
0074         TLorentzVector(fitsol.getTaum().px(), fitsol.getTaum().py(), fitsol.getTaum().pz(), fitsol.getTaum().energy());
0075   }
0076 
0077   // first jet
0078   LV_b = TLorentzVector(
0079       fitsol.getCalJetB().px(), fitsol.getCalJetB().py(), fitsol.getCalJetB().pz(), fitsol.getCalJetB().energy());
0080 
0081   // second jet
0082   LV_b_ = TLorentzVector(fitsol.getCalJetBbar().px(),
0083                          fitsol.getCalJetBbar().py(),
0084                          fitsol.getCalJetBbar().pz(),
0085                          fitsol.getCalJetBbar().energy());
0086 
0087   //loop on top mass parameter
0088   double weightmax = -1e30;
0089   double mtmax = 0;
0090   for (double mt = topmass_begin; mt < topmass_end + 0.5 * topmass_step; mt += topmass_step) {
0091     //cout << "mt = " << mt << endl;
0092     double q_coeff[5], q_sol[4];
0093     FindCoeff(LV_e, LV_e_, LV_b, LV_b_, mt, mt, pxmiss_, pymiss_, q_coeff);
0094     int NSol = quartic(q_coeff, q_sol);
0095 
0096     //loop on all solutions
0097     for (int isol = 0; isol < NSol; isol++) {
0098       TopRec(LV_e, LV_e_, LV_b, LV_b_, q_sol[isol]);
0099       double weight = useMCforBest_ ? WeightSolfromMC() : WeightSolfromShape();
0100       if (weight > weightmax) {
0101         weightmax = weight;
0102         mtmax = mt;
0103       }
0104     }
0105 
0106     //for (int i=0;i<5;i++) cout << " q_coeff["<<i<< "]= " << q_coeff[i];
0107     //cout << endl;
0108 
0109     //for (int i=0;i<4;i++) cout << " q_sol["<<i<< "]= " << q_sol[i];
0110     //cout << endl;
0111     //cout << "NSol_" << NSol << endl;
0112   }
0113 
0114   fitsol.setRecTopMass(mtmax);
0115   fitsol.setRecWeightMax(weightmax);
0116 
0117   return fitsol;
0118 }
0119 
0120 void TtFullLepKinSolver::SetConstraints(const double xx, const double yy) {
0121   pxmiss_ = xx;
0122   pymiss_ = yy;
0123 }
0124 
0125 TtFullLepKinSolver::NeutrinoSolution TtFullLepKinSolver::getNuSolution(const TLorentzVector& LV_l,
0126                                                                        const TLorentzVector& LV_l_,
0127                                                                        const TLorentzVector& LV_b,
0128                                                                        const TLorentzVector& LV_b_) {
0129   math::XYZTLorentzVector maxLV_n = math::XYZTLorentzVector(0, 0, 0, 0);
0130   math::XYZTLorentzVector maxLV_n_ = math::XYZTLorentzVector(0, 0, 0, 0);
0131 
0132   //loop on top mass parameter
0133   double weightmax = -1;
0134   for (double mt = topmass_begin; mt < topmass_end + 0.5 * topmass_step; mt += topmass_step) {
0135     double q_coeff[5], q_sol[4];
0136     FindCoeff(LV_l, LV_l_, LV_b, LV_b_, mt, mt, pxmiss_, pymiss_, q_coeff);
0137     int NSol = quartic(q_coeff, q_sol);
0138 
0139     //loop on all solutions
0140     for (int isol = 0; isol < NSol; isol++) {
0141       TopRec(LV_l, LV_l_, LV_b, LV_b_, q_sol[isol]);
0142       double weight = WeightSolfromShape();
0143       if (weight > weightmax) {
0144         weightmax = weight;
0145         maxLV_n.SetPxPyPzE(LV_n.Px(), LV_n.Py(), LV_n.Pz(), LV_n.E());
0146         maxLV_n_.SetPxPyPzE(LV_n_.Px(), LV_n_.Py(), LV_n_.Pz(), LV_n_.E());
0147       }
0148     }
0149   }
0150   TtFullLepKinSolver::NeutrinoSolution nuSol;
0151   nuSol.neutrino = reco::LeafCandidate(0, maxLV_n);
0152   nuSol.neutrinoBar = reco::LeafCandidate(0, maxLV_n_);
0153   nuSol.weight = weightmax;
0154   return nuSol;
0155 }
0156 
0157 void TtFullLepKinSolver::FindCoeff(const TLorentzVector& al,
0158                                    const TLorentzVector& l,
0159                                    const TLorentzVector& b_al,
0160                                    const TLorentzVector& b_l,
0161                                    const double mt,
0162                                    const double mat,
0163                                    const double px_miss,
0164                                    const double py_miss,
0165                                    double* koeficienty) {
0166   double E, apom1, apom2, apom3;
0167   double k11, k21, k31, k41, cpom1, cpom2, cpom3, l11, l21, l31, l41, l51, l61, k1, k2, k3, k4, k5, k6;
0168   double l1, l2, l3, l4, l5, l6, k15, k25, k35, k45;
0169 
0170   C = -al.Px() - b_al.Px() - l.Px() - b_l.Px() + px_miss;
0171   D = -al.Py() - b_al.Py() - l.Py() - b_l.Py() + py_miss;
0172 
0173   // right side of first two linear equations - missing pT
0174 
0175   E = (sqr(mt) - sqr(mw) - sqr(mb)) / (2 * b_al.E()) - sqr(mw) / (2 * al.E()) - al.E() +
0176       al.Px() * b_al.Px() / b_al.E() + al.Py() * b_al.Py() / b_al.E() + al.Pz() * b_al.Pz() / b_al.E();
0177   F = (sqr(mat) - sqr(mw) - sqr(mb)) / (2 * b_l.E()) - sqr(mw) / (2 * l.E()) - l.E() + l.Px() * b_l.Px() / b_l.E() +
0178       l.Py() * b_l.Py() / b_l.E() + l.Pz() * b_l.Pz() / b_l.E();
0179 
0180   m1 = al.Px() / al.E() - b_al.Px() / b_al.E();
0181   m2 = al.Py() / al.E() - b_al.Py() / b_al.E();
0182   m3 = al.Pz() / al.E() - b_al.Pz() / b_al.E();
0183 
0184   n1 = l.Px() / l.E() - b_l.Px() / b_l.E();
0185   n2 = l.Py() / l.E() - b_l.Py() / b_l.E();
0186   n3 = l.Pz() / l.E() - b_l.Pz() / b_l.E();
0187 
0188   pom = E - m1 * C - m2 * D;
0189   apom1 = sqr(al.Px()) - sqr(al.E());
0190   apom2 = sqr(al.Py()) - sqr(al.E());
0191   apom3 = sqr(al.Pz()) - sqr(al.E());
0192 
0193   k11 = 1 / sqr(al.E()) *
0194         (pow(mw, 4) / 4 + sqr(C) * apom1 + sqr(D) * apom2 + apom3 * sqr(pom) / sqr(m3) +
0195          sqr(mw) * (al.Px() * C + al.Py() * D + al.Pz() * pom / m3) + 2 * al.Px() * al.Py() * C * D +
0196          2 * al.Px() * al.Pz() * C * pom / m3 + 2 * al.Py() * al.Pz() * D * pom / m3);
0197   k21 = 1 / sqr(al.E()) *
0198         (-2 * C * m3 * n3 * apom1 + 2 * apom3 * n3 * m1 * pom / m3 - sqr(mw) * m3 * n3 * al.Px() +
0199          sqr(mw) * m1 * n3 * al.Pz() - 2 * al.Px() * al.Py() * D * m3 * n3 + 2 * al.Px() * al.Pz() * C * m1 * n3 -
0200          2 * al.Px() * al.Pz() * n3 * pom + 2 * al.Py() * al.Pz() * D * m1 * n3);
0201   k31 = 1 / sqr(al.E()) *
0202         (-2 * D * m3 * n3 * apom2 + 2 * apom3 * n3 * m2 * pom / m3 - sqr(mw) * m3 * n3 * al.Py() +
0203          sqr(mw) * m2 * n3 * al.Pz() - 2 * al.Px() * al.Py() * C * m3 * n3 + 2 * al.Px() * al.Pz() * C * m2 * n3 -
0204          2 * al.Py() * al.Pz() * n3 * pom + 2 * al.Py() * al.Pz() * D * m2 * n3);
0205   k41 = 1 / sqr(al.E()) *
0206         (2 * apom3 * m1 * m2 * sqr(n3) + 2 * al.Px() * al.Py() * sqr(m3) * sqr(n3) -
0207          2 * al.Px() * al.Pz() * m2 * m3 * sqr(n3) - 2 * al.Py() * al.Pz() * m1 * m3 * sqr(n3));
0208   k51 = 1 / sqr(al.E()) *
0209         (apom1 * sqr(m3) * sqr(n3) + apom3 * sqr(m1) * sqr(n3) - 2 * al.Px() * al.Pz() * m1 * m3 * sqr(n3));
0210   k61 = 1 / sqr(al.E()) *
0211         (apom2 * sqr(m3) * sqr(n3) + apom3 * sqr(m2) * sqr(n3) - 2 * al.Py() * al.Pz() * m2 * m3 * sqr(n3));
0212 
0213   cpom1 = sqr(l.Px()) - sqr(l.E());
0214   cpom2 = sqr(l.Py()) - sqr(l.E());
0215   cpom3 = sqr(l.Pz()) - sqr(l.E());
0216 
0217   l11 = 1 / sqr(l.E()) * (pow(mw, 4) / 4 + cpom3 * sqr(F) / sqr(n3) + sqr(mw) * l.Pz() * F / n3);
0218   l21 =
0219       1 / sqr(l.E()) *
0220       (-2 * cpom3 * F * m3 * n1 / n3 + sqr(mw) * (l.Px() * m3 * n3 - l.Pz() * n1 * m3) + 2 * l.Px() * l.Pz() * F * m3);
0221   l31 =
0222       1 / sqr(l.E()) *
0223       (-2 * cpom3 * F * m3 * n2 / n3 + sqr(mw) * (l.Py() * m3 * n3 - l.Pz() * n2 * m3) + 2 * l.Py() * l.Pz() * F * m3);
0224   l41 = 1 / sqr(l.E()) *
0225         (2 * cpom3 * n1 * n2 * sqr(m3) + 2 * l.Px() * l.Py() * sqr(m3) * sqr(n3) -
0226          2 * l.Px() * l.Pz() * n2 * n3 * sqr(m3) - 2 * l.Py() * l.Pz() * n1 * n3 * sqr(m3));
0227   l51 = 1 / sqr(l.E()) *
0228         (cpom1 * sqr(m3) * sqr(n3) + cpom3 * sqr(n1) * sqr(m3) - 2 * l.Px() * l.Pz() * n1 * n3 * sqr(m3));
0229   l61 = 1 / sqr(l.E()) *
0230         (cpom2 * sqr(m3) * sqr(n3) + cpom3 * sqr(n2) * sqr(m3) - 2 * l.Py() * l.Pz() * n2 * n3 * sqr(m3));
0231 
0232   k1 = k11 * k61;
0233   k2 = k61 * k21 / k51;
0234   k3 = k31;
0235   k4 = k41 / k51;
0236   k5 = k61 / k51;
0237   k6 = 1;
0238 
0239   l1 = l11 * k61;
0240   l2 = l21 * k61 / k51;
0241   l3 = l31;
0242   l4 = l41 / k51;
0243   l5 = l51 * k61 / (sqr(k51));
0244   l6 = l61 / k61;
0245 
0246   k15 = k1 * l5 - l1 * k5;
0247   k25 = k2 * l5 - l2 * k5;
0248   k35 = k3 * l5 - l3 * k5;
0249   k45 = k4 * l5 - l4 * k5;
0250 
0251   k16 = k1 * l6 - l1 * k6;
0252   k26 = k2 * l6 - l2 * k6;
0253   k36 = k3 * l6 - l3 * k6;
0254   k46 = k4 * l6 - l4 * k6;
0255   k56 = k5 * l6 - l5 * k6;
0256 
0257   koeficienty[0] = k15 * sqr(k36) - k35 * k36 * k16 - k56 * sqr(k16);
0258   koeficienty[1] =
0259       2 * k15 * k36 * k46 + k25 * sqr(k36) + k35 * (-k46 * k16 - k36 * k26) - k45 * k36 * k16 - 2 * k56 * k26 * k16;
0260   koeficienty[2] = k15 * sqr(k46) + 2 * k25 * k36 * k46 + k35 * (-k46 * k26 - k36 * k56) -
0261                    k56 * (sqr(k26) + 2 * k56 * k16) - k45 * (k46 * k16 + k36 * k26);
0262   koeficienty[3] = k25 * sqr(k46) - k35 * k46 * k56 - k45 * (k46 * k26 + k36 * k56) - 2 * sqr(k56) * k26;
0263   koeficienty[4] = -k45 * k46 * k56 - pow(k56, 3);
0264 
0265   // normalization of coefficients
0266   int moc = (int(log10(fabs(koeficienty[0]))) + int(log10(fabs(koeficienty[4])))) / 2;
0267 
0268   koeficienty[0] = koeficienty[0] / TMath::Power(10, moc);
0269   koeficienty[1] = koeficienty[1] / TMath::Power(10, moc);
0270   koeficienty[2] = koeficienty[2] / TMath::Power(10, moc);
0271   koeficienty[3] = koeficienty[3] / TMath::Power(10, moc);
0272   koeficienty[4] = koeficienty[4] / TMath::Power(10, moc);
0273 }
0274 
0275 void TtFullLepKinSolver::TopRec(const TLorentzVector& al,
0276                                 const TLorentzVector& l,
0277                                 const TLorentzVector& b_al,
0278                                 const TLorentzVector& b_l,
0279                                 const double sol) {
0280   TVector3 t_ttboost;
0281   TLorentzVector aux;
0282   double pxp, pyp, pzp, pup, pvp, pwp;
0283 
0284   pxp = sol * (m3 * n3 / k51);
0285   pyp = -(m3 * n3 / k61) * (k56 * pow(sol, 2) + k26 * sol + k16) / (k36 + k46 * sol);
0286   pzp = -1 / n3 * (n1 * pxp + n2 * pyp - F);
0287   pwp = 1 / m3 * (m1 * pxp + m2 * pyp + pom);
0288   pup = C - pxp;
0289   pvp = D - pyp;
0290 
0291   LV_n_.SetXYZM(pxp, pyp, pzp, 0.0);
0292   LV_n.SetXYZM(pup, pvp, pwp, 0.0);
0293 
0294   LV_t_ = b_l + l + LV_n_;
0295   LV_t = b_al + al + LV_n;
0296 
0297   aux = (LV_t_ + LV_t);
0298   t_ttboost = -aux.BoostVector();
0299   LV_tt_t_ = LV_t_;
0300   LV_tt_t = LV_t;
0301   LV_tt_t_.Boost(t_ttboost);
0302   LV_tt_t.Boost(t_ttboost);
0303 }
0304 
0305 double TtFullLepKinSolver::WeightSolfromMC() const {
0306   double weight = 1;
0307   weight = ((LV_n.E() > genLV_n.E()) ? genLV_n.E() / LV_n.E() : LV_n.E() / genLV_n.E()) *
0308            ((LV_n_.E() > genLV_n_.E()) ? genLV_n_.E() / LV_n_.E() : LV_n_.E() / genLV_n_.E());
0309   return weight;
0310 }
0311 
0312 double TtFullLepKinSolver::WeightSolfromShape() const { return EventShape_->Eval(LV_n.E(), LV_n_.E()); }
0313 
0314 int TtFullLepKinSolver::quartic(double* koeficienty, double* koreny) const {
0315   double w, b0, b1, b2;
0316   double c[4];
0317   double d0, d1, h, t, z;
0318   double* px;
0319 
0320   if (koeficienty[4] == 0.0) {
0321     return cubic(koeficienty, koreny);
0322   }
0323   /* quartic problem? */
0324   w = koeficienty[3] / (4 * koeficienty[4]);
0325   /* offset */
0326   b2 = -6 * sqr(w) + koeficienty[2] / koeficienty[4];
0327   /* koeficienty. of shifted polynomial */
0328   b1 = (8 * sqr(w) - 2 * koeficienty[2] / koeficienty[4]) * w + koeficienty[1] / koeficienty[4];
0329   b0 = ((-3 * sqr(w) + koeficienty[2] / koeficienty[4]) * w - koeficienty[1] / koeficienty[4]) * w +
0330        koeficienty[0] / koeficienty[4];
0331 
0332   c[3] = 1.0;
0333   /* cubic resolvent */
0334   c[2] = b2;
0335   c[1] = -4 * b0;
0336   c[0] = sqr(b1) - 4 * b0 * b2;
0337 
0338   if (cubic(c, koreny) == 0) {
0339     // No real solutions, returning zero
0340     return 0;
0341   } else {
0342     z = koreny[0];
0343   }
0344   //double z1=1.0,z2=2.0,z3=3.0;
0345   //TMath::RootsCubic(c,z1,z2,z3);
0346   //if (z2 !=0) z = z2;
0347   //if (z1 !=0) z = z1;
0348   /* only lowermost root needed */
0349 
0350   int nreal = 0;
0351   px = koreny;
0352   t = sqrt(0.25 * sqr(z) - b0);
0353   for (int i = -1; i <= 1; i += 2) {
0354     d0 = -0.5 * z + i * t;
0355     /* coeffs. of quadratic factor */
0356     d1 = (t != 0.0) ? -i * 0.5 * b1 / t : i * sqrt(-z - b2);
0357     h = 0.25 * sqr(d1) - d0;
0358     if (h >= 0.0) {
0359       h = sqrt(h);
0360       nreal += 2;
0361       *px++ = -0.5 * d1 - h - w;
0362       *px++ = -0.5 * d1 + h - w;
0363     }
0364   }
0365 
0366   //  if (nreal==4) {
0367   /* sort results */
0368   //    if (koreny[2]<koreny[0]) SWAP(koreny[0], koreny[2]);
0369   //    if (koreny[3]<koreny[1]) SWAP(koreny[1], koreny[3]);
0370   //    if (koreny[1]<koreny[0]) SWAP(koreny[0], koreny[1]);
0371   //    if (koreny[3]<koreny[2]) SWAP(koreny[2], koreny[3]);
0372   //    if (koreny[2]<koreny[1]) SWAP(koreny[1], koreny[2]);
0373   //  }
0374   return nreal;
0375 }
0376 
0377 unsigned int TtFullLepKinSolver::cubic(const double* coeffs, double* koreny) const {
0378   unsigned nreal;
0379   double w, p, q, dis, h, phi;
0380 
0381   if (coeffs[3] != 0.0) {
0382     /* cubic problem? */
0383     w = coeffs[2] / (3 * coeffs[3]);
0384     p = sqr(coeffs[1] / (3 * coeffs[3]) - sqr(w)) * (coeffs[1] / (3 * coeffs[3]) - sqr(w));
0385     q = -0.5 * (2 * sqr(w) * w - (coeffs[1] * w - coeffs[0]) / coeffs[3]);
0386     dis = sqr(q) + p;
0387     /* discriminant */
0388     if (dis < 0.0) {
0389       /* 3 real solutions */
0390       h = q / sqrt(-p);
0391       if (h > 1.0)
0392         h = 1.0;
0393       /* confine the argument of */
0394       if (h < -1.0)
0395         h = -1.0;
0396       /* acos to [-1;+1] */
0397       phi = acos(h);
0398       p = 2 * TMath::Power(-p, 1.0 / 6.0);
0399       for (unsigned i = 0; i < 3; i++)
0400         koreny[i] = p * cos((phi + 2 * i * TMath::Pi()) / 3.0) - w;
0401       if (koreny[1] < koreny[0])
0402         SWAP(koreny[0], koreny[1]);
0403       /* sort results */
0404       if (koreny[2] < koreny[1])
0405         SWAP(koreny[1], koreny[2]);
0406       if (koreny[1] < koreny[0])
0407         SWAP(koreny[0], koreny[1]);
0408       nreal = 3;
0409     } else {
0410       /* only one real solution */
0411       dis = sqrt(dis);
0412       h = TMath::Power(fabs(q + dis), 1.0 / 3.0);
0413       p = TMath::Power(fabs(q - dis), 1.0 / 3.0);
0414       koreny[0] = ((q + dis > 0.0) ? h : -h) + ((q - dis > 0.0) ? p : -p) - w;
0415       nreal = 1;
0416     }
0417 
0418     /* Perform one step of a Newton iteration in order to minimize
0419        round-off errors */
0420     for (unsigned i = 0; i < nreal; i++) {
0421       h = coeffs[1] + koreny[i] * (2 * coeffs[2] + 3 * koreny[i] * coeffs[3]);
0422       if (h != 0.0)
0423         koreny[i] -= (coeffs[0] + koreny[i] * (coeffs[1] + koreny[i] * (coeffs[2] + koreny[i] * coeffs[3]))) / h;
0424     }
0425   }
0426 
0427   else if (coeffs[2] != 0.0) {
0428     /* quadratic problem? */
0429     p = 0.5 * coeffs[1] / coeffs[2];
0430     dis = sqr(p) - coeffs[0] / coeffs[2];
0431     if (dis >= 0.0) {
0432       /* two real solutions */
0433       dis = sqrt(dis);
0434       koreny[0] = -p - dis;
0435       koreny[1] = -p + dis;
0436       nreal = 2;
0437     } else
0438       /* no real solution */
0439       nreal = 0;
0440   }
0441 
0442   else if (coeffs[1] != 0.0) {
0443     /* linear problem? */
0444     koreny[0] = -coeffs[0] / coeffs[1];
0445     nreal = 1;
0446   }
0447 
0448   else {
0449     /* no equation */
0450     nreal = 0;
0451   }
0452   return nreal;
0453 }
0454 
0455 void TtFullLepKinSolver::SWAP(double& realone, double& realtwo) const {
0456   if (realtwo < realone) {
0457     double aux = realtwo;
0458     realtwo = realone;
0459     realone = aux;
0460   }
0461 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!