Reference Manual

Go to the documentation of this file.
00001 //based on a code by Jan Valenta
00002 #include "TopQuarkAnalysis/TopKinFitter/interface/TtFullLepKinSolver.h"
00003 #include "TF2.h"
00004 
00005 TtFullLepKinSolver::TtFullLepKinSolver() 
00006 {
00007   // That crude parametrisation has been obtained from a fit of O(1000) pythia events.
00008   // It is normalized to 1.
00009   EventShape_ = new TF2("landau2D","[0]*TMath::Landau(x,[1],[2],0)*TMath::Landau(y,[3],[4],0)",0,500,0,500);
00010   EventShape_->SetParameters(30.7137,56.2880,23.0744,59.1015,24.9145);
00011 }
00012 
00013 TtFullLepKinSolver::TtFullLepKinSolver(double b, double e, double s, std::vector<double> nupars) 
00014 {
00015   topmass_begin = b;
00016   topmass_end = e;
00017   topmass_step = s;
00018   pxmiss_ = 0;
00019   pymiss_ = 0;
00020   mw = 80.22;
00021   maw = 80.22;
00022   mab = 4.800;
00023   mb = 4.800;
00024   // That crude parametrisation has been obtained from a fit of O(1000) pythia events.
00025   // It is normalized to 1.
00026   EventShape_ = new TF2("landau2D","[0]*TMath::Landau(x,[1],[2],0)*TMath::Landau(y,[3],[4],0)",0,500,0,500);
00027   EventShape_->SetParameters(nupars[0],nupars[1],nupars[2],nupars[3],nupars[4]);
00028 }
00029 
00030 //
00031 // destructor
00032 //
00033 TtFullLepKinSolver::~TtFullLepKinSolver() 
00034 {
00035   delete EventShape_;
00036 }
00037 
00038 
00039 
00040 TtDilepEvtSolution TtFullLepKinSolver::addKinSolInfo(TtDilepEvtSolution * asol) 
00041 {
00042   TtDilepEvtSolution fitsol(*asol);
00043   
00044   //antilepton and lepton
00045   TLorentzVector LV_e, LV_e_;
00046   //b and bbar quark
00047   TLorentzVector LV_b, LV_b_;
00048   
00049   bool hasMCinfo = true;
00050   if(fitsol.getGenN()) { // protect against non-dilept genevents
00051     genLV_n = TLorentzVector(fitsol.getGenN()->px(), fitsol.getGenN()->py(),
00052                              fitsol.getGenN()->pz(), fitsol.getGenN()->energy());
00053   } else hasMCinfo = false;
00054 
00055   if(fitsol.getGenNbar()) { // protect against non-dilept genevents
00056     genLV_n_ = TLorentzVector(fitsol.getGenNbar()->px(), fitsol.getGenNbar()->py(),
00057                               fitsol.getGenNbar()->pz(), fitsol.getGenNbar()->energy());
00058   } else hasMCinfo = false;
00059   // if MC is to be used to select the best top mass and is not available,
00060   // then nothing can be done. Stop here.
00061   if(useMCforBest_&&!hasMCinfo) return fitsol;
00062 
00063   // first lepton
00064   if (fitsol.getWpDecay() == "muon") {
00065     LV_e = TLorentzVector(fitsol.getMuonp().px(), fitsol.getMuonp().py(),
00066                           fitsol.getMuonp().pz(), fitsol.getMuonp().energy());
00067   } else if (fitsol.getWpDecay() == "electron") {
00068     LV_e = TLorentzVector(fitsol.getElectronp().px(), fitsol.getElectronp().py(),
00069                           fitsol.getElectronp().pz(), fitsol.getElectronp().energy());
00070   } else if (fitsol.getWpDecay() == "tau") {
00071     LV_e = TLorentzVector(fitsol.getTaup().px(), fitsol.getTaup().py(),
00072                           fitsol.getTaup().pz(), fitsol.getTaup().energy());
00073   }
00074     
00075   // second lepton
00076   if (fitsol.getWmDecay() == "muon") {
00077     LV_e_ = TLorentzVector(fitsol.getMuonm().px(), fitsol.getMuonm().py(),
00078                           fitsol.getMuonm().pz(), fitsol.getMuonm().energy());
00079   } else if (fitsol.getWmDecay() == "electron") {
00080     LV_e_ = TLorentzVector(fitsol.getElectronm().px(), fitsol.getElectronm().py(),
00081                            fitsol.getElectronm().pz(), fitsol.getElectronm().energy());
00082   } else if (fitsol.getWmDecay() == "tau") {
00083     LV_e_ = TLorentzVector(fitsol.getTaum().px(), fitsol.getTaum().py(),
00084                            fitsol.getTaum().pz(), fitsol.getTaum().energy());
00085   }
00086 
00087   // first jet
00088   LV_b = TLorentzVector(fitsol.getCalJetB().px(), fitsol.getCalJetB().py(),
00089                         fitsol.getCalJetB().pz(), fitsol.getCalJetB().energy());
00090 
00091   // second jet
00092   LV_b_ = TLorentzVector(fitsol.getCalJetBbar().px(), fitsol.getCalJetBbar().py(),
00093                          fitsol.getCalJetBbar().pz(), fitsol.getCalJetBbar().energy());
00094   
00095   //loop on top mass parameter
00096   double weightmax = -1e30;
00097   double mtmax = 0;
00098   for (double mt = topmass_begin; 
00099               mt < topmass_end + 0.5*topmass_step; 
00100               mt += topmass_step) {
00101     //cout << "mt = " << mt << endl;
00102     double q_coeff[5], q_sol[4];
00103     FindCoeff(LV_e, LV_e_, LV_b, LV_b_, mt, mt, pxmiss_, pymiss_, q_coeff);
00104     int NSol = quartic(q_coeff, q_sol);
00105     
00106     //loop on all solutions
00107     for (int isol = 0; isol < NSol; isol++) {
00108       TopRec(LV_e, LV_e_, LV_b, LV_b_, q_sol[isol]);
00109       double weight = useMCforBest_ ? WeightSolfromMC() : WeightSolfromShape();
00110       if (weight > weightmax) {
00111         weightmax =weight;
00112         mtmax = mt;
00113       }
00114     }
00115     
00116     //for (int i=0;i<5;i++) cout << " q_coeff["<<i<< "]= " << q_coeff[i];
00117     //cout << endl;
00118     
00119     //for (int i=0;i<4;i++) cout << " q_sol["<<i<< "]= " << q_sol[i];
00120     //cout << endl;
00121     //cout << "NSol_" << NSol << endl;
00122   }
00123   
00124   fitsol.setRecTopMass(mtmax);
00125   fitsol.setRecWeightMax(weightmax);
00126   
00127   return fitsol;
00128 }
00129 
00130 void TtFullLepKinSolver::SetConstraints(double xx, double yy)
00131 {
00132   pxmiss_ = xx;
00133   pymiss_ = yy;
00134 }
00135 
00136 TtFullLepKinSolver::NeutrinoSolution TtFullLepKinSolver::getNuSolution(TLorentzVector LV_l, 
00137                                                                        TLorentzVector LV_l_, 
00138                                                                        TLorentzVector LV_b, 
00139                                                                        TLorentzVector LV_b_) {                                                                                                                      
00140     
00141   math::XYZTLorentzVector maxLV_n  = math::XYZTLorentzVector(0,0,0,0); 
00142   math::XYZTLorentzVector maxLV_n_ = math::XYZTLorentzVector(0,0,0,0);   
00143 
00144   //loop on top mass parameter
00145   double weightmax = -1;
00146   for(double mt = topmass_begin; 
00147              mt < topmass_end + 0.5*topmass_step; 
00148              mt += topmass_step) {
00149     double q_coeff[5], q_sol[4];
00150     FindCoeff(LV_l, LV_l_, LV_b, LV_b_, mt, mt, pxmiss_, pymiss_, q_coeff);
00151     int NSol = quartic(q_coeff, q_sol);
00152     
00153     //loop on all solutions
00154     for (int isol = 0; isol < NSol; isol++) {
00155       TopRec(LV_l, LV_l_, LV_b, LV_b_, q_sol[isol]);
00156       double weight = WeightSolfromShape();
00157       if (weight > weightmax) {
00158         weightmax =weight;
00159         maxLV_n.SetPxPyPzE(LV_n.Px(), LV_n.Py(), LV_n.Pz(), LV_n.E());
00160         maxLV_n_.SetPxPyPzE(LV_n_.Px(), LV_n_.Py(), LV_n_.Pz(), LV_n_.E());
00161       }
00162     }
00163   }
00164   TtFullLepKinSolver::NeutrinoSolution nuSol;
00165   nuSol.neutrino    = reco::LeafCandidate(0, maxLV_n  );
00166   nuSol.neutrinoBar = reco::LeafCandidate(0, maxLV_n_ ); 
00167   nuSol.weight = weightmax; 
00168   return nuSol;
00169 }
00170 
00171 void TtFullLepKinSolver::FindCoeff(const TLorentzVector al, 
00172                                    const TLorentzVector l,
00173                                    const TLorentzVector b_al,
00174                                    const TLorentzVector b_l,
00175                                    double mt, 
00176                                    double mat, 
00177                                    double px_miss, 
00178                                    double py_miss,
00179                                    double* koeficienty) {
00180 
00181 double E, apom1, apom2, apom3;
00182 double k11, k21, k31, k41,  cpom1, cpom2, cpom3, l11, l21, l31, l41, l51, l61, k1, k2, k3, k4, k5,k6;
00183 double l1, l2, l3, l4, l5, l6, k15, k25, k35, k45;
00184 
00185   C = -al.Px()-b_al.Px()-l.Px()-b_l.Px() + px_miss;
00186   D = -al.Py()-b_al.Py()-l.Py()-b_l.Py() + py_miss;
00187 
00188   // right side offirst two linear equations - missing pT
00189 
00190   
00191   E = (sqr(mt)-sqr(mw)-sqr(mb))/(2*b_al.E())-sqr(mw)/(2*al.E())-al.E()+al.Px()*b_al.Px()/b_al.E()+al.Py()*b_al.Py()/b_al.E()+al.Pz()*b_al.Pz()/b_al.E();
00192   F = (sqr(mat)-sqr(maw)-sqr(mab))/(2*b_l.E())-sqr(maw)/(2*l.E())-l.E()+l.Px()*b_l.Px()/b_l.E()+l.Py()*b_l.Py()/b_l.E()+l.Pz()*b_l.Pz()/b_l.E();
00193   
00194   m1 = al.Px()/al.E()-b_al.Px()/b_al.E();
00195   m2 = al.Py()/al.E()-b_al.Py()/b_al.E();
00196   m3 = al.Pz()/al.E()-b_al.Pz()/b_al.E();
00197   
00198   n1 = l.Px()/l.E()-b_l.Px()/b_l.E();
00199   n2 = l.Py()/l.E()-b_l.Py()/b_l.E();
00200   n3 = l.Pz()/l.E()-b_l.Pz()/b_l.E();
00201   
00202   pom = E-m1*C-m2*D;
00203   apom1 = sqr(al.Px())-sqr(al.E());
00204   apom2 = sqr(al.Py())-sqr(al.E());
00205   apom3 = sqr(al.Pz())-sqr(al.E());
00206   
00207   k11 = 1/sqr(al.E())*(pow(mw,4)/4+sqr(C)*apom1+sqr(D)*apom2+apom3*sqr(pom)/sqr(m3)+sqr(mw)*(al.Px()*C+al.Py()*D+al.Pz()*pom/m3)+2*al.Px()*al.Py()*C*D+2*al.Px()*al.Pz()*C*pom/m3+2*al.Py()*al.Pz()*D*pom/m3);
00208   k21 = 1/sqr(al.E())*(-2*C*m3*n3*apom1+2*apom3*n3*m1*pom/m3-sqr(mw)*m3*n3*al.Px()+sqr(mw)*m1*n3*al.Pz()-2*al.Px()*al.Py()*D*m3*n3+2*al.Px()*al.Pz()*C*m1*n3-2*al.Px()*al.Pz()*n3*pom+2*al.Py()*al.Pz()*D*m1*n3);
00209   k31 = 1/sqr(al.E())*(-2*D*m3*n3*apom2+2*apom3*n3*m2*pom/m3-sqr(mw)*m3*n3*al.Py()+sqr(mw)*m2*n3*al.Pz()-2*al.Px()*al.Py()*C*m3*n3+2*al.Px()*al.Pz()*C*m2*n3-2*al.Py()*al.Pz()*n3*pom+2*al.Py()*al.Pz()*D*m2*n3);
00210   k41 = 1/sqr(al.E())*(2*apom3*m1*m2*sqr(n3)+2*al.Px()*al.Py()*sqr(m3)*sqr(n3)-2*al.Px()*al.Pz()*m2*m3*sqr(n3)-2*al.Py()*al.Pz()*m1*m3*sqr(n3));
00211   k51 = 1/sqr(al.E())*(apom1*sqr(m3)*sqr(n3)+apom3*sqr(m1)*sqr(n3)-2*al.Px()*al.Pz()*m1*m3*sqr(n3));
00212   k61 = 1/sqr(al.E())*(apom2*sqr(m3)*sqr(n3)+apom3*sqr(m2)*sqr(n3)-2*al.Py()*al.Pz()*m2*m3*sqr(n3));
00213   
00214   cpom1 = sqr(l.Px())-sqr(l.E());
00215   cpom2 = sqr(l.Py())-sqr(l.E());
00216   cpom3 = sqr(l.Pz())-sqr(l.E());
00217   
00218   l11 = 1/sqr(l.E())*(pow(maw,4)/4+cpom3*sqr(F)/sqr(n3)+sqr(maw)*l.Pz()*F/n3);
00219   l21 = 1/sqr(l.E())*(-2*cpom3*F*m3*n1/n3+sqr(maw)*(l.Px()*m3*n3-l.Pz()*n1*m3)+2*l.Px()*l.Pz()*F*m3);
00220   l31 = 1/sqr(l.E())*(-2*cpom3*F*m3*n2/n3+sqr(maw)*(l.Py()*m3*n3-l.Pz()*n2*m3)+2*l.Py()*l.Pz()*F*m3);
00221   l41 = 1/sqr(l.E())*(2*cpom3*n1*n2*sqr(m3)+2*l.Px()*l.Py()*sqr(m3)*sqr(n3)-2*l.Px()*l.Pz()*n2*n3*sqr(m3)-2*l.Py()*l.Pz()*n1*n3*sqr(m3));
00222   l51 = 1/sqr(l.E())*(cpom1*sqr(m3)*sqr(n3)+cpom3*sqr(n1)*sqr(m3)-2*l.Px()*l.Pz()*n1*n3*sqr(m3));
00223   l61 = 1/sqr(l.E())*(cpom2*sqr(m3)*sqr(n3)+cpom3*sqr(n2)*sqr(m3)-2*l.Py()*l.Pz()*n2*n3*sqr(m3));
00224   
00225   k1 = k11*k61;
00226   k2 = k61*k21/k51;
00227   k3 = k31;
00228   k4 = k41/k51;
00229   k5 = k61/k51;
00230   k6 = 1;
00231   
00232   l1 = l11*k61;
00233   l2 = l21*k61/k51;
00234   l3 = l31;
00235   l4 = l41/k51;
00236   l5 = l51*k61/(sqr(k51));
00237   l6 = l61/k61;
00238   
00239   k15 = k1*l5-l1*k5;
00240   k25 = k2*l5-l2*k5;
00241   k35 = k3*l5-l3*k5;
00242   k45 = k4*l5-l4*k5;
00243   
00244   k16 = k1*l6-l1*k6;
00245   k26 = k2*l6-l2*k6;
00246   k36 = k3*l6-l3*k6;
00247   k46 = k4*l6-l4*k6;
00248   k56 = k5*l6-l5*k6;
00249   
00250 
00251   koeficienty[0] = k15*sqr(k36)-k35*k36*k16-k56*sqr(k16);
00252   koeficienty[1] = 2*k15*k36*k46+k25*sqr(k36)+k35*(-k46*k16-k36*k26)-k45*k36*k16-2*k56*k26*k16;
00253   koeficienty[2] = k15*sqr(k46)+2*k25*k36*k46+k35*(-k46*k26-k36*k56)-k56*(sqr(k26)+2*k56*k16)-k45*(k46*k16+k36*k26);
00254   koeficienty[3] = k25*sqr(k46)-k35*k46*k56-k45*(k46*k26+k36*k56)-2*sqr(k56)*k26;
00255   koeficienty[4] = -k45*k46*k56-pow(k56,3);
00256   
00257   // normalization of coefficients
00258   int moc=(int(log10(fabs(koeficienty[0])))+int(log10(fabs(koeficienty[4]))))/2;
00259   
00260   koeficienty[0]=koeficienty[0]/TMath::Power(10,moc);
00261   koeficienty[1]=koeficienty[1]/TMath::Power(10,moc);
00262   koeficienty[2]=koeficienty[2]/TMath::Power(10,moc);
00263   koeficienty[3]=koeficienty[3]/TMath::Power(10,moc);
00264   koeficienty[4]=koeficienty[4]/TMath::Power(10,moc);
00265 }
00266 
00267 void TtFullLepKinSolver::TopRec(const TLorentzVector al, 
00268                                 const TLorentzVector l,
00269                                 const TLorentzVector b_al,
00270                                 const TLorentzVector b_l, 
00271                                 double sol) {
00272 
00273   TVector3 t_ttboost;
00274   TLorentzVector aux;
00275   double pxp, pyp, pzp, pup, pvp, pwp;
00276   
00277     
00278   pxp = sol*(m3*n3/k51);   
00279   pyp = -(m3*n3/k61)*(k56*pow(sol,2) + k26*sol + k16)/(k36 + k46*sol);
00280   pzp = -1/n3*(n1*pxp + n2*pyp - F);
00281   pwp = 1/m3*(m1*pxp + m2*pyp + pom);
00282   pup = C - pxp;
00283   pvp = D - pyp;
00284      
00285   LV_n_.SetXYZM(pxp, pyp, pzp, 0.0);
00286   LV_n.SetXYZM(pup, pvp, pwp, 0.0);
00287       
00288   
00289   LV_t_ = b_l + l + LV_n_;
00290   LV_t = b_al + al + LV_n;  
00291   
00292   
00293   aux = (LV_t_ + LV_t);
00294   t_ttboost = -aux.BoostVector();
00295   LV_tt_t_ = LV_t_;
00296   LV_tt_t = LV_t;
00297   LV_tt_t_.Boost(t_ttboost);
00298   LV_tt_t.Boost(t_ttboost); 
00299 }
00300 
00301 double TtFullLepKinSolver::WeightSolfromMC() {
00302 
00303   double weight = 1;
00304   weight = ((LV_n.E() > genLV_n.E())? genLV_n.E()/LV_n.E(): LV_n.E()/genLV_n.E())
00305            *((LV_n_.E() > genLV_n_.E())? genLV_n_.E()/LV_n_.E(): LV_n_.E()/genLV_n_.E());
00306   return weight;
00307 
00308 }
00309 
00310 double TtFullLepKinSolver::WeightSolfromShape() {
00311   
00312   // Use the parametrized event shape to obtain the solution weight.
00313   return EventShape_->Eval(LV_n.E(),LV_n_.E());
00314 
00315 }
00316                      
00317 int TtFullLepKinSolver::quartic(double *koeficienty, double* koreny) {
00318   int i, nreal;
00319   double w, b0, b1, b2;
00320   double c[4];
00321   double d0, d1, h, t, z;
00322   double *px;
00323 
00324   
00325   if (koeficienty[4]==0.0) 
00326     return cubic(koeficienty, koreny);
00327   /* quartic problem? */
00328   w = koeficienty[3]/(4*koeficienty[4]);
00329   /* offset */
00330   b2 = -6*sqr(w) + koeficienty[2]/koeficienty[4];
00331   /* koeficienty. of shifted polynomial */
00332   b1 = (8*sqr(w) - 2*koeficienty[2]/koeficienty[4])*w + koeficienty[1]/koeficienty[4];
00333   b0 = ((-3*sqr(w) + koeficienty[2]/koeficienty[4])*w - koeficienty[1]/koeficienty[4])*w + koeficienty[0]/koeficienty[4];
00334 
00335   c[3] = 1.0;
00336   /* cubic resolvent */
00337   c[2] = b2;
00338   c[1] = -4*b0;
00339   c[0] = sqr(b1) - 4*b0*b2;
00340   
00341   
00342   i = cubic(c, koreny);
00343   z = koreny[0];
00344   //double z1=1.0,z2=2.0,z3=3.0;
00345   //TMath::RootsCubic(c,z1,z2,z3);
00346   //if (z2 !=0) z = z2;
00347   //if (z1 !=0) z = z1;
00348   /* only lowermost root needed */
00349 
00350   nreal = 0;
00351   px = koreny;
00352   t = sqrt(0.25*sqr(z) - b0);
00353   for (i=-1; i<=1; i+=2) {
00354     d0 = -0.5*z + i*t;
00355     /* coeffs. of quadratic factor */
00356     d1 = (t!=0.0)? -i*0.5*b1/t : i*sqrt(-z - b2);
00357     h = 0.25*sqr(d1) - d0;
00358     if (h>=0.0) {
00359       h = sqrt(h);
00360       nreal += 2;
00361       *px++ = -0.5*d1 - h - w;
00362       *px++ = -0.5*d1 + h - w;
00363     }
00364   }
00365 
00366   if (nreal==4) {
00367     /* sort results */
00368 //    if (koreny[2]<koreny[0]) SWAP(koreny[0], koreny[2]);
00369 //    if (koreny[3]<koreny[1]) SWAP(koreny[1], koreny[3]);
00370 //    if (koreny[1]<koreny[0]) SWAP(koreny[0], koreny[1]);
00371 //    if (koreny[3]<koreny[2]) SWAP(koreny[2], koreny[3]);
00372 //    if (koreny[2]<koreny[1]) SWAP(koreny[1], koreny[2]);
00373   }
00374   return nreal;
00375 
00376 }
00377 
00378 int TtFullLepKinSolver::cubic(double *coeffs, double* koreny) {
00379   int i, nreal;
00380   double w, p, q, dis, h, phi;
00381   
00382   if (coeffs[3]!=0.0) {
00383     /* cubic problem? */
00384     w = coeffs[2]/(3*coeffs[3]);
00385     p = sqr(coeffs[1]/(3*coeffs[3])-sqr(w))*(coeffs[1]/(3*coeffs[3])-sqr(w));
00386     q = -0.5*(2*sqr(w)*w-(coeffs[1]*w-coeffs[0])/coeffs[3]);
00387     dis = sqr(q)+p;
00388     /* discriminant */
00389     if (dis<0.0) {
00390       /* 3 real solutions */
00391       h = q/sqrt(-p);
00392       if (h>1.0) h = 1.0;
00393       /* confine the argument of */
00394       if (h<-1.0) h = -1.0;
00395       /* acos to [-1;+1] */
00396       phi = acos(h);
00397       p = 2*TMath::Power(-p, 1.0/6.0);
00398       for (i=0; i<3; i++) 
00399         koreny[i] = p*cos((phi+2*i*TMath::Pi())/3.0) - w;
00400       if (koreny[1]<koreny[0]) SWAP(koreny[0], koreny[1]);
00401       /* sort results */
00402       if (koreny[2]<koreny[1]) SWAP(koreny[1], koreny[2]);
00403       if (koreny[1]<koreny[0]) SWAP(koreny[0], koreny[1]);
00404       nreal = 3;
00405     }
00406     else {
00407       /* only one real solution */
00408       dis = sqrt(dis);
00409       h = TMath::Power(fabs(q+dis), 1.0/3.0);
00410       p = TMath::Power(fabs(q-dis), 1.0/3.0);
00411       koreny[0] = ((q+dis>0.0)? h : -h) + ((q-dis>0.0)? p : -p) -  w;
00412       nreal = 1;
00413     }
00414 
00415     /* Perform one step of a Newton iteration in order to minimize
00416        round-off errors */
00417     for (i=0; i<nreal; i++) {
00418       h = coeffs[1] + koreny[i] * (2 * coeffs[2] + 3 * koreny[i] * coeffs[3]);
00419       if (h != 0.0)
00420         koreny[i] -= (coeffs[0] + koreny[i] * (coeffs[1] + koreny[i] * (coeffs[2] + koreny[i] * coeffs[3])))/h;
00421     }
00422   }
00423 
00424   else if (coeffs[2]!=0.0) {
00425     /* quadratic problem? */
00426     p = 0.5*coeffs[1]/coeffs[2];
00427     dis = sqr(p) - coeffs[0]/coeffs[2];
00428     if (dis>=0.0) {
00429       /* two real solutions */
00430       dis = sqrt(dis);
00431       koreny[0] = -p - dis;
00432       koreny[1] = -p + dis;
00433       nreal = 2;
00434     }
00435     else
00436       /* no real solution */
00437       nreal = 0;
00438   }
00439 
00440   else if (coeffs[1]!=0.0) {
00441     /* linear problem? */
00442     koreny[0] = -coeffs[0]/coeffs[1];
00443     nreal = 1;
00444   }
00445 
00446   else
00447     /* no equation */
00448     nreal = 0;
00449   
00450   return nreal;
00451 }
00452 
00453 
00454 void  TtFullLepKinSolver::SWAP(double& realone, double& realtwo) {
00455   if (realtwo < realone) {
00456     double aux = realtwo;
00457     realtwo = realone;
00458     realone = aux;
00459   }
00460 }
/data/refman/pasoursint/CMSSW_5_3_1/src/TopQuarkAnalysis/TopKinFitter/src/TtFullLepKinSolver.cc
