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

TLorentzVector leptonic_fitter_algebraic::lv_with_mass( const TMatrixD& mat, double mass )
{
  assert( mat.GetNcols() == 1 || mat.GetNrows() == 1 );
  TLorentzVector lv;
  const double *array = mat.GetMatrixArray();
  lv.SetXYZM( array[0], array[1], array[2], mass );
  return lv;
}

int addInQuadrature(const TMatrixD &sf, const TMatrixD &sfErr, const TMatrixD &sfSystRelErr, TMatrixD &totErr) {
  for (int ir=0; ir<sf.GetNrows(); ++ir) {
    for (int ic=0; ic<sf.GetNcols(); ++ic) {
      double err2=pow(sfErr(ir,ic),2) + pow(sf(ir,ic)*sfSystRelErr(ir,ic),2);
      totErr(ir,ic)=sqrt(err2);
    }
  }
  return 1;
}

void calculateInvertedMatrixErrors(TMatrixD &T, TMatrixD &TErrPos, TMatrixD &TErrNeg,
				   TMatrixD &TinvErr){

  // Calculate errors on the inverted matrix by the Monte Carlo
  // method

  Double_t det;
  int nRow = T.GetNrows();
  int nCol = T.GetNcols();
  TMatrixD TinvSum(nRow,nCol);
  TMatrixD TinvSumSquares(nRow,nCol);

  // Reset Matrix where we will be accumulating RMS/sigma:
  TinvSum        = 0;
  TinvSumSquares = 0;
  TinvErr        = 0;

  // Do many tries, accumulate RMS
  int N = 100000;
  for(int iTry = 0; iTry<N; iTry++){
    // Find the smeared matrix
    TMatrixD Tsmeared = T;
    for(int i = 0; i<nRow; i++){
      for(int j = 0; j<nCol; j++){
	double central = T(i,j);
	double sigPos = TErrPos(i,j);
	double sigNeg = TErrNeg(i,j);
	// Switch to symmetric errors: approximation, but much simpler
	double sig = (sigPos+sigNeg)/2.0;
	Tsmeared(i,j) = gRandom->Gaus(central,sig);
      }
    }
    // Find the inverted to smeared matrix
    TMatrixD TinvSmeared = Tsmeared;
    TinvSmeared.Invert(&det);
    // Accumulate sum and sum of squares for each element
    for(int i2 = 0; i2<nRow; i2++){
      for(int j2 = 0; j2<nCol; j2++){
	TinvSum       (i2,j2) += TinvSmeared(i2,j2);
	TinvSumSquares(i2,j2) += TinvSmeared(i2,j2)*TinvSmeared(i2,j2);
      }
    }
  }

  // Calculate the error matrix
  TMatrixD TinvAverage = TinvSum;
  for(int i = 0; i<nRow; i++){
    for(int j = 0; j<nCol; j++){
      TinvErr(i,j) = sqrt( TinvSumSquares(i,j)/double(N) 
			   - (TinvSum(i,j)/double(N))*(TinvSum(i,j)/double(N)) );
      TinvAverage(i,j) = TinvSum(i,j)/double(N);
    }
  }

  return;
}

//==========================================================================
//
//  Dump element content of TMatrixD and TVectorD
//
//==========================================================================
void dumpElements(TMatrixD& a)
{
  cout << endl << endl;
  const int nrows = a.GetNrows();
  const int ncols = a.GetNcols();
  if(nrows==ncols)
    cout << "determinent = " << a.Determinant() << endl;
  a.Print();
  cout << endl << endl;

  return;
}

double leptonic_fitter_algebraic::cofactor( const TMatrixD& mat, int i, int j )
{
  assert( mat.GetNcols() == 3 && mat.GetNrows() == 3 );
  assert( i >= 0 && i < 3 && j >= 0 && j < 3 );

  int i0 = ( i      ) ? 0 : 1;
  int i1 = ( i == 2 ) ? 1 : 2;
  int j0 = ( j      ) ? 0 : 1;
  int j1 = ( j == 2 ) ? 1 : 2;

  int sign = ( ( i+j ) % 2 ) ? -1 : 1;

  return sign * ( mat[i0][j0] * mat[i1][j1] - mat[i0][j1] * mat[i1][j0] );
}

void leptonic_fitter_algebraic::real_eigenvalues( const TMatrixD& mat )
{
  Int_t ncol = mat.GetNcols();
  assert( ncol < 10 && mat.GetNrows() == ncol );

  _real_eigs.clear();

  TMatrixDEigen eigen( mat );
  const TVectorD& e_res = eigen.GetEigenValuesRe();
  const TVectorD& e_ims = eigen.GetEigenValuesIm();
  for( int ir=0; ir < ncol; ++ir ) {
    if( e_ims[ ir ] == 0 ) _real_eigs.push_back( e_res[ ir ] );
  }
}

int LoadMatrix(const TString &fname, TMatrixD &M, TMatrixD &Merr, const TString &field, const TString &fieldErr) {
  TFile f(fname,"read");
  TMatrixD *Mptr=(TMatrixD*)f.Get(field);
  TMatrixD *MerrPtr=(TMatrixD*)f.Get(fieldErr);
  f.Close();
  if (!Mptr || !MerrPtr) {
    std::cout << "error in loading from <" << fname << ">\n";
    if (!Mptr) std::cout << " failed to load <" << field << ">\n";
    if (!MerrPtr) std::cout << " failed to load <" << fieldErr << ">\n";
    return 0;
  }
  if ((Mptr->GetNrows() != M.GetNrows()) ||
      (Mptr->GetNcols() != M.GetNcols()) ) {
    std::cout << "dim mismatch in <" << fname << ">\n";
    return 0;
  }
  M = *Mptr;
  Merr= *MerrPtr;
  delete Mptr;
  delete MerrPtr;
  return 1;
}

//==========================================================================
//
//  Dump element content of TMatrixD and TVectorD
//
//==========================================================================
void dumpElements(TMatrixD& a)
{
  cout << endl << endl;
  const int nrows = a.GetNrows();
  const int ncols = a.GetNcols();
  if(nrows==ncols)
    cout << "determinent = " << a.Determinant() << endl;
  cout << "------------------------------------------------------------------"
       << endl;
  for(int i=0; i< nrows; i++)
    {
      for(int j=0; j< ncols; j++)
	cout << a(i,j) << ", \t";    
	
      cout << endl;
    }
  cout << "------------------------------------------------------------------"
       << endl;
  cout << endl << endl;

  return;
}

//.........这里部分代码省略.........
      } // end loop over dielectrons

    } // end loop over events 
    delete infile;
    infile=0, eventTree=0;
  } // end loop over files
  delete gen;


  //finish reweighting of mass spectra

  
  // Find bin by bin corrections and the errors
  double tCentral, tErrNeg, tErrPos;
  for(int i=0; i<DYTools::nMassBins; i++){
    if ( DetCorrFactorDenominator(i) != 0 ){
      // This method does not take into account correlation between numerator
      // and denominator in calculation of errors. This is a flaw to be corrected
      // in the future.
      SimpleDivide( DetCorrFactorNumerator(i), DetCorrFactorDenominator(i), tCentral, tErrNeg, tErrPos);
      DetCorrFactor(i) = tCentral;
      DetCorrFactorErrPos(i) = tErrPos;
      DetCorrFactorErrNeg(i) = tErrNeg;
    } else {
      DetCorrFactor(i) = 0;
      DetCorrFactorErrPos(i) = 0;
      DetCorrFactorErrNeg(i) = 0;
    }
  }
  
  // Find response matrix, which is simply the normalized migration matrix
  for(int ifsr = 0; ifsr < DetMigration.GetNrows(); ifsr++){
    double nEventsInFsrMassBin = 0;
    for(int ireco = 0; ireco < DetMigration.GetNcols(); ireco++)
      nEventsInFsrMassBin += DetMigration(ifsr,ireco);
    // Now normalize each element and find errors
    for(int ireco = 0; ireco < DetMigration.GetNcols(); ireco++){
      tCentral = 0;
      tErrPos  = 0;
      tErrNeg  = 0;
      // Note: the event counts supposedly are dominated with events with weight "1"
      // coming from the primary MC sample, so the error is assumed Poissonian in
      // the call for efficiency-calculating function below.
      if( nEventsInFsrMassBin != 0 ){
	EfficiencyDivide(DetMigration(ifsr,ireco), nEventsInFsrMassBin, tCentral, tErrNeg, tErrPos);
      // BayesEfficiency does not seem to be working in newer ROOT versions, 
      // so it is replaced by simler method
//         BayesEfficiency( DetMigration(ifsr,ireco), nEventsInFsrMassBin, tCentral, tErrNeg, tErrPos);
      }
      DetResponse      (ifsr,ireco) = tCentral;
      DetResponseErrPos(ifsr,ireco) = tErrPos;
      DetResponseErrNeg(ifsr,ireco) = tErrNeg;
    }
  }

  // Find inverted response matrix
  TMatrixD DetInvertedResponse = DetResponse;
  Double_t det;
  DetInvertedResponse.Invert(&det);
  TMatrixD DetInvertedResponseErr(DetInvertedResponse.GetNrows(), DetInvertedResponse.GetNcols());
  calculateInvertedMatrixErrors(DetResponse, DetResponseErrPos, DetResponseErrNeg, DetInvertedResponseErr);

  // Package bin limits into TVectorD for storing in a file
  TVectorD BinLimitsArray(DYTools::nMassBins+1);
  for(int i=0; i<DYTools::nMassBins+1; i++)
    BinLimitsArray(i) = DYTools::massBinLimits[i];