BOSS6.5.5: KalmanKinematicFit Class Reference

00515                              {
00516     bool okfit = false;
00517     TStopwatch timer;
00518     m_nktrk = numberWTrack();
00519     m_p0 = HepVector(numberone(), 0);
00520     m_p  = HepVector(numberone(), 0);
00521     m_q0 = HepVector(numbertwo(), 0);
00522     m_q  = HepVector(numbertwo(), 0);
00523     m_C0 = HepSymMatrix(numberone(), 0);
00524     m_C  = HepSymMatrix(numberone(), 0);
00525     m_D0inv = HepSymMatrix(numbertwo(), 0);
00526     m_D  = HepSymMatrix(numbertwo(), 0);
00527 
00528     HepVector  m_p_tmp = HepVector(numberone(), 0);
00529     HepVector  m_q_tmp = HepVector(numbertwo(), 0);
00530     HepMatrix  m_A_tmp;
00531     HepSymMatrix m_W_tmp;
00532     HepMatrix  m_KQ_tmp;
00533     HepSymMatrix m_Dinv_tmp;
00534     for(int i = 0; i < numberWTrack(); i++) {
00535         setWTrackInfit(i, wTrackOrigin(i));
00536         int pa = mappositionA()[i] + 1;
00537         int pb = mappositionB()[i] + 1;
00538         int ptype = mapkinematic()[i];
00539         switch(ptype) {
00540             case 0 : {
00541                          HepVector w1(4, 0);
00542                          HepSymMatrix Ew1(4, 0);
00543                          for(int j = 0; j < 3; j++) {
00544                              w1[j] = wTrackOrigin(i).w()[j];
00545                          }
00546                          w1[3] = wTrackOrigin(i).mass();
00547 
00548                          for(int m = 0; m < 3; m++) {
00549                              for(int n = 0; n < 3; n++) {
00550                                  Ew1[m][n] = wTrackOrigin(i).Ew()[m][n];
00551                              }
00552                          }
00553                          setPOrigin(pa, w1);
00554                          setPInfit(pa, w1);
00555                          setCOrigin(pa, Ew1);
00556                          break;
00557                      }
00558             case 1 : {
00559                          setPOrigin(pa, (wTrackOrigin(i).plmp()).sub(1, 4));
00560                          setPInfit(pa, (wTrackOrigin(i).plmp()).sub(1, 4));
00561                          setCOrigin(pa, (wTrackOrigin(i).Vplm()).sub(1, 4));
00562                          break;
00563                      }
00564             case 2 : {
00565                          setQOrigin(pb, (wTrackOrigin(i).w()).sub(1, NTRKPAR));
00566                          setQInfit(pb, (wTrackOrigin(i).w()).sub(1, NTRKPAR));
00567                          HepSymMatrix Dinv(4,0);
00568                          setDOriginInv(pb, Dinv);
00569                          break;
00570                      }
00571             case 3 : {
00572                          setPOrigin(pa, (wTrackOrigin(i).plmp()).sub(3, 3));
00573                          setPInfit(pa, (wTrackOrigin(i).plmp()).sub(3, 3));
00574                          setCOrigin(pa, (wTrackOrigin(i).Vplm()).sub(3, 3));
00575                          setQOrigin(pb, (wTrackOrigin(i).plmp()).sub(1, 2));
00576                          setQInfit(pb, (wTrackOrigin(i).plmp()).sub(1, 2));
00577                          setQOrigin(pb+2, (wTrackOrigin(i).plmp()).sub(4, 4));
00578                          setQInfit(pb+2, (wTrackOrigin(i).plmp()).sub(4, 4));
00579                          HepSymMatrix Dinv(3,0);
00580                          setDOriginInv(pb, Dinv);
00581                          break;
00582                      }
00583             case 4 : {
00584                          setPOrigin(pa, (wTrackOrigin(i).plmp()).sub(1, 2));
00585                          setPInfit(pa, (wTrackOrigin(i).plmp()).sub(1, 2));
00586                          setCOrigin(pa, (wTrackOrigin(i).Vplm()).sub(1, 2));
00587                          setQOrigin(pb, (wTrackOrigin(i).plmp()).sub(3, 4));
00588                          setQInfit(pb, (wTrackOrigin(i).plmp()).sub(3, 4));
00589                          HepSymMatrix Dinv(2,0);
00590                          setDOriginInv(pb, Dinv);
00591                          break;
00592                      }
00593             case 5 : {
00594                          setPOrigin(pa, (wTrackOrigin(i).plmp()).sub(1, 3));
00595                          setPInfit(pa, (wTrackOrigin(i).plmp()).sub(1, 3));
00596                          setCOrigin(pa, (wTrackOrigin(i).Vplm()).sub(1, 3));
00597                          setQOrigin(pb, (wTrackOrigin(i).plmp()).sub(4, 4));
00598                          setQInfit(pb, (wTrackOrigin(i).plmp()).sub(4, 4));
00599                          HepSymMatrix Dinv(1,0);
00600                          setDOriginInv(pb, Dinv);
00601                          break;
00602                      }
00603             case 6 : {
00604                          setPOrigin(pa, (wTrackOrigin(i).w()).sub(5, 7));
00605                          setPOrigin(pa+3, (wTrackOrigin(i).w()).sub(4, 4));
00606                          setPInfit(pa, (wTrackOrigin(i).w()).sub(5, 7));
00607                          setPInfit(pa+3, (wTrackOrigin(i).w()).sub(4, 4));
00608                          setCOrigin(pa, (wTrackOrigin(i).Ew()).sub(5,7));
00609                          setCOrigin(pa+3, (wTrackOrigin(i).Ew()).sub(4,4));
00610                          HepVector beam(3,0);
00611                          beam[0] = getBeamPosition().x();
00612                          beam[1] = getBeamPosition().y();
00613                          beam[2] = getBeamPosition().z();
00614                          setQOrigin(pb, beam);
00615                          setQInfit(pb, beam);
00616                          HepSymMatrix Dinv(3, 0);
00617                          int ifail;
00618                          Dinv = getVBeamPosition().inverse(ifail);
00619                          setDOriginInv(pb, Dinv);
00620                          break;
00621                      }
00622             case 7 : {
00623                          HepVector w1(4, 0);
00624                          HepSymMatrix Ew1(4, 0);
00625                          for(int j = 0; j < 3; j++) {
00626                              w1[j] = wTrackOrigin(i).w()[j];
00627                          }
00628                          w1[3] = wTrackOrigin(i).mass();
00629 
00630                          for(int m = 0; m < 3; m++) {
00631                              for(int n = 0; n < 3; n++) {
00632                                  Ew1[m][n] = wTrackOrigin(i).Ew()[m][n];
00633                              }
00634                          }
00635                          setPOrigin(pa, w1);
00636                          setPInfit(pa, w1);
00637                          setCOrigin(pa, Ew1);
00638                          break;
00639                      }
00640 
00641         }
00642     }
00643 
00644 
00645     //
00646     // iteration
00647     //
00648 
00649     std::vector<double> chisq;
00650     chisq.clear();
00651     int nc = 0;
00652     for(int i = 0; i < m_kc.size(); i++)
00653         nc += m_kc[i].nc();
00654 
00655     m_A = HepMatrix(nc, numberone(), 0);
00656     m_AT = HepMatrix(numberone(), nc, 0);
00657     m_B = HepMatrix(nc, numbertwo(), 0);
00658     m_BT = HepMatrix(numbertwo(), nc, 0);
00659     m_G = HepVector(nc, 0);
00660     m_Vm = HepSymMatrix(nc, 0);
00661     int num1 = 0; 
00662     for(unsigned int i = 0; i < m_kc.size(); i++) {
00663         KinematicConstraints kc = m_kc[i];
00664         m_Vm.sub(num1+1,kc.Vmeasure());
00665         num1 = num1 + kc.nc();
00666     }
00667 
00668     double tmp_chisq = 999;
00669     bool flag_break = 0;
00670     for(int it = 0; it < m_niter; it++) {
00671         timer.Start();
00672         m_nc = 0;
00673         for(unsigned int i = 0; i < m_kc.size(); i++) {
00674             KinematicConstraints kc = m_kc[i];
00675             updateConstraints(kc);
00676         }
00677         timer.Stop();
00678         m_cpu[0] += timer.CpuTime();
00679         fit();
00680         //              
00681         //reset origin parameters for virtual particle
00682         //
00683         //                if(it == 0){
00684         //                    for(int i = 0; i < numberWTrack(); i++) {
00685         //                      //  setWTrackInfit(i, wTrackOrigin(i));
00686         //                        int pa = mappositionA()[i] + 1;
00687         //                        int pb = mappositionB()[i] + 1;
00688         //                        int ptype = mapkinematic()[i];
00689         //                        switch(ptype) {
00690         //                           case 3 : {
00691         //                                         setPOrigin(pa, m_p.sub(pa, pa));
00692         //                                         setPInfit(pa,  m_p.sub(pa, pa));
00693         //                                         setQOrigin(pb, m_q.sub(pb, pb+2));
00694         //                                         setQInfit(pb,  m_q.sub(pb, pb+2));
00695         //                                         break;
00696         //                                    }
00697         //                           case 4 : {
00698         //                                        setPOrigin(pa, (wTrackOrigin(i).plmp()).sub(1, 2));
00699         //                                        setPInfit(pa, (wTrackOrigin(i).plmp()).sub(1, 2));
00700         //                                        setCOrigin(pa, (wTrackOrigin(i).Vplm()).sub(1, 2));
00701         //                                        setQOrigin(pb, (wTrackOrigin(i).plmp()).sub(3, 4));
00702         //                                        setQInfit(pb, (wTrackOrigin(i).plmp()).sub(3, 4));
00703         //                                        HepSymMatrix Dinv(2,0);
00704         //                                        setDOriginInv(pb, Dinv);
00705         //                                        break;
00706         //                                    }
00707         //                           case 5 : {
00708         //                                        setPOrigin(pa, (wTrackOrigin(i).plmp()).sub(1, 3));
00709         //                                        setPInfit(pa, (wTrackOrigin(i).plmp()).sub(1, 3));
00710         //                                        setCOrigin(pa, (wTrackOrigin(i).Vplm()).sub(1, 3));
00711         //                                        setQOrigin(pb, (wTrackOrigin(i).plmp()).sub(4, 4));
00712         //                                        setQInfit(pb, (wTrackOrigin(i).plmp()).sub(4, 4));
00713         //                                        HepSymMatrix Dinv(1,0);
00714         //                                        setDOriginInv(pb, Dinv);
00715         //                                        break;
00716         //                                    }
00717         //                           case 6 : {
00718         //                                        setPOrigin(pa, (wTrackOrigin(i).w()).sub(5, 7));
00719         //                                        setPOrigin(pa+3, (wTrackOrigin(i).w()).sub(4, 4));
00720         //                                        setPInfit(pa, (wTrackOrigin(i).w()).sub(5, 7));
00721         //                                        setPInfit(pa+3, (wTrackOrigin(i).w()).sub(4, 4));
00722         //                                        setCOrigin(pa, (wTrackOrigin(i).Ew()).sub(5,7));
00723         //                                        setCOrigin(pa+3, (wTrackOrigin(i).Ew()).sub(4,4));
00724         //                                        HepVector beam(3,0);
00725         //                                        beam[0] = getBeamPosition().x();
00726         //                                        beam[1] = getBeamPosition().y();
00727         //                                        beam[2] = getBeamPosition().z();
00728         //                                        setQOrigin(pb, beam);
00729         //                                        setQInfit(pb, beam);
00730         //                                        HepSymMatrix Dinv(3, 0);
00731         //                                        int ifail;
00732         //                                        Dinv = getVBeamPosition().inverse(ifail);
00733         //                                        setDOriginInv(pb, Dinv);
00734         //                                        break;
00735         //                                    }
00736         //
00737         //                        }
00738         //                    }
00739         //
00740         //                }
00741         //
00742         //
00743         //===================reset over=============================
00744         //
00745         chisq.push_back(m_chi);
00746         if(it > 0) {
00747             double delchi = chisq[it]- chisq[it-1];  
00748             if(fabs(delchi) < m_chiter){
00749                 flag_break =1; 
00750                 break;
00751             }
00752         }
00753     }
00754     if(!flag_break) {
00755         return okfit;
00756     }
00757     if(m_chi > m_chicut) return okfit;
00758     timer.Start();
00759     timer.Stop();
00760     m_cpu[5] += timer.CpuTime();
00761     okfit = true;
00762     return okfit;
00763 
00764 }

01450                                                                   {
01451     KinematicConstraints kc = k;
01452 
01453 
01454     int type = kc.Type();
01455     switch(type) {
01456         case  Resonance: {
01457                              //
01458                              //  E^2 - px^2 - py^2 - pz^2 = mres^2
01459                              //
01460                              double mres = kc.mres();
01461                              HepLorentzVector pmis;
01462                              for(unsigned int j = 0; j< (kc.Ltrk()).size(); j++){
01463                                  int n = (kc.Ltrk())[j];
01464                                  pmis = pmis + p4Infit(n);
01465                              }
01466                              for(unsigned int j = 0; j< (kc.Ltrk()).size(); j++){
01467                                  int n = (kc.Ltrk())[j];
01468                                  double lambda = p4Infit(n).pz()/p4Infit(n).perp();
01469                                  double a1 = 1 + lambda*lambda;
01470                                  int pa = mappositionA()[n] + 1;
01471                                  int pb = mappositionB()[n] + 1;
01472                                  int ptype = mapkinematic()[n];
01473                                  switch(ptype) {
01474                                      case 0 :{
01475                                                  HepMatrix ta(1, NTRKPAR, 0);
01476                                                  ta[0][0] = -2 * pmis.px() + 2 * pmis.e() * p4Infit(n).px() / p4Infit(n).e();
01477                                                  ta[0][1] = -2 * pmis.py() + 2 * pmis.e() * p4Infit(n).py() / p4Infit(n).e();
01478                                                  ta[0][2] = -2 * pmis.pz() + 2 * pmis.e() * p4Infit(n).pz() / p4Infit(n).e();
01479                                                  ta[0][3] = 2 * pmis.e() * p4Infit(n).m() / p4Infit(n).e();
01480                                                  setA(m_nc,pa,ta);
01481                                                  setAT(pa, m_nc, ta.T());
01482                                                  break;
01483                                              }
01484                                      case 1 :  {
01485                                                    HepMatrix ta(1, NTRKPAR, 0);
01486                                                    double a1 = lambda * lambda + 1;
01487                                                    ta[0][0] = 2 * pmis.px() * p4Infit(n).py() - 2 * pmis.py() * p4Infit(n).px();
01488                                                    ta[0][1] = 2 * pmis.px() * p4Infit(n).px() * lambda/sqrt(a1) + 2 * pmis.py() * (p4Infit(n)).py() * lambda/sqrt(a1) - 2 * pmis.pz() * (p4Infit(n)).pz() /(sqrt(a1) * lambda);
01489                                                    ta[0][2] =  2 * pmis.px() * (p4Infit(n)).px() * p4Infit(n).m() /((p4Infit(n)).rho() * p4Infit(n).rho()) + 2 * pmis.py() * (p4Infit(n)).py() * p4Infit(n).m()/((p4Infit(n)).rho() * p4Infit(n).rho())+ 2 * pmis.pz() * (p4Infit(n)).pz() * p4Infit(n).m()/((p4Infit(n)).rho() * p4Infit(n).rho());
01490                                                    ta[0][3] = 2 * pmis.e()  - 2 * pmis.px() * (p4Infit(n)).px() * p4Infit(n).e() /((p4Infit(n)).rho() * p4Infit(n).rho())- 2 * pmis.py() * (p4Infit(n)).py() * p4Infit(n).e()/((p4Infit(n)).rho() * p4Infit(n).rho())- 2 * pmis.pz() * (p4Infit(n)).pz() * p4Infit(n).e()/((p4Infit(n)).rho() * p4Infit(n).rho());
01491                                                    setA(m_nc,pa,ta);
01492                                                    setAT(pa, m_nc, ta.T());
01493                                                    break;
01494                                                }
01495                                      case 2 :  {
01496                                                    HepMatrix tb(1, 4, 0);
01497                                                    tb[0][0] = -2 * pmis.px();
01498                                                    tb[0][1] = -2 * pmis.py();
01499                                                    tb[0][2] = -2 * pmis.pz();
01500                                                    tb[0][3] = 2 * pmis.e();
01501                                                    setB(m_nc,pb,tb);
01502                                                    setBT(pb,m_nc,tb.T());
01503                                                    break;
01504                                                }
01505                                      case 3 :  {
01506                                                    HepMatrix ta(1, 1, 0);
01507                                                    double a1 = lambda * lambda + 1;
01508                                                    ta[0][0] =  2 * pmis.px() * (p4Infit(n)).px() * p4Infit(n).m() /((p4Infit(n)).rho() * p4Infit(n).rho()) + 2 * pmis.py() * (p4Infit(n)).py() * p4Infit(n).m()/((p4Infit(n)).rho() * p4Infit(n).rho())+ 2 * pmis.pz() * (p4Infit(n)).pz() * p4Infit(n).m()/((p4Infit(n)).rho() * p4Infit(n).rho());
01509                                                    setA(m_nc,pa,ta);
01510                                                    setAT(pa, m_nc, ta.T());
01511                                                    HepMatrix tb(1, 3, 0);
01512                                                    tb[0][0] = 2 * pmis.px() * p4Infit(n).py() - 2 * pmis.py() * p4Infit(n).px();
01513                                                    tb[0][1] = 2 * pmis.px() * p4Infit(n).px() * lambda/sqrt(a1) + 2 * pmis.py() * (p4Infit(n)).py() * lambda/sqrt(a1) - 2 * pmis.pz() * (p4Infit(n)).pz() /(sqrt(a1) * lambda);
01514                                                    tb[0][2] = 2 * pmis.e()  - 2 * pmis.px() * (p4Infit(n)).px() * p4Infit(n).e() /((p4Infit(n)).rho() * p4Infit(n).rho())- 2 * pmis.py() * (p4Infit(n)).py() * p4Infit(n).e()/((p4Infit(n)).rho() * p4Infit(n).rho())- 2 * pmis.pz() * (p4Infit(n)).pz() * p4Infit(n).e()/((p4Infit(n)).rho() * p4Infit(n).rho());
01515                                                    setB(m_nc,pb,tb);
01516                                                    setBT(pb,m_nc,tb.T());
01517                                                    break;
01518                                                }
01519                                      case 4 :  {
01520                                                    HepMatrix ta(1, 2, 0);
01521                                                    double a1 = lambda * lambda + 1;
01522                                                    ta[0][0] = 2 * pmis.px() * p4Infit(n).py() - 2 * pmis.py() * p4Infit(n).px();
01523                                                    ta[0][1] = 2 * pmis.px() * p4Infit(n).px() * lambda/sqrt(a1) + 2 * pmis.py() * (p4Infit(n)).py() * lambda/sqrt(a1) - 2 * pmis.pz() * (p4Infit(n)).pz() /(sqrt(a1) * lambda);
01524                                                    setA(m_nc,pa,ta);
01525                                                    setAT(pa, m_nc, ta.T());
01526 
01527                                                    HepMatrix tb(1, 2, 0);
01528                                                    tb[0][0] =  2 * pmis.px() * (p4Infit(n)).px() * p4Infit(n).m() /((p4Infit(n)).rho() * p4Infit(n).rho()) + 2 * pmis.py() * (p4Infit(n)).py() * p4Infit(n).m()/((p4Infit(n)).rho() * p4Infit(n).rho())+ 2 * pmis.pz() * (p4Infit(n)).pz() * p4Infit(n).m()/((p4Infit(n)).rho() * p4Infit(n).rho());
01529                                                    tb[0][1] = 2 * pmis.e()  - 2 * pmis.px() * (p4Infit(n)).px() * p4Infit(n).e() /((p4Infit(n)).rho() * p4Infit(n).rho())- 2 * pmis.py() * (p4Infit(n)).py() * p4Infit(n).e()/((p4Infit(n)).rho() * p4Infit(n).rho())- 2 * pmis.pz() * (p4Infit(n)).pz() * p4Infit(n).e()/((p4Infit(n)).rho() * p4Infit(n).rho());
01530                                                    setB(m_nc,pb,tb);
01531                                                    setBT(pb, m_nc, tb.T());
01532                                                    break;
01533                                                }
01534                                      case 5 :  {
01535                                                    HepMatrix ta(1, 3, 0);
01536                                                    ta[0][0] = 2 * pmis.px() * p4Infit(n).py() - 2 * pmis.py() * p4Infit(n).px();
01537                                                    ta[0][1] = 2 * pmis.px() * p4Infit(n).px() * lambda/sqrt(a1)+ 2 * pmis.py() * (p4Infit(n)).py() * lambda/sqrt(a1) - 2 * pmis.pz() * (p4Infit(n)).pz() /(sqrt(a1) * lambda);
01538                                                    ta[0][2] = 2 * pmis.e() * (p4Infit(n)).m()/(p4Infit(n)).e();
01539                                                    setA(m_nc,pa,ta);
01540                                                    setAT(pa, m_nc, ta.T());
01541                                                    HepMatrix tb(1, 1, 0);
01542                                                    tb[0][0] = 2 * pmis.e() * (p4Infit(n)).rho()/(p4Infit(n)).e() - 2 * pmis.px() * (p4Infit(n)).px()/(p4Infit(n)).rho() - 2 * pmis.py() * (p4Infit(n)).py()/(p4Infit(n)).rho() - 2 * pmis.pz() * (p4Infit(n)).pz()/(p4Infit(n)).rho();
01543                                                    setB(m_nc,pb,tb);
01544                                                    setBT(pb, m_nc, tb.T());
01545                                                    break;
01546                                                }       
01547                                      case 7 :{
01548                                                  HepMatrix ta(1, NTRKPAR, 0);
01549                                                  ta[0][0] = -2 * pmis.px() + 2 * pmis.e() * p4Infit(n).px() / p4Infit(n).e();
01550                                                  ta[0][1] = -2 * pmis.py() + 2 * pmis.e() * p4Infit(n).py() / p4Infit(n).e();
01551                                                  ta[0][2] = -2 * pmis.pz() + 2 * pmis.e() * p4Infit(n).pz() / p4Infit(n).e();
01552                                                  ta[0][3] = 2 * pmis.e() * p4Infit(n).m() / p4Infit(n).e();
01553                                                  setA(m_nc,pa,ta);
01554                                                  setAT(pa, m_nc, ta.T());
01555                                                  break;
01556                                              }
01557                                  }   
01558                              }
01559 
01560                              HepVector dc(1, 0);
01561                              dc[0] = pmis.m2() - mres * mres;
01562                              m_G[m_nc] = dc[0];
01563                              m_nc+=1;
01564 
01565                              break;
01566                          }
01567         case TotalEnergy: {
01568                               //
01569                               //  E - Etot = 0
01570                               //
01571                               double etot = kc.etot();
01572                               HepLorentzVector pmis;
01573                               for(unsigned int j = 0; j < (kc.Ltrk()).size(); j++){
01574                                   int n = (kc.Ltrk())[j];
01575                                   pmis = pmis + p4Infit(n);
01576                               }
01577 
01578                               for(unsigned int j = 0; j < (kc.Ltrk()).size(); j++) {
01579                                   int n = (kc.Ltrk())[j];
01580                                   int pa = mappositionA()[n] + 1;
01581                                   int pb = mappositionB()[n] + 1;
01582                                   int ptype = mapkinematic()[n];
01583                                   switch(ptype) {
01584                                       case 0 :{
01585                                                   HepMatrix ta(1, NTRKPAR, 0);
01586                                                   ta[0][0] = p4Infit(n).px()/p4Infit(n).e();
01587                                                   ta[0][1] = p4Infit(n).py()/p4Infit(n).e();
01588                                                   ta[0][2] = p4Infit(n).pz()/p4Infit(n).e();
01589                                                   ta[0][3] = p4Infit(n).m()/p4Infit(n).e();
01590                                                   setA(m_nc,pa,ta);
01591                                                   setAT(pa, m_nc, ta.T());
01592                                                   break;
01593                                               }
01594                                       case 1: {
01595                                                   HepMatrix ta(1, NTRKPAR, 0);
01596                                                   ta[0][3] = 1.0;
01597                                                   setA(m_nc,pa,ta);
01598                                                   setAT(pa, m_nc, ta.T());
01599                                                   break;
01600                                               }
01601                                       case 2: {
01602                                                   HepMatrix tb(1, 4, 0);
01603                                                   tb[0][3] = 1.0;
01604                                                   setA(m_nc,pb,tb);
01605                                                   setAT(pb, m_nc, tb.T());
01606                                                   break;
01607                                               }
01608                                       case 3: {
01609                                                   HepMatrix ta(1, 1, 0);
01610                                                   ta[0][0] = p4Infit(n).m()/p4Infit(n).e();
01611                                                   setA(m_nc,pa,ta);
01612                                                   setAT(pa, m_nc, ta.T());
01613 
01614                                                   HepMatrix tb(1, 3, 0);
01615                                                   setB(m_nc,pb,tb);
01616                                                   setBT(pb, m_nc, tb.T());
01617                                                   break;
01618                                               }
01619                                       case 4: {
01620                                                   HepMatrix ta(1, 2, 0);
01621                                                   setA(m_nc,pa,ta);
01622                                                   setAT(pa, m_nc, ta.T());
01623 
01624                                                   HepMatrix tb(1, 2, 0);
01625                                                   tb[0][0] = 0.0;
01626                                                   tb[0][1] = 1.0;
01627                                                   //   tb[0][0] = p4Infit(n).m()/p4Infit(n).e();
01628                                                   //   tb[0][1] = p4Infit(n).rho()/p4Infit(n).e();
01629                                                   setB(m_nc,pb,tb);
01630                                                   setBT(pb, m_nc, tb.T());
01631                                                   break;
01632                                               }
01633                                       case 5: {
01634                                                   HepMatrix ta(1, 3, 0);
01635                                                   ta[0][2] = p4Infit(n).m()/p4Infit(n).e();
01636                                                   setA(m_nc,pa,ta);
01637                                                   setAT(pa, m_nc, ta.T());
01638 
01639                                                   HepMatrix tb(1, 1, 0);
01640                                                   tb[0][0] = p4Infit(n).rho()/p4Infit(n).e();
01641                                                   setB(m_nc,pb,tb);
01642                                                   setBT(pb, m_nc, tb.T());
01643                                                   break;
01644                                               } 
01645                                       case 7 :{
01646                                                   HepMatrix ta(1, NTRKPAR, 0);
01647                                                   ta[0][0] = p4Infit(n).px()/p4Infit(n).e();
01648                                                   ta[0][1] = p4Infit(n).py()/p4Infit(n).e();
01649                                                   ta[0][2] = p4Infit(n).pz()/p4Infit(n).e();
01650                                                   ta[0][3] = p4Infit(n).m()/p4Infit(n).e();
01651                                                   setA(m_nc,pa,ta);
01652                                                   setAT(pa, m_nc, ta.T());
01653                                                   break;
01654                                               }
01655                                   }     
01656 
01657                               }
01658 
01659                               HepVector dc(1, 0);
01660                               dc[0] = pmis.e() - etot;
01661                               m_G[m_nc] = dc[0];
01662                               m_nc+=1;
01663                               break;
01664                           }
01665         case TotalMomentum: {
01666                                 //
01667                                 //  sqrt(px^2+py^2+pz^2) - ptot = 0
01668                                 //
01669                                 double ptot = kc.ptot();
01670                                 HepLorentzVector pmis;
01671                                 for(unsigned int j = 0; j< (kc.Ltrk()).size(); j++){
01672                                     int n = (kc.Ltrk())[j];
01673                                     pmis = pmis + p4Infit(n);
01674                                 }
01675 
01676                                 for(unsigned int j = 0; j < (kc.Ltrk()).size(); j++) {
01677                                     int n = (kc.Ltrk())[j];
01678                                     int pa = mappositionA()[n] + 1;
01679                                     int pb = mappositionB()[n] + 1;
01680                                     int ptype = mapkinematic()[n];
01681                                     double lambda = p4Infit(n).pz()/p4Infit(n).perp();
01682                                     switch(ptype) {
01683                                         case 0 : {
01684                                                      HepMatrix ta(1, NTRKPAR, 0);
01685                                                      ta[0][0] = pmis.px()/pmis.rho();
01686                                                      ta[0][1] = pmis.py()/pmis.rho();
01687                                                      ta[0][2] = pmis.pz()/pmis.rho();
01688                                                      setA(m_nc,pa,ta);
01689                                                      setAT(pa, m_nc, ta.T());
01690                                                      break;
01691                                                  }
01692                                         case 1 : {
01693                                                      HepMatrix ta(1, NTRKPAR, 0);
01694                                                      ta[0][0] = - (pmis.px()/pmis.rho()) * p4Infit(n).py() + (pmis.px()/pmis.rho())*p4Infit(n).px();
01695                                                      ta[0][1] = - (pmis.px()/pmis.rho()) * p4Infit(n).px() * (lambda/(1 + lambda*lambda)) - (pmis.py()/pmis.rho()) * p4Infit(n).py() * (lambda/(1 + lambda*lambda)) + (pmis.pz()/pmis.rho()) * p4Infit(n).pz() * (1/(lambda * (1 + lambda*lambda)));
01696                                                      ta[0][2] =  -((pmis.px()/pmis.rho()) * p4Infit(n).m()/(p4Infit(n).rho()*p4Infit(n).rho())) * (p4Infit(n).px() + p4Infit(n).py() +p4Infit(n).pz());
01697                                                      ta[0][3] = ((pmis.px()/pmis.rho()) * p4Infit(n).e()/(p4Infit(n).rho()*p4Infit(n).rho())) * (p4Infit(n).px() + p4Infit(n).py() +
01698                                                              p4Infit(n).pz()); 
01699 
01700                                                      setA(m_nc,pa,ta);
01701                                                      setAT(pa, m_nc, ta.T());
01702                                                      break;
01703                                                  }
01704                                         case 2 : {
01705                                                      HepMatrix tb(1, 4, 0);
01706                                                      tb[0][0] = pmis.px()/pmis.rho();
01707                                                      tb[0][1] = pmis.py()/pmis.rho();
01708                                                      tb[0][2] = pmis.pz()/pmis.rho();
01709                                                      setB(m_nc,pb,tb);
01710                                                      setBT(pb, m_nc, tb.T());
01711                                                      break;
01712                                                  }
01713                                         case 3 : {
01714                                                      HepMatrix ta(1, 1, 0); 
01715                                                      setA(m_nc,pa,ta);
01716                                                      setAT(pa, m_nc, ta.T());
01717                                                      HepMatrix tb(1, 3, 0);
01718                                                      tb[0][0] = pmis.px()/pmis.rho();
01719                                                      tb[0][1] = pmis.py()/pmis.rho();
01720                                                      tb[0][2] = pmis.pz()/pmis.rho();
01721                                                      setB(m_nc,pb,tb);
01722                                                      setBT(pb, m_nc, tb.T());
01723                                                      break;
01724                                                  }
01725                                         case 4 : {
01726                                                      HepMatrix ta(1, 2, 0);
01727                                                      ta[0][0] = - (pmis.px()/pmis.rho()) * p4Infit(n).py() + (pmis.px()/pmis.rho())*p4Infit(n).px();
01728                                                      ta[0][1] = - (pmis.px()/pmis.rho()) * p4Infit(n).px() * (lambda/(1 + lambda*lambda)) - (pmis.py()/pmis.rho()) * p4Infit(n).py() * (lambda/(1 + lambda*lambda)) + (pmis.pz()/pmis.rho()) * p4Infit(n).pz() * (1/(lambda * (1 + lambda*lambda)));
01729                                                      setA(m_nc,pa,ta);
01730                                                      setAT(pa, m_nc, ta.T());
01731 
01732                                                      HepMatrix tb(1, 2, 0);
01733                                                      tb[0][0] =  -((pmis.px()/pmis.rho()) * p4Infit(n).m()/(p4Infit(n).rho()*p4Infit(n).rho())) * (p4Infit(n).px() + p4Infit(n).py() +p4Infit(n).pz());
01734                                                      tb[0][1] = ((pmis.px()/pmis.rho()) * p4Infit(n).e()/(p4Infit(n).rho()*p4Infit(n).rho())) * (p4Infit(n).px() + p4Infit(n).py() +
01735                                                              p4Infit(n).pz());
01736                                                      //             tb[0][1] =  (pmis.px()/pmis.rho()) * (p4Infit(n).px()/p4Infit(n).rho()) + (pmis.py()/pmis.rho()) * (p4Infit(n).py()/p4Infit(n).rho()) + (pmis.pz()/pmis.rho()) * (p4Infit(n).pz()/p4Infit(n).rho());
01737                                                      setB(m_nc,pb,tb);
01738                                                      setBT(pb, m_nc, tb.T());                
01739                                                      break;
01740                                                  }
01741                                         case 5 : {
01742                                                      HepMatrix ta(1, 3, 0);
01743                                                      ta[0][0] = - (pmis.px()/pmis.rho()) * p4Infit(n).py() + (pmis.px()/pmis.rho())*p4Infit(n).px();
01744                                                      ta[0][1] = - (pmis.px()/pmis.rho()) * p4Infit(n).px() * (lambda/(1 + lambda*lambda)) - (pmis.py()/pmis.rho()) * p4Infit(n).py() * (lambda/(1 + lambda*lambda)) + (pmis.pz()/pmis.rho()) * p4Infit(n).pz() * (1/(lambda * (1 + lambda*lambda)));   
01745                                                      setA(m_nc,pa,ta);
01746                                                      setAT(pa, m_nc, ta.T());
01747 
01748                                                      HepMatrix tb(1, 1, 0);
01749                                                      tb[0][0] =  (pmis.px()/pmis.rho()) * (p4Infit(n).px()/p4Infit(n).rho()) + (pmis.py()/pmis.rho()) * (p4Infit(n).py()/p4Infit(n).rho()) + (pmis.pz()/pmis.rho()) * (p4Infit(n).pz()/p4Infit(n).rho());
01750                                                      setB(m_nc,pb,tb);
01751                                                      setBT(pb, m_nc, tb.T()); 
01752                                                      break;
01753                                                  }
01754                                         case 7 : {
01755                                                      HepMatrix ta(1, NTRKPAR, 0);
01756                                                      ta[0][0] = pmis.px()/pmis.rho();
01757                                                      ta[0][1] = pmis.py()/pmis.rho();
01758                                                      ta[0][2] = pmis.pz()/pmis.rho();
01759                                                      setA(m_nc,pa,ta);
01760                                                      setAT(pa, m_nc, ta.T());
01761                                                      break;
01762                                                  }
01763                                     }    
01764                                 }
01765 
01766 
01767                                 HepVector dc(1, 0);
01768                                 dc[0] = pmis.rho() - ptot;
01769                                 m_G[m_nc] = dc[0];
01770                                 m_nc+=1;     
01771                                 break;
01772                             }
01773 
01774 
01775         case ThreeMomentum: {
01776                                 //
01777                                 //  px - p3x = 0
01778                                 //  py - p3y = 0
01779                                 //  pz - p3z = 0
01780                                 //
01781                                 Hep3Vector p3 = kc.p3();
01782                                 HepLorentzVector pmis;
01783                                 for(unsigned int j = 0; j< (kc.Ltrk()).size(); j++){
01784                                     int n = (kc.Ltrk())[j];
01785                                     pmis = pmis + p4Infit(n);
01786                                 }
01787                                 for(unsigned int j = 0; j < (kc.Ltrk()).size(); j++) {
01788                                     int n = (kc.Ltrk())[j];
01789                                     int pa = mappositionA()[n] + 1;
01790                                     int pb = mappositionB()[n] + 1;
01791                                     int ptype = mapkinematic()[n];
01792                                     double lambda = p4Infit(n).pz()/p4Infit(n).perp();
01793                                     switch(ptype) {
01794                                         case 0 : {
01795                                                      HepMatrix ta(kc.nc(), NTRKPAR, 0);
01796                                                      ta[0][0] = 1.0;
01797                                                      ta[1][1] = 1.0;
01798                                                      ta[2][2] = 1.0;
01799                                                      setA(m_nc, pa, ta);
01800                                                      setAT(pa, m_nc, ta.T());
01801                                                      break;
01802                                                  }
01803                                         case 1 : {
01804                                                      HepMatrix ta(kc.nc(), NTRKPAR, 0);
01805                                                      ta[0][0] = -p4Infit(n).py();
01806                                                      ta[0][1] = -p4Infit(n).px()*(lambda/sqrt(1 + lambda*lambda));
01807                                                      ta[0][2] = -p4Infit(n).m()*p4Infit(n).px()/(p4Infit(n).rho()*p4Infit(n).rho());
01808                                                      ta[0][3] = p4Infit(n).e()*p4Infit(n).px()/(p4Infit(n).rho()*p4Infit(n).rho());
01809                                                      ta[1][0] = p4Infit(n).px();
01810                                                      ta[1][1] = -p4Infit(n).py()*(lambda/sqrt(1 + lambda*lambda));
01811                                                      ta[1][2] = -p4Infit(n).m()*p4Infit(n).py()/(p4Infit(n).rho()*p4Infit(n).rho());
01812                                                      ta[1][3] = p4Infit(n).e()*p4Infit(n).py()/(p4Infit(n).rho()*p4Infit(n).rho());
01813                                                      ta[2][0] = 0;
01814                                                      ta[2][1] = p4Infit(n).pz()*(1/(lambda * (1 + lambda*lambda)));
01815                                                      ta[2][2] = -p4Infit(n).m()*p4Infit(n).pz()/(p4Infit(n).rho()*p4Infit(n).rho());
01816                                                      ta[2][3] = p4Infit(n).e()*p4Infit(n).pz()/(p4Infit(n).rho()*p4Infit(n).rho());
01817                                                      //             ta[0][0] = 1.0;
01818                                                      //             ta[1][1] = 1.0;
01819                                                      //             ta[2][2] = 1.0;
01820                                                      setA(m_nc, pa, ta);
01821                                                      setAT(pa, m_nc, ta.T()); 
01822                                                      break;
01823                                                  }
01824                                         case 2 : {
01825                                                      HepMatrix tb(kc.nc(), 4, 0);
01826                                                      tb[0][0] = 1.0;
01827                                                      tb[1][1] = 1.0;
01828                                                      tb[2][2] = 1.0;
01829                                                      setB(m_nc, pb, tb);
01830                                                      setBT(pb, m_nc, tb.T());
01831                                                      break;
01832                                                  }
01833                                         case 3 : {
01834                                                      HepMatrix ta(kc.nc(), 1, 0);
01835                                                      setA(m_nc, pa, ta);
01836                                                      setAT(pa, m_nc, ta.T()); 
01837 
01838                                                      HepMatrix tb(kc.nc(), 3, 0);
01839                                                      tb[0][0] = 1.0;
01840                                                      tb[1][1] = 1.0;
01841                                                      tb[2][2] = 1.0;
01842                                                      setB(m_nc, pb, tb);
01843                                                      setBT(pb, m_nc, tb.T());
01844 
01845                                                      break;
01846                                                  }
01847                                         case 4 : {
01848                                                      HepMatrix ta(kc.nc(), 2, 0);
01849                                                      ta[0][0] = -p4Infit(n).py(); 
01850                                                      ta[0][1] = -p4Infit(n).px()*(lambda/sqrt(1 + lambda*lambda));
01851                                                      ta[1][0] =  p4Infit(n).px(); 
01852                                                      ta[1][1] = -p4Infit(n).py()*(lambda/sqrt(1 + lambda*lambda));
01853                                                      ta[2][0] =  0;
01854                                                      ta[2][1] =  p4Infit(n).pz()*(1/(lambda * (1 + lambda*lambda)));
01855                                                      setA(m_nc, pa, ta);
01856                                                      setAT(pa, m_nc, ta.T()); 
01857 
01858                                                      HepMatrix tb(kc.nc(), 2, 0);
01859                                                      tb[0][1] = -p4Infit(n).m()*p4Infit(n).px()/(p4Infit(n).rho()*p4Infit(n).rho());
01860                                                      tb[1][1] = -p4Infit(n).m()*p4Infit(n).py()/(p4Infit(n).rho()*p4Infit(n).rho());
01861                                                      tb[2][1] = -p4Infit(n).m()*p4Infit(n).pz()/(p4Infit(n).rho()*p4Infit(n).rho());
01862                                                      setB(m_nc, pb, tb);
01863                                                      setBT(pb, m_nc, tb.T());
01864                                                      break;
01865                                                  }
01866                                         case 5 : {
01867                                                      HepMatrix ta(kc.nc(), 3, 0);
01868                                                      ta[0][0] = -p4Infit(n).py(); 
01869                                                      ta[0][1] = -p4Infit(n).px()*(lambda/sqrt(1 + lambda*lambda));
01870                                                      ta[1][0] =  p4Infit(n).px(); 
01871                                                      ta[1][1] = -p4Infit(n).py()*(lambda/sqrt(1 + lambda*lambda));
01872                                                      ta[2][0] =  0;
01873                                                      ta[2][1] =  p4Infit(n).pz()*(1/(lambda * (1 + lambda*lambda)));
01874                                                      setA(m_nc, pa, ta);
01875                                                      setAT(pa, m_nc, ta.T()); 
01876 
01877                                                      HepMatrix tb(kc.nc(), 1, 0);
01878                                                      tb[0][0] = p4Infit(n).px()/p4Infit(n).rho();
01879                                                      tb[1][0] = p4Infit(n).py()/p4Infit(n).rho();
01880                                                      tb[2][0] = p4Infit(n).pz()/p4Infit(n).rho();
01881                                                      setB(m_nc, pb, tb);
01882                                                      setBT(pb, m_nc, tb.T());
01883                                                      break;
01884                                                  }
01885                                         case 7 : {
01886                                                      HepMatrix ta(kc.nc(), NTRKPAR, 0);
01887                                                      ta[0][0] = 1.0;
01888                                                      ta[1][1] = 1.0;
01889                                                      ta[2][2] = 1.0;
01890                                                      setA(m_nc, pa, ta);
01891                                                      setAT(pa, m_nc, ta.T());
01892                                                      break;
01893                                                  }
01894                                     }
01895                                 }    
01896                                 HepVector dc(kc.nc(), 0);
01897                                 dc[0] = pmis.px() - p3.x();
01898                                 dc[1] = pmis.py() - p3.y();
01899                                 dc[2] = pmis.pz() - p3.z();
01900                                 for(int i = 0; i < kc.nc(); i++) {
01901                                     m_G[m_nc+i] = dc[i];
01902                                 }        
01903 
01904                                 m_nc += 3; 
01905 
01906                                 break;
01907                             }
01908         case EqualMass: {
01909                             //
01910                             //  (E_1 ^2 - Px_1 ^2 - Py_1 ^2 - Pz_1 ^2) - (E_2 ^2 - Px_2 ^2 - Py_2 ^2 - Pz_2 ^2) = 0
01911                             //
01912 
01913                             int isiz = (kc.numEqual())[0];
01914                             HepLorentzVector pmis1, pmis2;
01915                             for(int n = 0; n < isiz; n++) {
01916                                 int n1 = (kc.Ltrk())[n]; 
01917                                 pmis1 = pmis1 + p4Infit(n1);
01918                             }
01919 
01920                             int jsiz = (kc.numEqual())[1];
01921                             for(int n = 0; n < jsiz; n++){
01922                                 int n2 = (kc.Ltrk())[n+isiz];
01923                                 pmis2 = pmis2 + p4Infit(n2);
01924                             }
01925 
01926                             for(int i = 0; i < isiz; i++) {
01927                                 int n1 = (kc.Ltrk())[i];
01928                                 double lambda = p4Infit(n1).pz()/p4Infit(n1).perp();
01929                                 int pa = mappositionA()[n1] + 1;
01930                                 int pb = mappositionB()[n1] + 1;
01931                                 int ptype = mapkinematic()[n1];
01932                                 switch(ptype) {
01933                                     case 0: {
01934                                                 HepMatrix ta(1, NTRKPAR, 0);
01935                                                 ta[0][0] = -2 * pmis1.px() + 2 * pmis1.e() * p4Infit(n1).px() / p4Infit(n1).e();
01936                                                 ta[0][1] = -2 * pmis1.py() + 2 * pmis1.e() * p4Infit(n1).py() / p4Infit(n1).e();
01937                                                 ta[0][2] = -2 * pmis1.pz() + 2 * pmis1.e() * p4Infit(n1).pz() / p4Infit(n1).e();
01938                                                 ta[0][3] = 2 * pmis1.e() *  p4Infit(n1).m() / p4Infit(n1).e();
01939                                                 setA(m_nc,pa,ta);
01940                                                 setAT(pa, m_nc, ta.T());
01941                                                 break;
01942                                             }
01943                                     case 1 :  {
01944                                                   HepMatrix ta(1, NTRKPAR, 0);
01945                                                   double a1 = lambda * lambda + 1;
01946                                                   ta[0][0] = 2 * pmis1.px() * p4Infit(n1).py() - 2 * pmis1.py() * p4Infit(n1).px();
01947                                                   ta[0][1] = 2 * pmis1.px() * p4Infit(n1).px() * lambda/sqrt(a1) + 2 * pmis1.py() * (p4Infit(n1)).py() * lambda/sqrt(a1) - 2 * pmis1.pz() * (p4Infit(n1)).pz() /(sqrt(a1) * lambda);
01948                                                   ta[0][2] = 2 * pmis1.px() * (p4Infit(n1)).px() * p4Infit(n1).m() /((p4Infit(n1)).rho() * p4Infit(n1).rho()) + 2 * pmis1.py() * (p4Infit(n1)).py() * p4Infit(n1).m()/((p4Infit(n1)).rho() * p4Infit(n1).rho())+ 2 * pmis1.pz() * (p4Infit(n1)).pz() * p4Infit(n1).m()/((p4Infit(n1)).rho() * p4Infit(n1).rho());
01949                                                   ta[0][3] = 2 * pmis1.e()  - 2 * pmis1.px() * (p4Infit(n1)).px() * p4Infit(n1).e() /((p4Infit(n1)).rho() * p4Infit(n1).rho())- 2 * pmis1.py() * (p4Infit(n1)).py() * p4Infit(n1).e()/((p4Infit(n1)).rho() * p4Infit(n1).rho())- 2 * pmis1.pz() * (p4Infit(n1)).pz() * p4Infit(n1).e()/((p4Infit(n1)).rho() * p4Infit(n1).rho());
01950                                                   setA(m_nc,pa,ta);
01951                                                   setAT(pa, m_nc, ta.T());
01952                                                   break;
01953                                               }
01954                                     case 2 :  {
01955                                                   HepMatrix tb(1, 4, 0);
01956                                                   tb[0][0] = -2 * pmis1.px();
01957                                                   tb[0][1] = -2 * pmis1.py();
01958                                                   tb[0][2] = -2 * pmis1.pz();
01959                                                   tb[0][3] = 2 * pmis1.e();
01960                                                   setB(m_nc,pb,tb);
01961                                                   setBT(pb,m_nc,tb.T());
01962                                                   break;
01963                                               }
01964                                     case 3 :  {
01965                                                   HepMatrix ta(1, 1, 0);
01966                                                   ta[0][0] = 2 * pmis1.e() * (p4Infit(n1).m()/p4Infit(n1).e());
01967                                                   setA(m_nc,pa,ta);
01968                                                   setAT(pa, m_nc, ta.T());
01969                                                   HepMatrix tb(1, 3, 0);
01970                                                   tb[0][0] = -2 * pmis1.px();
01971                                                   tb[0][1] = -2 * pmis1.py();
01972                                                   tb[0][2] = -2 * pmis1.pz();
01973                                                   setB(m_nc,pb,tb);
01974                                                   setBT(pb,m_nc,tb.T());
01975                                                   break;
01976                                               }
01977                                     case 4 :  {
01978                                                   HepMatrix ta(1, 2, 0);
01979                                                   double a1 = lambda * lambda + 1;
01980                                                   ta[0][0] = 2 * pmis1.px() * p4Infit(n1).py() - 2 * pmis1.py() * p4Infit(n1).px();
01981                                                   ta[0][1] = 2 * pmis1.px() * p4Infit(n1).px() * lambda/sqrt(a1) + 2 * pmis1.py() * (p4Infit(n1)).py() * lambda/sqrt(a1) - 2 * pmis1.pz() * (p4Infit(n1)).pz() /(sqrt(a1) * lambda);
01982                                                   setA(m_nc,pa,ta);
01983                                                   setAT(pa, m_nc, ta.T());
01984 
01985                                                   HepMatrix tb(1, 2, 0);
01986                                                   tb[0][0] =  2 * pmis1.px() * (p4Infit(n1)).px() * p4Infit(n1).m() /((p4Infit(n1)).rho() * p4Infit(n1).rho()) + 2 * pmis1.py() * (p4Infit(n1)).py() * p4Infit(n1).m()/((p4Infit(n1)).rho() * p4Infit(n1).rho())+ 2 * pmis1.pz() * (p4Infit(n1)).pz() * p4Infit(n1).m()/((p4Infit(n1)).rho() * p4Infit(n1).rho());
01987                                                   tb[0][1] = 2 * pmis1.e()  - 2 * pmis1.px() * (p4Infit(n1)).px() * p4Infit(n1).e() /((p4Infit(n1)).rho() * p4Infit(n1).rho())- 2 * pmis1.py() * (p4Infit(n1)).py() * p4Infit(n1).e()/((p4Infit(n1)).rho() * p4Infit(n1).rho())- 2 * pmis1.pz() * (p4Infit(n1)).pz() * p4Infit(n1).e()/((p4Infit(n1)).rho() * p4Infit(n1).rho());
01988                                                   setB(m_nc,pb,tb);
01989                                                   setBT(pb, m_nc, tb.T());
01990                                                   break;
01991                                               }
01992                                     case 5 :  {
01993                                                   HepMatrix ta(1, 3, 0);
01994                                                   double a1 = lambda * lambda + 1;
01995                                                   ta[0][0] = 2 * pmis1.px() * p4Infit(n1).py() - 2 * pmis1.py() * p4Infit(n1).px();
01996                                                   ta[0][1] = 2 * pmis1.px() * p4Infit(n1).px() * lambda/sqrt(a1)+ 2 * pmis1.py() * (p4Infit(n1)).py() * lambda/sqrt(a1) - 2 * pmis1.pz() * (p4Infit(n1)).pz() /(sqrt(a1) * lambda);
01997                                                   ta[0][2] = 2 * pmis1.e() * (p4Infit(n1)).m()/(p4Infit(n1)).e();
01998                                                   setA(m_nc,pa,ta);
01999                                                   setAT(pa, m_nc, ta.T());
02000                                                   HepMatrix tb(1, 1, 0);
02001                                                   tb[0][0] = 2 * pmis1.e() * (p4Infit(n1)).rho()/(p4Infit(n1)).e() - 2 * pmis1.px() * (p4Infit(n1)).px()/(p4Infit(n1)).rho() - 2 * pmis1.py() * (p4Infit(n1)).py()/(p4Infit(n1)).rho() - 2 * pmis1.pz() * (p4Infit(n1)).pz()/(p4Infit(n1)).rho();
02002                                                   setB(m_nc,pb,tb);
02003                                                   setBT(pb, m_nc, tb.T());
02004                                                   break;
02005                                               }       
02006                                     case 7: {
02007                                                 HepMatrix ta(1, NTRKPAR, 0);
02008                                                 ta[0][0] = -2 * pmis1.px() + 2 * pmis1.e() * p4Infit(n1).px() / p4Infit(n1).e();
02009                                                 ta[0][1] = -2 * pmis1.py() + 2 * pmis1.e() * p4Infit(n1).py() / p4Infit(n1).e();
02010                                                 ta[0][2] = -2 * pmis1.pz() + 2 * pmis1.e() * p4Infit(n1).pz() / p4Infit(n1).e();
02011                                                 ta[0][3] = 2 * pmis1.e() *  p4Infit(n1).m() / p4Infit(n1).e();
02012                                                 setA(m_nc,pa,ta);
02013                                                 setAT(pa, m_nc, ta.T());
02014                                                 break;
02015                                             }
02016                                 }   
02017                             }
02018 
02019                             for(int i = isiz; i < isiz+jsiz; i++) {
02020                                 int n2 = (kc.Ltrk())[i];
02021                                 double lambda = p4Infit(n2).pz()/p4Infit(n2).perp();
02022                                 int pa = mappositionA()[n2] + 1;
02023                                 int pb = mappositionB()[n2] + 1;
02024                                 int ptype = mapkinematic()[n2];
02025                                 switch(ptype) {
02026                                     case 0 : {
02027                                                  HepMatrix ta(1, NTRKPAR, 0);
02028                                                  ta[0][0] = -2 * pmis2.px() + 2 * pmis2.e() * p4Infit(n2).px() / p4Infit(n2).e();
02029                                                  ta[0][1] = -2 * pmis2.py() + 2 * pmis2.e() * p4Infit(n2).py() / p4Infit(n2).e();
02030                                                  ta[0][2] = -2 * pmis2.pz() + 2 * pmis2.e() * p4Infit(n2).pz() / p4Infit(n2).e();
02031                                                  ta[0][3] = 2 * pmis2.e() * p4Infit(n2).m() / p4Infit(n2).e();
02032                                                  setA(m_nc,pa,-ta);
02033                                                  setAT(pa, m_nc, -ta.T());
02034                                                  break;
02035                                              }
02036                                     case 1 :  {
02037                                                   HepMatrix ta(1, NTRKPAR, 0);
02038                                                   double a1 = lambda * lambda + 1;
02039                                                   ta[0][0] = 2 * pmis2.px() * p4Infit(n2).py() - 2 * pmis2.py() * p4Infit(n2).px();
02040                                                   ta[0][1] = 2 * pmis2.px() * p4Infit(n2).px() * lambda/sqrt(a1) + 2 * pmis2.py() * (p4Infit(n2)).py() * lambda/sqrt(a1) - 2 * pmis2.pz() * (p4Infit(n2)).pz() /(sqrt(a1) * lambda);
02041                                                   ta[0][2] =  2 * pmis2.px() * (p4Infit(n2)).px() * p4Infit(n2).m() /((p4Infit(n2)).rho() * p4Infit(n2).rho()) + 2 * pmis2.py() * (p4Infit(n2)).py() * p4Infit(n2).m()/((p4Infit(n2)).rho() * p4Infit(n2).rho())+ 2 * pmis2.pz() * (p4Infit(n2)).pz() * p4Infit(n2).m()/((p4Infit(n2)).rho() * p4Infit(n2).rho());
02042                                                   ta[0][3] = 2 * pmis2.e()  - 2 * pmis2.px() * (p4Infit(n2)).px() * p4Infit(n2).e() /((p4Infit(n2)).rho() * p4Infit(n2).rho())- 2 * pmis2.py() * (p4Infit(n2)).py() * p4Infit(n2).e()/((p4Infit(n2)).rho() * p4Infit(n2).rho())- 2 * pmis2.pz() * (p4Infit(n2)).pz() * p4Infit(n2).e()/((p4Infit(n2)).rho() * p4Infit(n2).rho());
02043                                                   setA(m_nc,pa,-ta);
02044                                                   setAT(pa, m_nc, -ta.T());
02045                                                   break;
02046                                               }
02047                                     case 2 :  {
02048                                                   HepMatrix tb(1, 4, 0);
02049                                                   tb[0][0] = -2 * pmis2.px();
02050                                                   tb[0][1] = -2 * pmis2.py();
02051                                                   tb[0][2] = -2 * pmis2.pz();
02052                                                   tb[0][3] = 2 * pmis2.e();
02053                                                   setB(m_nc,pb,-tb);
02054                                                   setBT(pb,m_nc,-tb.T());
02055                                                   break;
02056                                               }
02057                                     case 3 :  {
02058                                                   HepMatrix ta(1, 1, 0);
02059                                                   ta[0][0] = 2 * pmis2.e() * (p4Infit(n2).m()/p4Infit(n2).e());
02060                                                   setA(m_nc,pa,-ta);
02061                                                   setAT(pa, m_nc, -ta.T());
02062                                                   HepMatrix tb(1, 3, 0);
02063                                                   tb[0][0] = -2 * pmis2.px();
02064                                                   tb[0][1] = -2 * pmis2.py();
02065                                                   tb[0][2] = -2 * pmis2.pz();
02066                                                   setB(m_nc,pb,-tb);
02067                                                   setBT(pb,m_nc,-tb.T());
02068                                                   break;
02069                                               }
02070                                     case 4 :  {
02071                                                   HepMatrix ta(1, 2, 0);
02072                                                   double a1 = lambda * lambda + 1;
02073                                                   ta[0][0] = 2 * pmis2.px() * p4Infit(n2).py() - 2 * pmis2.py() * p4Infit(n2).px();
02074                                                   ta[0][1] = 2 * pmis2.px() * p4Infit(n2).px() * lambda/sqrt(a1) + 2 * pmis2.py() * (p4Infit(n2)).py() * lambda/sqrt(a1) - 2 * pmis2.pz() * (p4Infit(n2)).pz() /(sqrt(a1) * lambda);
02075                                                   setA(m_nc,pa,-ta);
02076                                                   setAT(pa, m_nc, -ta.T());
02077 
02078                                                   HepMatrix tb(1, 2, 0);
02079                                                   tb[0][0] =  2 * pmis2.px() * (p4Infit(n2)).px() * p4Infit(n2).m() /((p4Infit(n2)).rho() * p4Infit(n2).rho()) + 2 * pmis2.py() * (p4Infit(n2)).py() * p4Infit(n2).m()/((p4Infit(n2)).rho() * p4Infit(n2).rho())+ 2 * pmis2.pz() * (p4Infit(n2)).pz() * p4Infit(n2).m()/((p4Infit(n2)).rho() * p4Infit(n2).rho());
02080                                                   tb[0][1] = 2 * pmis2.e()  - 2 * pmis2.px() * (p4Infit(n2)).px() * p4Infit(n2).e() /((p4Infit(n2)).rho() * p4Infit(n2).rho())- 2 * pmis2.py() * (p4Infit(n2)).py() * p4Infit(n2).e()/((p4Infit(n2)).rho() * p4Infit(n2).rho())- 2 * pmis2.pz() * (p4Infit(n2)).pz() * p4Infit(n2).e()/((p4Infit(n2)).rho() * p4Infit(n2).rho());
02081                                                   setB(m_nc,pb,-tb);
02082                                                   setBT(pb, m_nc, -tb.T());
02083                                                   break;
02084                                               }
02085                                     case 5 :  {
02086                                                   HepMatrix ta(1, 3, 0);
02087                                                   double a1 = lambda * lambda + 1;
02088                                                   ta[0][0] = 2 * pmis2.px() * p4Infit(n2).py() - 2 * pmis2.py() * p4Infit(n2).px();
02089                                                   ta[0][1] = 2 * pmis2.px() * p4Infit(n2).px() * lambda/sqrt(a1) + 2 * pmis2.py() * (p4Infit(n2)).py() * lambda/sqrt(a1) - 2 * pmis2.pz() * (p4Infit(n2)).pz() /(sqrt(a1) * lambda);
02090                                                   ta[0][2] = 2 * pmis2.e() * (p4Infit(n2)).m()/(p4Infit(n2)).e();
02091                                                   setA(m_nc,pa,-ta);
02092                                                   setAT(pa, m_nc, -ta.T());
02093                                                   HepMatrix tb(1, 1, 0);
02094                                                   tb[0][0] = 2 * pmis2.e() * (p4Infit(n2)).rho()/(p4Infit(n2)).e() - 2 * pmis2.px() * (p4Infit(n2)).px()/(p4Infit(n2)).rho() - 2 * pmis2.py() * (p4Infit(n2)).py()/(p4Infit(n2)).rho() - 2 * pmis2.pz() * (p4Infit(n2)).pz()/(p4Infit(n2)).rho();
02095                                                   setB(m_nc,pb,-tb);
02096                                                   setBT(pb, m_nc, -tb.T());
02097                                                   break;
02098                                               }       
02099                                     case 7 : {
02100                                                  HepMatrix ta(1, NTRKPAR, 0);
02101                                                  ta[0][0] = -2 * pmis2.px() + 2 * pmis2.e() * p4Infit(n2).px() / p4Infit(n2).e();
02102                                                  ta[0][1] = -2 * pmis2.py() + 2 * pmis2.e() * p4Infit(n2).py() / p4Infit(n2).e();
02103                                                  ta[0][2] = -2 * pmis2.pz() + 2 * pmis2.e() * p4Infit(n2).pz() / p4Infit(n2).e();
02104                                                  ta[0][3] = 2 * pmis2.e() * p4Infit(n2).m() / p4Infit(n2).e();
02105                                                  setA(m_nc,pa,-ta);
02106                                                  setAT(pa, m_nc, -ta.T());
02107                                                  break;
02108                                              }
02109                                 }   
02110                             } 
02111 
02112 
02113 
02114 
02115                             HepVector dc(1, 0);
02116                             dc[0] = pmis1.m2() - pmis2.m2();
02117                             m_G[m_nc] = dc[0];
02118 
02119                             m_nc+=1;
02120 
02121                             break;
02122                         }
02123 
02124 
02125         case FourMomentum: 
02126         default: {
02127                      //
02128                      //  px - p4x = 0
02129                      //  py - p4y = 0
02130                      //  pz - p4z = 0
02131                      //  e  - p4e = 0
02132                      //
02133                      HepLorentzVector p4 = kc.p4();
02134                      HepLorentzVector pmis;
02135                      for(unsigned int j = 0; j< (kc.Ltrk()).size(); j++){
02136                          int n = (kc.Ltrk())[j];
02137                          pmis = pmis + p4Infit(n);
02138                      }
02139                      for(unsigned int j = 0; j < (kc.Ltrk()).size(); j++) {
02140                          int n = (kc.Ltrk())[j];
02141                          double lambda = p4Infit(n).pz()/p4Infit(n).perp();
02142                          int pa = mappositionA()[n] + 1;
02143                          int pb = mappositionB()[n] + 1;
02144                          int ptype = mapkinematic()[n];
02145                          switch(ptype) {
02146                              case 0 : {
02147                                           HepMatrix ta(kc.nc(), NTRKPAR, 0);
02148                                           ta[0][0] = 1.0;
02149                                           ta[1][1] = 1.0;
02150                                           ta[2][2] = 1.0;
02151                                           ta[3][0] = p4Infit(n).px() / p4Infit(n).e();
02152                                           ta[3][1] = p4Infit(n).py() / p4Infit(n).e();
02153                                           ta[3][2] = p4Infit(n).pz() / p4Infit(n).e();
02154                                           ta[3][3] = p4Infit(n).m() / p4Infit(n).e();
02155                                           setA(m_nc, pa, ta);
02156                                           setAT(pa, m_nc, ta.T());
02157                                           break;
02158                                       }
02159                              case 1 : {
02160                                           HepMatrix ta(kc.nc(), NTRKPAR, 0);
02161                                           ta[0][0] = -p4Infit(n).py();
02162                                           ta[0][1] = -p4Infit(n).px()*(lambda/(1 + lambda*lambda));
02163                                           ta[0][2] = -p4Infit(n).m()*p4Infit(n).px()/(p4Infit(n).rho()*p4Infit(n).rho());
02164                                           ta[0][3] = p4Infit(n).e()*p4Infit(n).px()/(p4Infit(n).rho()*p4Infit(n).rho());
02165                                           ta[1][0] = p4Infit(n).px();
02166                                           ta[1][1] = -p4Infit(n).py()*(lambda/(1 + lambda*lambda));
02167                                           ta[1][2] = -p4Infit(n).m()*p4Infit(n).py()/(p4Infit(n).rho()*p4Infit(n).rho());
02168                                           ta[1][3] = p4Infit(n).e()*p4Infit(n).py()/(p4Infit(n).rho()*p4Infit(n).rho());
02169                                           ta[2][0] = 0;
02170                                           ta[2][1] = p4Infit(n).pz()*(1/(lambda * (1 + lambda*lambda)));
02171                                           ta[2][2] = -p4Infit(n).m()*p4Infit(n).pz()/(p4Infit(n).rho()*p4Infit(n).rho());
02172                                           ta[2][3] = p4Infit(n).e()*p4Infit(n).pz()/(p4Infit(n).rho()*p4Infit(n).rho());
02173                                           ta[3][0] = 0;
02174                                           ta[3][1] = 0;
02175                                           ta[3][2] = 0;
02176                                           ta[3][3] = 1.0;
02177                                           setA(m_nc, pa, ta);
02178                                           setAT(pa, m_nc, ta.T()); 
02179                                           break;
02180                                       }
02181                              case 2 : {
02182                                           HepMatrix tb(kc.nc(), 4, 0);
02183                                           tb[0][0] = 1.0;
02184                                           tb[1][1] = 1.0;
02185                                           tb[2][2] = 1.0;
02186                                           tb[3][3] = 1.0;
02187                                           setB(m_nc, pb, tb);
02188                                           setBT(pb, m_nc, tb.T());
02189                                           break;
02190                                       }
02191                              case 3 : {
02192                                           HepMatrix ta(kc.nc(), 1, 0);
02193                                           ta[0][0] = -p4Infit(n).m()*p4Infit(n).px()/(p4Infit(n).rho()*p4Infit(n).rho());
02194                                           ta[1][0] = -p4Infit(n).m()*p4Infit(n).py()/(p4Infit(n).rho()*p4Infit(n).rho());
02195                                           ta[2][0] = -p4Infit(n).m()*p4Infit(n).pz()/(p4Infit(n).rho()*p4Infit(n).rho());
02196                                           ta[3][0] = 0;
02197                                           setA(m_nc, pa, ta);
02198                                           setAT(pa, m_nc, ta.T()); 
02199 
02200                                           HepMatrix tb(kc.nc(), 3, 0);
02201                                           tb[0][0] = -p4Infit(n).py();
02202                                           tb[0][1] = -p4Infit(n).px()*(lambda/(1 + lambda*lambda));
02203                                           tb[0][2] = p4Infit(n).e()*p4Infit(n).px()/(p4Infit(n).rho()*p4Infit(n).rho());
02204                                           tb[1][0] = p4Infit(n).px();
02205                                           tb[1][1] = -p4Infit(n).py()*(lambda/(1 + lambda*lambda));
02206                                           tb[1][2] = p4Infit(n).e()*p4Infit(n).py()/(p4Infit(n).rho()*p4Infit(n).rho());
02207                                           tb[2][0] = 0;
02208                                           tb[2][1] = p4Infit(n).pz()*(1/(lambda * (1 + lambda*lambda)));
02209                                           tb[2][2] = p4Infit(n).e()*p4Infit(n).pz()/(p4Infit(n).rho()*p4Infit(n).rho());
02210                                           tb[3][0] = 0;
02211                                           tb[3][1] = 0;
02212                                           tb[3][2] = 1.0;
02213 
02214                                           setB(m_nc, pb, tb);
02215                                           setBT(pb, m_nc, tb.T());
02216 
02217                                           break;
02218                                       }
02219                              case 4 : {
02220                                           HepMatrix ta(kc.nc(), 2, 0);
02221                                           ta[0][0] = -p4Infit(n).py(); 
02222                                           ta[0][1] = -p4Infit(n).px()*(lambda/sqrt(1 + lambda*lambda));
02223                                           ta[1][0] =  p4Infit(n).px(); 
02224                                           ta[1][1] = -p4Infit(n).py()*(lambda/sqrt(1 + lambda*lambda));
02225                                           ta[2][0] =  0;
02226                                           ta[2][1] =  p4Infit(n).pz()*(1/(lambda * (1 + lambda*lambda)));
02227 
02228                                           setA(m_nc, pa, ta);
02229                                           setAT(pa, m_nc, ta.T()); 
02230 
02231                                           HepMatrix tb(kc.nc(), 2, 0);
02232                                           //  tb[3][0] = p4Infit(n).m()/p4Infit(n).e();
02233                                           //  tb[0][1] = p4Infit(n).px()/p4Infit(n).rho();
02234                                           //  tb[1][1] = p4Infit(n).py()/p4Infit(n).rho();
02235                                           //  tb[2][1] = p4Infit(n).pz()/p4Infit(n).rho();
02236                                           //  tb[3][1] = p4Infit(n).rho()/p4Infit(n).e();
02237                                           tb[0][0] = -p4Infit(n).m()*p4Infit(n).px()/(p4Infit(n).rho()*p4Infit(n).rho());
02238                                           tb[1][0] = -p4Infit(n).m()*p4Infit(n).py()/(p4Infit(n).rho()*p4Infit(n).rho());
02239                                           tb[2][0] = -p4Infit(n).m()*p4Infit(n).pz()/(p4Infit(n).rho()*p4Infit(n).rho());
02240                                           tb[0][1] = p4Infit(n).e()*p4Infit(n).px()/(p4Infit(n).rho()*p4Infit(n).rho());
02241                                           tb[1][1] = p4Infit(n).e()*p4Infit(n).py()/(p4Infit(n).rho()*p4Infit(n).rho());
02242                                           tb[2][1] = p4Infit(n).e()*p4Infit(n).pz()/(p4Infit(n).rho()*p4Infit(n).rho());
02243                                           tb[3][1] = 1;
02244                                           setB(m_nc, pb, tb);
02245                                           setBT(pb, m_nc, tb.T());
02246                                           break;
02247                                       }
02248                              case 5 : {
02249                                           HepMatrix ta(kc.nc(), 3, 0);
02250                                           ta[0][0] = -p4Infit(n).py();
02251                                           ta[0][1] = -p4Infit(n).px()*(lambda/sqrt(1 + lambda*lambda));
02252                                           ta[1][0] =  p4Infit(n).px(); 
02253                                           ta[1][1] = -p4Infit(n).py()*(lambda/sqrt(1 + lambda*lambda));
02254                                           ta[2][0] =  0;
02255                                           ta[2][1] =  p4Infit(n).pz()*(1/(lambda * (1 + lambda*lambda)));
02256                                           ta[3][2] =  p4Infit(n).m()/p4Infit(n).e();
02257                                           setA(m_nc, pa, ta);
02258                                           setAT(pa, m_nc, ta.T()); 
02259 
02260                                           HepMatrix tb(kc.nc(), 1, 0);
02261                                           tb[0][0] = p4Infit(n).px()/p4Infit(n).rho();
02262                                           tb[1][0] = p4Infit(n).py()/p4Infit(n).rho();
02263                                           tb[2][0] = p4Infit(n).pz()/p4Infit(n).rho();
02264                                           tb[3][0] = p4Infit(n).rho()/p4Infit(n).e();
02265                                           setB(m_nc, pb, tb);
02266                                           setBT(pb, m_nc, tb.T());
02267                                           break;
02268                                       }
02269                              case 6 : {
02270                                           double ptrk = m_p.sub(pa+3, pa+3)[0];
02271                                           double x = m_p.sub(pa, pa)[0] - m_q.sub(pb, pb)[0];
02272                                           double y = m_p.sub(pa+1, pa+1)[0] - m_q.sub(pb+1, pb+1)[0];
02273                                           double z = m_p.sub(pa+2, pa+2)[0] - m_q.sub(pb+2, pb+2)[0];
02274                                           double R = sqrt(x*x + y*y + z*z);
02275                                           HepMatrix ta(kc.nc(), 4, 0);
02276                                           ta[0][0] = ptrk*(y*y + z*z)/(R*R*R);
02277                                           ta[0][1] = -ptrk*x*y/(R*R*R);
02278                                           ta[0][2] = -ptrk*x*z/(R*R*R);
02279                                           ta[0][3] = x/R;
02280                                           ta[1][0] = -ptrk*y*x/(R*R*R);
02281                                           ta[1][1] = ptrk*(x*x + z*z)/(R*R*R);
02282                                           ta[1][2] = -ptrk*y*z/(R*R*R);
02283                                           ta[1][3] = y/R;
02284                                           ta[2][0] = -ptrk*z*x/(R*R*R);
02285                                           ta[2][1] = -ptrk*z*y/(R*R*R);
02286                                           ta[2][2] = ptrk*(x*x + y*y)/(R*R*R);
02287                                           ta[2][3] = z/R;
02288                                           ta[3][3] = 1;
02289                                           setA(m_nc, pa, ta);
02290                                           setAT(pa, m_nc, ta.T());
02291 
02292                                           HepMatrix tb(kc.nc(), 3, 0);
02293                                           tb[0][0] = -ptrk*(y*y + z*z)/(R*R*R);
02294                                           tb[0][1] = ptrk*x*y/(R*R*R);
02295                                           tb[0][2] = ptrk*x*z/(R*R*R);
02296                                           tb[1][0] = ptrk*y*x/(R*R*R);
02297                                           tb[1][1] = -ptrk*(x*x + z*z)/(R*R*R);
02298                                           tb[1][2] = ptrk*y*z/(R*R*R);
02299                                           tb[2][0] = ptrk*z*x/(R*R*R);
02300                                           tb[2][1] = ptrk*z*y/(R*R*R);
02301                                           tb[2][2] = -ptrk*(x*x + y*y)/(R*R*R);
02302                                           HepMatrix tbp = m_B.sub(1, 4, 1, 3);
02303                                           tb = tbp + tb;
02304                                           setB(m_nc, pb, tb);
02305                                           setBT(pb, m_nc, tb.T());
02306                                           break;
02307                                       }
02308                              case 7 : {
02309                                           HepMatrix ta(kc.nc(), NTRKPAR, 0);
02310                                           ta[0][0] = 1.0;
02311                                           ta[1][1] = 1.0;
02312                                           ta[2][2] = 1.0;
02313                                           ta[3][0] = p4Infit(n).px() / p4Infit(n).e();
02314                                           ta[3][1] = p4Infit(n).py() / p4Infit(n).e();
02315                                           ta[3][2] = p4Infit(n).pz() / p4Infit(n).e();
02316                                           ta[3][3] = p4Infit(n).m() / p4Infit(n).e();
02317                                           setA(m_nc, pa, ta);
02318                                           setAT(pa, m_nc, ta.T());
02319                                           break;
02320                                       }
02321                          }
02322                      }  
02323 
02324                      HepVector dc(kc.nc(), 0);
02325                      dc[0] = pmis.px()  - p4.px();
02326                      dc[1] = pmis.py()  - p4.py();
02327                      dc[2] = pmis.pz()  - p4.pz();
02328                      dc[3] = pmis.e()  - p4.e();
02329                      for(int i = 0; i < kc.nc(); i++) {
02330                          m_G[m_nc+i] = dc[i];
02331                      }           
02332 
02333                      m_nc += 4;
02334 
02335                      break;
02336                  }
02337     }
02338 }

KalmanKinematicFit Class Reference

Public Member Functions

Static Public Member Functions

Private Member Functions

Private Attributes

Static Private Attributes

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation


Public Member Functions
void	AddEqualMass (int number, std::vector< int > tlis1, std::vector< int > tlis2)
void	AddEqualMass (int number, std::vector< int > tlis1, std::vector< int > tlis2)
void	AddFourMomentum (int number, double etot)
void	AddFourMomentum (int number, HepLorentzVector p4)
void	AddFourMomentum (int number, double etot)
void	AddFourMomentum (int number, HepLorentzVector p4)
std::vector< int >	AddList (int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8, int n9, int n10, int n11, int n12)
std::vector< int >	AddList (int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8, int n9, int n10, int n11)
std::vector< int >	AddList (int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8, int n9, int n10)
std::vector< int >	AddList (int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8, int n9)
std::vector< int >	AddList (int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8)
std::vector< int >	AddList (int n1, int n2, int n3, int n4, int n5, int n6, int n7)
std::vector< int >	AddList (int n1, int n2, int n3, int n4, int n5, int n6)
std::vector< int >	AddList (int n1, int n2, int n3, int n4, int n5)
std::vector< int >	AddList (int n1, int n2, int n3, int n4)
std::vector< int >	AddList (int n1, int n2, int n3)
std::vector< int >	AddList (int n1, int n2)
std::vector< int >	AddList (int n1)
std::vector< int >	AddList (int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8, int n9, int n10, int n11, int n12)
std::vector< int >	AddList (int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8, int n9, int n10, int n11)
std::vector< int >	AddList (int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8, int n9, int n10)
std::vector< int >	AddList (int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8, int n9)
std::vector< int >	AddList (int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8)
std::vector< int >	AddList (int n1, int n2, int n3, int n4, int n5, int n6, int n7)
std::vector< int >	AddList (int n1, int n2, int n3, int n4, int n5, int n6)
std::vector< int >	AddList (int n1, int n2, int n3, int n4, int n5)
std::vector< int >	AddList (int n1, int n2, int n3, int n4)
std::vector< int >	AddList (int n1, int n2, int n3)
std::vector< int >	AddList (int n1, int n2)
std::vector< int >	AddList (int n1)
void	AddMissTrack (const int number, const HepLorentzVector p4)
void	AddMissTrack (const int number, const RecEmcShower *trk)
void	AddMissTrack (const int number, const double mass, const RecEmcShower *trk)
void	AddMissTrack (const int number, const double mass, const HepLorentzVector p4)
void	AddMissTrack (const int number, const double mass)
void	AddMissTrack (const int number, const HepLorentzVector p4)
void	AddMissTrack (const int number, const RecEmcShower *trk)
void	AddMissTrack (const int number, const double mass, const RecEmcShower *trk)
void	AddMissTrack (const int number, const double mass, const HepLorentzVector p4)
void	AddMissTrack (const int number, const double mass)
void	AddResonance (int number, double mres, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8, int n9, int n10, int n11, int n12)
void	AddResonance (int number, double mres, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8, int n9, int n10, int n11)
void	AddResonance (int number, double mres, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8, int n9, int n10)
void	AddResonance (int number, double mres, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8, int n9)
void	AddResonance (int number, double mres, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8)
void	AddResonance (int number, double mres, int n1, int n2, int n3, int n4, int n5, int n6, int n7)
void	AddResonance (int number, double mres, int n1, int n2, int n3, int n4, int n5, int n6)
void	AddResonance (int number, double mres, int n1, int n2, int n3, int n4, int n5)
void	AddResonance (int number, double mres, int n1, int n2, int n3, int n4)
void	AddResonance (int number, double mres, int n1, int n2, int n3)
void	AddResonance (int number, double mres, int n1, int n2)
void	AddResonance (int number, double mres, int n1)
void	AddResonance (int number, double mres, std::vector< int > tlis)
void	AddResonance (int number, double mres, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8, int n9, int n10, int n11, int n12)
void	AddResonance (int number, double mres, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8, int n9, int n10, int n11)
void	AddResonance (int number, double mres, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8, int n9, int n10)
void	AddResonance (int number, double mres, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8, int n9)
void	AddResonance (int number, double mres, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8)
void	AddResonance (int number, double mres, int n1, int n2, int n3, int n4, int n5, int n6, int n7)
void	AddResonance (int number, double mres, int n1, int n2, int n3, int n4, int n5, int n6)
void	AddResonance (int number, double mres, int n1, int n2, int n3, int n4, int n5)
void	AddResonance (int number, double mres, int n1, int n2, int n3, int n4)
void	AddResonance (int number, double mres, int n1, int n2, int n3)
void	AddResonance (int number, double mres, int n1, int n2)
void	AddResonance (int number, double mres, int n1)
void	AddResonance (int number, double mres, std::vector< int > tlis)
void	AddThreeMomentum (int number, Hep3Vector p3)
void	AddThreeMomentum (int number, Hep3Vector p3)
void	AddTotalEnergy (int number, double etot, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8, int n9, int n10, int n11, int n12)
void	AddTotalEnergy (int number, double etot, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8, int n9, int n10, int n11)
void	AddTotalEnergy (int number, double etot, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8, int n9, int n10)
void	AddTotalEnergy (int number, double etot, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8, int n9)
void	AddTotalEnergy (int number, double etot, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8)
void	AddTotalEnergy (int number, double etot, int n1, int n2, int n3, int n4, int n5, int n6, int n7)
void	AddTotalEnergy (int number, double etot, int n1, int n2, int n3, int n4, int n5, int n6)
void	AddTotalEnergy (int number, double etot, int n1, int n2, int n3, int n4, int n5)
void	AddTotalEnergy (int number, double etot, int n1, int n2, int n3, int n4)
void	AddTotalEnergy (int number, double etot, int n1, int n2, int n3)
void	AddTotalEnergy (int number, double etot, int n1, int n2)
void	AddTotalEnergy (int number, double etot, int n1)
void	AddTotalEnergy (int number, double etot, std::vector< int > lis)
void	AddTotalEnergy (int number, double etot, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8, int n9, int n10, int n11, int n12)
void	AddTotalEnergy (int number, double etot, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8, int n9, int n10, int n11)
void	AddTotalEnergy (int number, double etot, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8, int n9, int n10)
void	AddTotalEnergy (int number, double etot, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8, int n9)
void	AddTotalEnergy (int number, double etot, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8)
void	AddTotalEnergy (int number, double etot, int n1, int n2, int n3, int n4, int n5, int n6, int n7)
void	AddTotalEnergy (int number, double etot, int n1, int n2, int n3, int n4, int n5, int n6)
void	AddTotalEnergy (int number, double etot, int n1, int n2, int n3, int n4, int n5)
void	AddTotalEnergy (int number, double etot, int n1, int n2, int n3, int n4)
void	AddTotalEnergy (int number, double etot, int n1, int n2, int n3)
void	AddTotalEnergy (int number, double etot, int n1, int n2)
void	AddTotalEnergy (int number, double etot, int n1)
void	AddTotalEnergy (int number, double etot, std::vector< int > lis)
void	AddTotalMomentum (int number, double ptot, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8, int n9, int n10, int n11, int n12)
void	AddTotalMomentum (int number, double ptot, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8, int n9, int n10, int n11)
void	AddTotalMomentum (int number, double ptot, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8, int n9, int n10)
void	AddTotalMomentum (int number, double ptot, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8, int n9)
void	AddTotalMomentum (int number, double ptot, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8)
void	AddTotalMomentum (int number, double ptot, int n1, int n2, int n3, int n4, int n5, int n6, int n7)
void	AddTotalMomentum (int number, double ptot, int n1, int n2, int n3, int n4, int n5, int n6)
void	AddTotalMomentum (int number, double ptot, int n1, int n2, int n3, int n4, int n5)
void	AddTotalMomentum (int number, double ptot, int n1, int n2, int n3, int n4)
void	AddTotalMomentum (int number, double ptot, int n1, int n2, int n3)
void	AddTotalMomentum (int number, double ptot, int n1, int n2)
void	AddTotalMomentum (int number, double ptot, int n1)
void	AddTotalMomentum (int number, double ptot, std::vector< int > lis)
void	AddTotalMomentum (int number, double ptot, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8, int n9, int n10, int n11, int n12)
void	AddTotalMomentum (int number, double ptot, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8, int n9, int n10, int n11)
void	AddTotalMomentum (int number, double ptot, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8, int n9, int n10)
void	AddTotalMomentum (int number, double ptot, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8, int n9)
void	AddTotalMomentum (int number, double ptot, int n1, int n2, int n3, int n4, int n5, int n6, int n7, int n8)
void	AddTotalMomentum (int number, double ptot, int n1, int n2, int n3, int n4, int n5, int n6, int n7)
void	AddTotalMomentum (int number, double ptot, int n1, int n2, int n3, int n4, int n5, int n6)
void	AddTotalMomentum (int number, double ptot, int n1, int n2, int n3, int n4, int n5)
void	AddTotalMomentum (int number, double ptot, int n1, int n2, int n3, int n4)
void	AddTotalMomentum (int number, double ptot, int n1, int n2, int n3)
void	AddTotalMomentum (int number, double ptot, int n1, int n2)
void	AddTotalMomentum (int number, double ptot, int n1)
void	AddTotalMomentum (int number, double ptot, std::vector< int > lis)
void	AddTrack (const int number, const WTrackParameter wtrk)
void	AddTrack (const int number, const double mass, const RecEmcShower *trk)
void	AddTrack (const int number, const double mass, const RecMdcKalTrack *trk)
void	AddTrack (const int number, const double mass, const RecMdcTrack *trk)
void	AddTrack (const int number, const WTrackParameter wtrk)
void	AddTrack (const int number, const double mass, const RecEmcShower *trk)
void	AddTrack (const int number, const double mass, const RecMdcKalTrack *trk)
void	AddTrack (const int number, const double mass, const RecMdcTrack *trk)
void	AddTrackVertex (const int number, const double mass, const RecEmcShower *trk)
void	AddTrackVertex (const int number, const double mass, const RecEmcShower *trk)
void	BuildVirtualParticle (int number)
void	BuildVirtualParticle (int number)
double	chisq (int n) const
double	chisq () const
double	chisq (int n) const
double	chisq () const
void	clearGammaShape ()
void	clearGammaShape ()
void	clearGammaShapeList ()
void	clearGammaShapeList ()
void	clearMapkinematic ()
void	clearMapkinematic ()
void	clearMappositionA ()
void	clearMappositionA ()
void	clearMappositionB ()
void	clearMappositionB ()
void	clearone ()
void	clearone ()
void	cleartwo ()
void	cleartwo ()
void	clearWTrackInfit ()
void	clearWTrackInfit ()
void	clearWTrackList ()
void	clearWTrackList ()
void	clearWTrackOrigin ()
void	clearWTrackOrigin ()
double	collideangle () const
double	collideangle () const
HepVector	cpu () const
HepVector	cpu () const
bool	dynamicerror () const
bool	dynamicerror () const
double	espread () const
double	espread () const
bool	Fit (int n)
bool	Fit ()
bool	Fit (int n)
bool	Fit ()
int	GammaShapeList (int n) const
std::vector< int >	GammaShapeList () const
int	GammaShapeList (int n) const
std::vector< int >	GammaShapeList () const
GammaShape	GammaShapeValue (int n) const
std::vector< GammaShape >	GammaShapeValue () const
GammaShape	GammaShapeValue (int n) const
std::vector< GammaShape >	GammaShapeValue () const
HepPoint3D	getBeamPosition () const
HepPoint3D	getBeamPosition () const
HepSymMatrix	getCInfit (int i) const
HepSymMatrix	getCInfit (int i) const
HepSymMatrix	getCOrigin (int i) const
HepSymMatrix	getCOrigin (int i) const
HepSymMatrix	getVBeamPosition () const
HepSymMatrix	getVBeamPosition () const
WTrackParameter	infit (int n) const
WTrackParameter	infit (int n) const
void	init ()
void	init ()
vector< int >	mapkinematic () const
vector< int >	mapkinematic () const
vector< int >	mappositionA () const
vector< int >	mappositionA () const
vector< int >	mappositionB () const
vector< int >	mappositionB () const
int	numberGammaShape () const
int	numberGammaShape () const
int	numberone () const
int	numberone () const
int	numbertwo () const
int	numbertwo () const
int	numberWTrack () const
int	numberWTrack () const
WTrackParameter	origin (int n) const
WTrackParameter	origin (int n) const
HepLorentzVector	pfit (int n) const
HepLorentzVector	pfit (int n) const
HepLorentzVector	pfit1 (int n)
HepLorentzVector	pfit1 (int n)
HepVector	pull (int n)
HepVector	pull (int n)
void	setBeamPosition (const HepPoint3D BeamPosition)
void	setBeamPosition (const HepPoint3D BeamPosition)
void	setChisqCut (const double chicut=200, const double chiter=0.05)
void	setChisqCut (const double chicut=200, const double chiter=0.05)
void	setCollideangle (const double collideangle=11e-3)
void	setCollideangle (const double collideangle=11e-3)
void	setDynamicerror (const bool dynamicerror=1)
void	setDynamicerror (const bool dynamicerror=1)
void	setEspread (const double espread=0.0009)
void	setEspread (const double espread=0.0009)
void	setFlag (const bool flag=1)
void	setFlag (const bool flag=1)
void	setGammaShape (const GammaShape gammashape)
void	setGammaShape (const int n, const GammaShape gammashape)
void	setGammaShape (const GammaShape gammashape)
void	setGammaShape (const int n, const GammaShape gammashape)
void	setGammaShapeList (const int n)
void	setGammaShapeList (const int n)
void	setIterNumber (const int niter=5)
void	setIterNumber (const int niter=5)
void	setMapkinematic (const int n)
void	setMapkinematic (const int n)
void	setMappositionA (const int n)
void	setMappositionA (const int n)
void	setMappositionB (const int n)
void	setMappositionB (const int n)
void	setTgraph (TGraph2D *graph2d)
void	setTgraph (TGraph2D *graph2d)
void	setVBeamPosition (const HepSymMatrix VBeamPosition)
void	setVBeamPosition (const HepSymMatrix VBeamPosition)
void	setWTrackInfit (const WTrackParameter wtrk)
void	setWTrackInfit (const int n, const WTrackParameter wtrk)
void	setWTrackInfit (const WTrackParameter wtrk)
void	setWTrackInfit (const int n, const WTrackParameter wtrk)
void	setWTrackList (const int n)
void	setWTrackList (const int n)
void	setWTrackOrigin (const WTrackParameter wtrk)
void	setWTrackOrigin (const int n, const WTrackParameter wtrk)
void	setWTrackOrigin (const WTrackParameter wtrk)
void	setWTrackOrigin (const int n, const WTrackParameter wtrk)
WTrackParameter	wTrackInfit (int n) const
std::vector< WTrackParameter >	wTrackInfit () const
WTrackParameter	wTrackInfit (int n) const
std::vector< WTrackParameter >	wTrackInfit () const
int	wTrackList (int n) const
std::vector< int >	wTrackList () const
int	wTrackList (int n) const
std::vector< int >	wTrackList () const
WTrackParameter	wTrackOrigin (int n) const
std::vector< WTrackParameter >	wTrackOrigin () const
WTrackParameter	wTrackOrigin (int n) const
std::vector< WTrackParameter >	wTrackOrigin () const
WTrackParameter	wVirtualTrack (int n) const
WTrackParameter	wVirtualTrack (int n) const
HepVector	xfit ()
HepVector	xfit ()
	~KalmanKinematicFit ()
	~KalmanKinematicFit ()
Static Public Member Functions
KalmanKinematicFit *	instance ()
KalmanKinematicFit *	instance ()
Private Member Functions
void	clearDDT ()
void	clearDDT ()
void	covMatrix (int n)
void	covMatrix (int n)
void	fit (int n)
void	fit ()
void	fit (int n)
void	fit ()
void	gda ()
void	gda ()
	KalmanKinematicFit ()
	KalmanKinematicFit ()
HepLorentzVector	p4Infit (int i) const
HepLorentzVector	p4Infit (int i) const
HepLorentzVector	p4Origin (int i) const
HepLorentzVector	p4Origin (int i) const
HepVector	pInfit (int i) const
HepVector	pInfit (int i) const
HepVector	pOrigin (int i) const
HepVector	pOrigin (int i) const
void	setA (int ic, int itk, const HepMatrix &p)
void	setA (int ic, int itk, const HepMatrix &p)
void	setAT (int itk, int ic, const HepMatrix &p)
void	setAT (int itk, int ic, const HepMatrix &p)
void	setB (int ic, int itk, const HepMatrix &p)
void	setB (int ic, int itk, const HepMatrix &p)
void	setBT (int itk, int ic, const HepMatrix &p)
void	setBT (int itk, int ic, const HepMatrix &p)
void	setCInfit (int i, const HepSymMatrix &D)
void	setCInfit (int i, const HepSymMatrix &D)
void	setCOrigin (int i, const HepSymMatrix &D)
void	setCOrigin (int i, const HepSymMatrix &D)
void	setDInfit (int i, const HepSymMatrix &D)
void	setDInfit (int i, const HepSymMatrix &D)
void	setDOrigin (int i, const HepSymMatrix &D)
void	setDOrigin (int i, const HepSymMatrix &D)
void	setDOriginInv (int i, const HepSymMatrix &Dinv)
void	setDOriginInv (int i, const HepSymMatrix &Dinv)
void	setPInfit (int i, const HepVector &p)
void	setPInfit (int i, const HepVector &p)
void	setPOrigin (int i, const HepVector &p)
void	setPOrigin (int i, const HepVector &p)
void	setQInfit (int i, const HepVector &q)
void	setQInfit (int i, const HepVector &q)
void	setQOrigin (int i, const HepVector &q)
void	setQOrigin (int i, const HepVector &q)
void	upCovmtx ()
void	upCovmtx ()
void	updateConstraints (KinematicConstraints kc)
void	updateConstraints (KinematicConstraints kc)
void	upTrkpar ()
void	upTrkpar ()
Private Attributes
HepMatrix	m_A
HepMatrix	m_AT
HepMatrix	m_B
HepMatrix	m_BT
HepSymMatrix	m_C
HepSymMatrix	m_C0
double	m_chi
double	m_chicut
std::vector< double >	m_chisq
std::vector< double >	m_chisq
double	m_chiter
double	m_collideangle
HepVector	m_cpu
HepSymMatrix	m_D
HepSymMatrix	m_D0
HepSymMatrix	m_D0inv
HepSymMatrix	m_Dinv
bool	m_dynamicerror
double	m_espread
bool	m_flag
HepVector	m_G
TGraph2D *	m_graph2d
TGraph2D *	m_graph2d
std::vector< KinematicConstraints >	m_kc
std::vector< KinematicConstraints >	m_kc
HepMatrix	m_KP
HepMatrix	m_KQ
int	m_nc
int	m_niter
int	m_nktrk
HepVector	m_p
HepVector	m_p0
HepVector	m_q
HepVector	m_q0
std::vector< WTrackParameter >	m_virtual_wtrk
std::vector< WTrackParameter >	m_virtual_wtrk
HepSymMatrix	m_Vm
HepSymMatrix	m_W
Static Private Attributes
const int	EqualMass = 64
const int	FourMomentum = 32
KalmanKinematicFit *	m_pointer
KalmanKinematicFit *	m_pointer = 0
const int	NKFPAR
const int	NTRKPAR = 4
const int	Position = 8
const int	Resonance = 1
const int	ThreeMomentum = 16
const int	TotalEnergy = 2
const int	TotalMomentum = 4