void plot_covariance(TMatrixD cov, string path, string prefix, string name, string legend_position = "top_left", bool detail = false)
{
//declaring the canvas
    if (detail) { cout << "Ploting " << name << endl; }
    TCanvas *c1 = new TCanvas("c1","Canvas",0,29,1450,870);
    gStyle->SetOptStat(0);
    gStyle->SetOptTitle(kFALSE);
    gStyle->SetPalette(1);
    gStyle->SetPaintTextFormat("4.2g");
    gPad->SetFillColor(0);
    gPad->SetBorderMode(0);
    gPad->SetBorderSize(2);
    gPad->SetLeftMargin(0.10);
    gPad->SetRightMargin(0.10);
    gPad->SetTopMargin(0.01);
    gPad->SetFrameBorderMode(0);
    gPad->SetLogz();

    cov.Draw("colz");
    cov.Draw("text same");


//setting the output files
   string fileout = prefix + name;
   string out_png = path + "png/" + fileout + ".png";
   string out_c = path + "c/" + fileout + ".C";
   string out_eps = path + "eps/" + fileout + ".eps";
    
//save the file and close the canvas
    c1->Print( out_png.c_str() );
    c1->Print( out_c.c_str() );
    c1->Print( out_eps.c_str() );
    c1->Close();

}

   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;  
   }

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;
}

TMatrix EUTelState::computePropagationJacobianFromLocalStateToNextLocalState(TVector3 positionEnd, TVector3 momentumEnd, float arcLength,float nextPlaneID) {
	streamlog_out(DEBUG2) << "-------------------------------EUTelState::computePropagationJacobianFromStateToThisZLocation()-------------------------BEGIN" << std::endl;
	if(arcLength == 0 or arcLength < 0 ){ 
		throw(lcio::Exception( Utility::outputColourString("The arc length is less than or equal to zero. ","RED"))); 
	}
	TMatrixD curvilinearJacobian = geo::gGeometry().getPropagationJacobianCurvilinear(arcLength,getOmega(), computeCartesianMomentum().Unit(),momentumEnd.Unit());
	streamlog_out(DEBUG0)<<"This is the curvilinear jacobian at sensor:" << std::endl; 
	streamlog_message( DEBUG0, curvilinearJacobian.Print();, std::endl; );

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;
}

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 ] );
  }
}

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 example() {
   ROOT::R::TRInterface &r = ROOT::R::TRInterface::Instance();
   // print R version
   r.Execute("print(version$version.string)");

   // compute standard deviation of 1000 vector normal numbers

   double std_dev_r  = r.Eval("sd(rnorm(10000))");
   std::vector<double> v = r.Eval("rnorm(10000)");
   double std_dev_root = TMath::StdDev(v.begin(),v.end());
   std::cout << "standard deviation from R    = " << std_dev_r << std::endl;
   std::cout << "standard deviation from ROOT = " <<  std_dev_root << std::endl;
   if (!TMath::AreEqualAbs(std_dev_r,std_dev_root,0.1))
      Error("ROOT-R-Example","Different std-dev found");

   // use << to execute the R command instead of Execute
   r << "mat<-matrix(c(1,2,3,4,5,6),2,3,byrow=TRUE)";
   TMatrixD m = r["mat"];
   std::array<double,6> a = r.Eval("seq(1:6)");
   TMatrixD m2(2,3,a.data());

   if  (!(m==m2)) {
      Error("ROOT-R-Example","Different matrix  found");
      m.Print();
      m2.Print();
   }

   // example on how to pass ROOT objects to R
   std::vector<double> v_root{1,2,3,4,5,6,7,8};
   r["v"] = v_root;
   r << "v2<-seq(1:8)";
   bool isEqual = r.Eval("prod(v==v2)");
   if (!isEqual) {
      Error("ROOT-R-Example","Different vector created");
      r << "print(v)";
      r << "print(v2)";
   }

   // example on how to pass functions to R

   r["gaus_pdf"] = ROOT::Math::normal_pdf;

   r << "y<-gaus_pdf(0,1,1)";
   double value_r = r["y"];
   double value_root = ROOT::Math::normal_pdf(0,1,1);
   std::cout << "Function gaussian(0,1,1) evaluated in  R    = " << value_r << std::endl;
   std::cout << "Function gaussian(0,1,1) evaluated in ROOT  = " <<  value_root << std::endl;
   if (value_r != value_root)
      Error("ROOT-R-Example","Different function value forund in r = %f and ROOT = %f", value_r, value_root);
}

void leptonic_fitter_algebraic::swap_x_y( const TMatrixD& in, TMatrixD& out )
{
  if( in.GetNoElements() != 9 ) 
    throw std::runtime_error( "ERROR! leptonic_fitter_algebraic::swap_x_y handles only 3 by 3 matrices");

  if( &in == &out )
    throw std::runtime_error( "ERROR! leptonic_fitter_algebraic::swap_x_y needs different in and out objects");

  static int remap[ 9 ]={ 4, 3, 5,   1, 0, 2,  7, 6, 8 };
  static double vec_out[9];

  const double *vec_in = in.GetMatrixArray();
  for( int ic=0;ic<9;++ic ) vec_out[ ic ] = vec_in[ remap[ ic ] ];
  out.Use( 3, 3, vec_out );
}

int GetMassProfile(const TMatrixD &m, int yIdx, TVectorD &v) {
  v=0;
  for (int i=0; i<m.GetNrows(); ++i) {
    v[i] = m(i,yIdx);
  }
  return 1;
}

void MultiGaus(const TVectorD& parMeans, const TMatrixDSym& covMatrix, TVectorD& genPars)
{

  TRandom3 rnd(0);

  int nPars = parMeans.GetNrows();
  if(nPars <= 0) {
    Error("MultiGaus", "Must have >0 pars");
    return;
  }
  if(covMatrix.GetNrows() != nPars) {
    Error("MultiGaus", "parMeans.GetNrows() != covMatrix.GetNrows()");
    return;
  }
 
  // Check that covMatrix is symmetric
  for(int iRow = 0; iRow < nPars; iRow++) {
    for(int iCol = iRow; iCol < nPars; iCol++) {
      if(covMatrix(iRow, iCol) != covMatrix(iCol, iRow)) {
        Error("MultiGaus", "malformed cov matrix at row %d, col %d", iRow, iCol);
        return;
      }
    }
  }

  genPars.ResizeTo(nPars);

  TMatrixDSymEigen eigenvariances(covMatrix);
  
  TMatrixD V = eigenvariances.GetEigenVectors();

  TVectorD rotParMeans = V * parMeans;

  for(int iPar = 0; iPar < nPars; iPar++) {
    double variance = eigenvariances.GetEigenValues()[iPar];
    // check for positive-definiteness of covMatrix
    if(variance < 0) {
      Error("MultiGaus", "Got a negative eigenvariance (%f) on iPar = %d", variance, iPar);
    }
    genPars[iPar] = rnd.Gaus(rotParMeans[iPar], sqrt(variance));
  }

  V.Invert();
  
  genPars = V * genPars;

}

//==========================================================================
//
//  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;
}

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;
}

// Open a FITS file and retrieve the first plane of the image array 
// as a TASImage object
void FITS_tutorial1()
{
   printf("\n\n--------------------------------\n");
   printf("WELCOME TO FITS tutorial #1 !!!!\n");
   printf("--------------------------------\n");
   printf("We're gonna open a FITS file that contains only the\n");
   printf("primary HDU, consisting on an image.\n");
   printf("The object you will see is a snapshot of the NGC7662 nebula,\n");
   printf("which was taken by the author on November 2009 in Barcelona (CATALONIA).\n\n");
      
   if (!gROOT->IsBatch()) {
      //printf("Press ENTER to start..."); getchar();
   }
   
   // Open primary HDU from file
   TFITSHDU *hdu = new TFITSHDU("sample1.fits");
   if (hdu == 0) {
      printf("ERROR: could not access the HDU\n"); return;
   }
   printf("File successfully open!\n");
   
   // Dump the HDUs within the FITS file
   // and also their metadata
   //printf("Press ENTER to see summary of all data stored in the file:"); getchar();
   
   hdu->Print("F+");
   
   printf("....................................\n");
   // Here we get the exposure time.
   //printf("Press ENTER to retrieve the exposure time from the HDU metadata..."); getchar();
   printf("Exposure time = %s\n", hdu->GetKeywordValue("EXPTIME").Data());

   
   // Read the primary array as a matrix,
   // selecting only layer 0.
   // This function may be useful to
   // do image processing.
   printf("....................................\n");
   printf("We can read the image as a matrix of values.\n");
   printf("This feature is useful to do image processing, e.g:\n");
   printf("histogram equalization, custom filtering, ...\n");
   //printf("Press ENTER to continue..."); getchar();
   
   TMatrixD *mat = hdu->ReadAsMatrix(0);
   mat->Print();
   delete mat;
   
   // Read the primary array as an image,
   // selecting only layer 0.
   printf("....................................\n");
   printf("Now the primary array will be read both as an image and as a histogram,\n");
   printf("and they will be shown in a canvas.\n");
   //printf("Press ENTER to continue..."); getchar();
   
   TASImage *im = hdu->ReadAsImage(0);
   
   // Read the primary array as a histogram.
   // Depending on array dimensions, returned
   // histogram will be 1D, 2D or 3D
   TH1 *hist = hdu->ReadAsHistogram();
   
   
   TCanvas *c = new TCanvas("c1", "FITS tutorial #1", 800, 300);
   c->Divide(2,1);
   c->cd(1);
   im->Draw();
   c->cd(2);
   hist->Draw("COL");

   // Clean up
   delete hdu;
}

//#############################################################################
int TrTrackA::StraightLineFit(){

    int _Nhit = _Hit.size();

  // Reset fit values
    Chi2 = FLT_MAX;
    Ndof = 2*_Nhit-4;
    Theta = 0.0;
    Phi = 0.0;
    U[0] = 0.0;
    U[1] = 0.0;
    U[2] = 1.0;
    P0[0] = 0.0;
    P0[1] = 0.0;
    P0[2] = 0.0;
    Rigidity = FLT_MAX;

  // Consistency check on the number of hits
    if (_Nhit<3) return -2;
    P0[2] = _Hit[0].Coo[2];

  // Get hit positions and uncertainties
  // Scale errors (we will use sigmas in microns)
    double hits[_Nhit][3];
    double sigma[_Nhit][3];
    for (int i=0;i<_Nhit;i++){
      for (int j=0;j<3;j++){
          hits[i][j] = _Hit[i].Coo[j];
          sigma[i][j] = 1.e4*_Hit[i].ECoo[j];
      }
    }

  // Lenghts
    double lenz[_Nhit];
    for (int i=0;i<_Nhit;i++) lenz[i] = hits[i][2] - P0[2];

  // F and G matrices
    const int idim = 4;
    double d[2*_Nhit][idim];
    for (int i=0;i<_Nhit;i++) {
      int ix = i;
      int iy = i+_Nhit;
      for (int j=0;j<idim;j++) { d[ix][j] = 0; d[iy][j] = 0;}
      d[ix][0] = 1.;
      d[iy][1] = 1.;
      d[ix][2] = lenz[i];
      d[iy][3] = lenz[i];
    }

  // F*S_x*x + G*S_y*y
    double dx[idim];
    for (int j=0;j<idim;j++) {
      dx[j] = 0.;
      for (int l=0;l<_Nhit;l++) {
        dx[j] += d[l][j]/sigma[l][0]/sigma[l][0]*hits[l][0];
        dx[j] += d[l+_Nhit][j]/sigma[l][1]/sigma[l][1]*hits[l][1];
      }
    }

  // (F*S_x*F + G*S_y*G)
    double Param[idim];
    double InvCov[idim][idim];
    for (int j=0;j<idim;j++) {
      for (int k=0;k<idim;k++) {
        InvCov[j][k] = 0.;
        for (int l=0;l<_Nhit;l++) {
          InvCov[j][k] += d[l][j]/sigma[l][0]/sigma[l][0]*d[l][k];
          InvCov[j][k] += d[l+_Nhit][j]/sigma[l][1]/sigma[l][1]*d[l+_Nhit][k];
        } 
      }
    }
      
  // (F*S_x*F + G*S_y*G)**{-1}
    double determ = 0.0;
    TMatrixD ParaCovariance = TMatrixD(idim,idim,(Double_t*)InvCov," ");
    ParaCovariance = ParaCovariance.Invert(&determ);
    if (determ<=0) return -1;

  // Solution
    for (int k=0;k<idim;k++) {
      Param[k] = 0.;
      for (int i=0;i<idim;i++) {
        Param[k] += ParaCovariance(k,i)*dx[i];
      }
    }

  // Chi2 (xl and yl in microns, since sigmas are in microns too)
    Chi2 = 0.;
    for (int l=0;l<_Nhit;l++) {
      double xl = hits[l][0]*1.e4;
      double yl = hits[l][1]*1.e4;
      for (int k=0;k<idim;k++) {
        xl -= d[l][k]*Param[k]*1.e4;
        yl -= d[l+_Nhit][k]*Param[k]*1.e4;
      }
      Chi2 += xl/sigma[l][0]/sigma[l][0]*xl + yl/sigma[l][1]/sigma[l][1]*yl;
    }

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

//.........这里部分代码省略.........
  TH1F hMassIdx("h_massIdx","h_massIdx",locMassBinCount-1,massRangeEdges);
  //delete massRangeEdges;

  if (0) {
    printHisto(&hMassIdx);
    if (0) {
      for (int i=0; i<int(massHigh); ++i) {
	double m=i+0.5;
	std::cout << "i=" << i << ", m=" << m << ", idx=" << (hMassIdx.FindBin(m)-1) << "\n";
      }
    }
    return;
  }

  //std::cout << "massHigh=" << massHigh << "\n";
  //std::cout << "locMassBinCount=" << locMassBinCount << "\n";

  if (!inpMgr.Load(mc_input)) {
    return;
  }
  
  //--------------------------------------------------------------------------------------------------------------
  // Main code 
  //==============================================================================================================
  
  //  
  // Set up histograms
  //

  //vector<TH1D*> hZMassv;
  
  Double_t   nZv = 0;

  TMatrixD nEvents (locMassBinCount,DYTools::nYBinsMax);    // number of weigthed events
  TMatrixD nEventsDET (locMassBinCount,DYTools::nYBinsMax); // number of weighted events in the detector acceptance
  TMatrixD nEventsDETrecoPostIdx (locMassBinCount,DYTools::nYBinsMax);    // number of weigthed events
  TMatrixD w2Events (locMassBinCount,DYTools::nYBinsMax);
  TMatrixD w2EventsDET (locMassBinCount,DYTools::nYBinsMax);
  TMatrixD w2EventsDETrecoPostIdx (locMassBinCount,DYTools::nYBinsMax);
  Double_t nZpeak=0, w2Zpeak=0;
  double nZpeakDET=0, w2ZpeakDET=0;
  double nZpeakDETrecoPostIdx=0, w2ZpeakDETrecoPostIdx=0;

  nEvents = 0;
  w2Events = 0;
  nEventsDET = 0;
  w2EventsDET  = 0;
  nEventsDETrecoPostIdx = 0;
  w2EventsDETrecoPostIdx  = 0;

  //char hname[100];
  //for(UInt_t ifile = 0; ifile<fnamev.size(); ifile++) {
  //  sprintf(hname,"hZMass_%i",ifile); hZMassv.push_back(new TH1F(hname,"",500,0,500)); hZMassv[ifile]->Sumw2();
  //}

  // 
  // Read weights from a file
  //
  const bool useFewzWeights = useFEWZ;
  const bool cutZPT100 = true;
  FEWZ_t fewz(useFewzWeights,cutZPT100);
  if (useFewzWeights && !fewz.isInitialized()) {
    std::cout << "failed to prepare FEWZ correction\n";
    throw 2;
  }

//#############################################################################
double* TrTrackA::Prediction(int index){

    static double pred[7];
#pragma omp threadprivate (pred)
    int _Nhit = _Hit.size();

  // Consistency check on the number of hits
    if (_Nhit<4) return NULL;

  // Scale errors (we will use sigmas in microns)
    double hits[_Nhit][3];
    double sigma[_Nhit][3];
    for (int i=0;i<_Nhit;i++){
      for (int j=0;j<3;j++){
          hits[i][j] = _Hit[i].Coo[j];
          sigma[i][j] = 1.e4*_Hit[i].ECoo[j];
          if (i==index) sigma[i][j] *= 1.e2;
      }
    }

    if (!_PathIntExist) SetPathInt();


  // F and G matrices
    const int idim = 5;
    double d[2*_Nhit][idim];
    for (int i=0;i<_Nhit;i++) {
      int ix = i;
      int iy = i+_Nhit;
      for (int j=0;j<idim;j++) { d[ix][j] = 0; d[iy][j] = 0;}
      d[ix][0] = 1.;
      d[iy][1] = 1.;
      for (int k=0;k<=i;k++) {
        d[ix][2] += _PathLength[k];
        d[iy][3] += _PathLength[k];
        d[ix][4] += _PathLength[k]*_PathLength[k]*_PathIntegral_x[0][k];
        d[iy][4] += _PathLength[k]*_PathLength[k]*_PathIntegral_x[1][k];
        for (int l=k+1;l<=i;l++) {
            d[ix][4] += _PathLength[k]*_PathLength[l]*_PathIntegral_u[0][k];
            d[iy][4] += _PathLength[k]*_PathLength[l]*_PathIntegral_u[1][k];
        }
      }
    }

  // F*S_x*x + G*S_y*y
    double dx[idim];
    for (int j=0;j<idim;j++) {
      dx[j] = 0.;
      for (int l=0;l<_Nhit;l++) {
        dx[j] += d[l][j]/sigma[l][0]/sigma[l][0]*hits[l][0];
        dx[j] += d[l+_Nhit][j]/sigma[l][1]/sigma[l][1]*hits[l][1];
      }
    }

  // (F*S_x*F + G*S_y*G)

    double Param[idim];
    double InvCov[idim][idim];
    for (int j=0;j<idim;j++) {
      for (int k=0;k<idim;k++) {
        InvCov[j][k] = 0.;
        for (int l=0;l<_Nhit;l++) {
          InvCov[j][k] += d[l][j]/sigma[l][0]/sigma[l][0]*d[l][k];
          InvCov[j][k] += d[l+_Nhit][j]/sigma[l][1]/sigma[l][1]*d[l+_Nhit][k];
        } 
      }
    }
      
  // (F*S_x*F + G*S_y*G)**{-1}
    double determ = 0.0;
    TMatrixD ParaCovariance = TMatrixD(idim,idim,(Double_t*)InvCov," ");
    ParaCovariance = ParaCovariance.Invert(&determ);
    if (determ<=0) return NULL;

  // Solution
    //printf(">>>>>>>>>>>>> Cov AFTER:\n");
    for (int k=0;k<idim;k++) {
      Param[k] = 0.;
      for (int i=0;i<idim;i++) {
        Param[k] += ParaCovariance[k][i]*dx[i];
        //printf(" %f", ParaCovariance[k][i]);
      }
      //printf("\n");
    }

  // Chi2 (xl and yl in microns, since sigmas are in microns too)
    pred[0] = 0;
    pred[1] = 0;
    pred[2] = hits[index][2];
    pred[3] = acos(-sqrt(1-Param[2]*Param[2]-Param[3]*Param[3]));
    pred[4] = atan2(Param[3],Param[2]);
    pred[5] = 0.;
    pred[6] = 0.;
    for (int k=0;k<idim;k++) {
      pred[0] += d[index][k]*Param[k];
      pred[1] += d[index+_Nhit][k]*Param[k];
      for (int l=0;l<idim;l++) {
        pred[5] += d[index][k]*d[index][l]*ParaCovariance(k,l);
        pred[6] += d[index+_Nhit][k]*d[index+_Nhit][l]*ParaCovariance(k,l);
//.........这里部分代码省略.........

void solveLinear(Double_t eps = 1.e-12)
{
   cout << "Perform the fit  y = c0 + c1 * x in four different ways" << endl;

   const Int_t nrVar  = 2;
   const Int_t nrPnts = 4;

   Double_t ax[] = {0.0,1.0,2.0,3.0};
   Double_t ay[] = {1.4,1.5,3.7,4.1};
   Double_t ae[] = {0.5,0.2,1.0,0.5};

   // Make the vectors 'Use" the data : they are not copied, the vector data
   // pointer is just set appropriately

   TVectorD x; x.Use(nrPnts,ax);
   TVectorD y; y.Use(nrPnts,ay);
   TVectorD e; e.Use(nrPnts,ae);

   TMatrixD A(nrPnts,nrVar);
   TMatrixDColumn(A,0) = 1.0;
   TMatrixDColumn(A,1) = x;

   cout << " - 1. solve through Normal Equations" << endl;

   const TVectorD c_norm = NormalEqn(A,y,e);

   cout << " - 2. solve through SVD" << endl;
   // numerically  preferred method

   // first bring the weights in place
   TMatrixD Aw = A;
   TVectorD yw = y;
   for (Int_t irow = 0; irow < A.GetNrows(); irow++) {
      TMatrixDRow(Aw,irow) *= 1/e(irow);
      yw(irow) /= e(irow);
   }

   TDecompSVD svd(Aw);
   Bool_t ok;
   const TVectorD c_svd = svd.Solve(yw,ok);

   cout << " - 3. solve with pseudo inverse" << endl;

   const TMatrixD pseudo1  = svd.Invert();
   TVectorD c_pseudo1 = yw;
   c_pseudo1 *= pseudo1;

   cout << " - 4. solve with pseudo inverse, calculated brute force" << endl;

   TMatrixDSym AtA(TMatrixDSym::kAtA,Aw);
   const TMatrixD pseudo2 = AtA.Invert() * Aw.T();
   TVectorD c_pseudo2 = yw;
   c_pseudo2 *= pseudo2;

   cout << " - 5. Minuit through TGraph" << endl;

   TGraphErrors *gr = new TGraphErrors(nrPnts,ax,ay,0,ae);
   TF1 *f1 = new TF1("f1","pol1",0,5);
   gr->Fit("f1","Q");
   TVectorD c_graph(nrVar);
   c_graph(0) = f1->GetParameter(0);
   c_graph(1) = f1->GetParameter(1);

   // Check that all 4 answers are identical within a certain
   // tolerance . The 1e-12 is somewhat arbitrary . It turns out that
   // the TGraph fit is different by a few times 1e-13.

   Bool_t same = kTRUE;
   same &= VerifyVectorIdentity(c_norm,c_svd,0,eps);
   same &= VerifyVectorIdentity(c_norm,c_pseudo1,0,eps);
   same &= VerifyVectorIdentity(c_norm,c_pseudo2,0,eps);
   same &= VerifyVectorIdentity(c_norm,c_graph,0,eps);
   if (same)
      cout << " All solutions are the same within tolerance of " << eps << endl;
   else
      cout << " Some solutions differ more than the allowed tolerance of " << eps << endl;
}