C++ TMatrixD::ResizeTo方法代码示例

   void CheckPoints(){
      std::vector<int> RemovedPoints;
      TMatrixD y_cov_new;
      int j = 0;
    
      for(int i = 0; i < y_val.size(); i++) {
         //         std::cout << "i: " << i << "   j: " << j << std::endl;
         if( y_val.at(i) == 0) {
            x_val.erase(x_val.begin()+i);
            y_val.erase(y_val.begin()+i);
            RemovedPoints.push_back(j);
            i = i-1;
         }
         j++;
      }    
    //   for(int i = 0; i < x_val.size(); i++) {
//           std::cout << "x: " << x_val.at(i) << std::endl;

//       }

      //  std::cout<< "Removed Points: " << RemovedPoints.size() << std::endl;
      //  std::cout<< "Remaining Points: " << x_val.size() << std::endl;

      y_cov_new.ResizeTo(x_val.size(),x_val.size());
      for( int i=0; i < x_val.size(); i++) {
         for(int k= 0; k < x_val.size(); k++) {
            y_cov_new(i,k) = y_cov(i+RemovedPoints.size(),k+RemovedPoints.size());
         }
      }
      y_cov.ResizeTo(0,0);
      y_cov.ResizeTo(x_val.size(),x_val.size());
      y_cov = y_cov_new;  
   }

double leptonic_fitter_algebraic::calc_MLL( const double* par, bool track_prob )
{
  double sB = par[0];
  if( _dbg > 219 ) cout<<"DBG220 lfa::calc for: "<<sB<<endl;
  update_non_nu_gen( sB );
  if( _dbg > 221 ) cout<<"DBG222 lfa::calc _genB: "<<lv2str(_genB)<<endl;

  double x1 = _x1_0 + 0.5*(_TmW_m2 - _b_m2*sB*sB) / (_denom*sB);
  double Z2 = _Z2_0 - 2*_x0*x1;
  double Z = (Z2 >= 0) ? TMath::Sqrt( Z2 ) : 0;
  if( _dbg > 224 ) cout<<"DBG225 lfa::calc Z2: "<<Z2<<endl;

  double entries[9] = { -Z*_y1/_x0,  0,  x1-_lep_p,     Z, 0, _y1,       0, Z, 0 };
  _Ebl.Use( 3, 3, entries );
  _Emat = _R_T;
  _Emat *= _Ebl;
    
  if( _dbg > 220 && Z2 >= 0 && _lep_p < 2 * _top_mass ) { // very energetic leptons can cause inaccuracies
    for( double iphi = 0; iphi < 4; ++iphi ) {
      double test_phi = iphi*TMath::Pi()/2;
      double tnums[3] = { TMath::Cos( test_phi ), TMath::Sin( test_phi ), 1. };
      TMatrixD test_tvec( 3, 1, tnums );
      TMatrixD test_nu( _Emat, TMatrixD::kMult, test_tvec );
      update_nu_and_decay_chain( test_nu );
      cout<<"DBG221 lfa::calc test_phi: "<<test_phi<<", test_nu: "<<lv2str(_genN)<<", -> W: "<<lv2str(_W)<<", t: "<<lv2str(_T)<<endl;
      assert( TMath::Abs( _W.M() - _W_mass ) < 1 && TMath::Abs( _T.M() - _top_mass ) < 1 );
    }
  }
  _dnu_mat = _Nu;
  _dnu_mat -= _Emat;

  bool tvec_set = false;
  if( Z > 0.002 ) { // ellipse is large enough to find the best point on it.
    _M_form.Transpose( _dnu_mat );
    _M_form *= _invFlatVar;
    _M_form *= _dnu_mat;
    if( _dbg > 221 ) cout<<"DBG222 lfa::calc _M_form: "<<m2str( _M_form )<<endl;

    // find the closest tvec
    std::vector< TMatrixD > extrimal_tvecs;
    _N = _M_form;
    _N *= _Q;
    _NT.Transpose( _N );
    _P = _N;
    _P += _NT;

    TMatrixD UP( _unit_circle, TMatrixD::kMult, _P );
    real_eigenvalues( UP );
    if( _dbg > 221 ) {
      cout<<"DBG222 lfa::calc _P: "<<m2str( _P )<<" -> the following "<<_real_eigs.size()<<" evals: ";
      for( unsigned int ire=0; ire < _real_eigs.size(); ++ire ) cout<<_real_eigs[ire]<<", ";
      cout<<endl;
    }
    if( _real_eigs.size() == 0 ) {
      cerr<<"ERROR: found no real eigenvectors of UP, which makes no sense."<<endl
	  <<"       Z: "<<Z<<" UP: "<<m2str( UP )<<endl;
      throw std::runtime_error( "leptonic_fitter_algebraic has numerical problems - found no real eigenvectors of UP.");
    }

    // find the degenrate matrix, if possible, choose one with intersecting lines
    double c22 = 1;
    for( unsigned int ire=0; ire < _real_eigs.size(); ++ire ) {
      _D = _unit_circle;
      _D *= - _real_eigs[ ire ];
      _D += _P;
    
      c22 = cofactor( _D, 2, 2 );
      if( _dbg > 221 ) cout<<"DBG222 lfa::calc e#: "<<ire<<" -> _D: "<<m2str( _D )<<" -> c22: "<<c22<<endl;
      if( c22 < 0 ) break;
    }

    // D is degenrate --> represents an equation which factorizes into two line equations.
    // Several cases which determine how we intersect these lines with the unit circle

    bool as_parallel_lines = c22 >= 0; // which means it's really 0, and positive value is due to numerical inaccuracies
                                       // so this is the first logical decision in the code. The actual first if is technical.
    if( _swapped ) { // will avoid D11=~0 in all but ~10-10 of the events
      static TMatrixD tmp;
      tmp.ResizeTo( _D );
      tmp =  _D;
      swap_x_y( tmp, _D );
    }
    double D11 = _D[1][1];
    if( D11 == 0 ) { // vertical lines, where y does not appear in the equations of the two lines
      cerr<<"ERROR: leptonic_fitter_algebraic does not handle the super-rare case of D00=D11==0..\n"
	  <<"ERROR: Estimated frequency is only 10^-10....\n"
	  <<"       Please use this as the first test case..."<<endl;

    } else { // normal case: D11 isn't 0, which means that y appears in the line equations

      if( as_parallel_lines ) {
	cout<<"Information: leptonic_fitter_algebraic encountered the \"0 probability\" case of parallel lines"<<endl;
	double c00 = cofactor( _D, 0, 0 );
	if( c00 > 0 ) {
	  cerr<<"ERROR: algorithm assumes for the \"0 probability\" case of positive c22, at least c00 is non-positive. D11!=0, Z: "
	      <<Z<<" c22: "<<c22<<" c00: "<<c00<<endl;
	  throw std::runtime_error( "leptonic_fitter_algebraic has numerical problems with c22 >= 0 and c00 > 0 (D11!=0)"); 
	}
	double sqrtNc00 = TMath::Sqrt( -c00 );
	double A = D11;
//.........这里部分代码省略.........

 void reset(){
    x_val.clear();
    y_val.clear();
    y_cov.ResizeTo(0,0);
    y_cov_inv.ResizeTo(0,0);
 }

//.........这里部分代码省略.........

   Double_t threshold = 0;
   Double_t Tthreshold = 1.;

   Double_t charge, Ccharge;
   Double_t Rm,Cm,Zcor,Zm;

   //Matrix for the charge map
   Double_t **padCharge=new Double_t*[numberOfRows];
   Double_t **padTime=new Double_t*[numberOfRows];
   Double_t **padHeight=new Double_t*[numberOfRows];
   for(Int_t j=0;j<numberOfRows;j++){
     padCharge[j]=new Double_t[numberOfColumns];
     padTime[j]=new Double_t[numberOfColumns];
   }
   
   //Double_t *padChargeTest=new Double_t[numberOfRows];
   //Double_t *padTimeTest=new Double_t[numberOfRows];

   //TMatrixD *padCharge=0;
   //TMatrixD *padTime=0;
   //TMatrixD *padHeight=0;

   //TMatrixD padCharge(32,64,2000);
   //TMatrixD padTime(32,64,2000);
   //TMatrixD padHeight(32,64,2000);

   //TMatrixD *padChargeTest=0;
   //TMatrixD *padTimeTest=0;
   TMatrixD padChargeTest;
   TMatrixD padTimeTest;

   TMatrixD a(32,64);
   padChargeTest.ResizeTo(a);
   padTimeTest.ResizeTo(a);

   TMatrixD *pointerCharge = &padChargeTest;
   TMatrixD *pointerTime = &padTimeTest;

   //Silicon charge
   //Double_t SilCharge[16];
   //Int_t SilID[16];

   TVectorD SilCharge;
   SilCharge.ResizeTo(16);
   TVectorD SilID;
   SilID.ResizeTo(16);

   TVectorD *pointerSilCharge=&SilCharge;
   TVectorD *pointerSilID=&SilID;
   //Double_t *SilChargePointer=&SilCharge;

   //Branching the out tree
   //out_tree->Branch("PadCharge","TMatrixD",&padCharge,64000,0);
   //out_tree->Branch("PadTime","TMatrixD",&padTime,64000,0);
   //out_tree->Branch("SilCharge",&SilCharge,64000,0);
   out_tree->Branch("PadCharge","TMatrixD",&pointerCharge,64000,0);
   out_tree->Branch("PadTime","TMatrixD",&pointerTime,64000,0);
   out_tree->Branch("SilCharge","TVectorD",&pointerSilCharge,64000,0);
   out_tree->Branch("SilID","TVectorD",&pointerSilID,64000,0);
   //out_tree->Branch("SilCharge",&SilCharge,"energy/D");


   // out_tree->Branch("Energy",&Qtot,"energy/D");
   // out_tree->Branch("HAngle",&HAngle,"angle/D");
   // out_tree->Branch("VAngle",&VAngle,"angle/D");