C++ TMinuit::DefineParameter方法代码示例

void SetParameter(const unsigned int i, const unsigned int num_params, TMinuit& minuit){
  if(i==num_params-3){
    minuit.DefineParameter(i, "chi_b", 10.0, 1.0, 0.0, 0.0);
  }else if(i==num_params-2){
    minuit.DefineParameter(i, "chi_c", 10.0, 1.0, 0.0, 0.0);
  }else if(i==num_params-1){
    minuit.DefineParameter(i, "sig str", 0.5, 0.05, 0.0, 0.0);
  }else if(i<num_params-3){
    char name[128];
    switch(i%5){
    case 0:
      sprintf(name, "mu (%d)", i/5);
      break;
    case 1:
      sprintf(name, "As (%d)", i/5);
      break;
    case 2:
      sprintf(name, "Bs (%d)", i/5);
      break;
    case 3:
      sprintf(name, "Cs (%d)", i/5);
      break;
    case 4:
      sprintf(name, "Ds (%d)", i/5);
      break;
    }
    minuit.DefineParameter(i, name, 10.0, 1.0, 0.0, 0.0);
  }
}

vector<double> fit()
{
  vector<double> out(2*nbeta+1);
  
  TMinuit minu;
  //minu.SetPrintLevel(-1);
  minu.SetFCN(ch2);
  //set tolerance
  double tol=1e-16;
  int iflag;
  //minu.SetMaxIterations(10000000);
  //minu.mnexcm("SET ERR",&tol,1,iflag);
  
  for(int ib=0;ib<nbeta;ib++)
    {
      minu.DefineParameter(2*ib+0,"infi",1,0.001,0,0);
      minu.DefineParameter(2*ib+1,"coef",1,0.001,0,0);
    }
  minu.DefineParameter(2*nbeta,"M",0,0.001,0,0);
  
  minu.Migrad();
  
  double dum;
  for(int ib=0;ib<=2*nbeta;ib++) minu.GetParameter(ib,out[ib],dum);
  
  return out;
}

void fit(boot &A,boot &B,boot &C,bvec &X,bvec &Y)
{
  //copy X
  X_fit=new double[nens];
  for(int iens=0;iens<nens;iens++) X_fit[iens]=X[iens].med();
  Y_fit=new double[nens];
  err_Y_fit=new double[nens];
  
  TMinuit minu;
  minu.SetPrintLevel(-1);
  
  minu.DefineParameter(0,"A",0.0,0.0001,0,0);
  minu.DefineParameter(1,"B",0.0,0.0001,0,0);
  minu.DefineParameter(2,"C",0.0,0.0001,0,0);
  
  minu.SetFCN(chi2_wr);
  
  double C2;
  for(int iboot=0;iboot<nboot+1;iboot++)
    {
      if(iboot>0)
        minu.SetPrintLevel(-1);
      
      minu.DefineParameter(3,"a380",lat[0][iboot],0.0001,0,0);
      minu.DefineParameter(4,"a390",lat[1][iboot],0.0001,0,0);
      minu.DefineParameter(5,"a405",lat[2][iboot],0.0001,0,0);
      minu.DefineParameter(6,"a420",lat[3][iboot],0.0001,0,0);
      minu.FixParameter(3);
      minu.FixParameter(4);
      minu.FixParameter(5);
      minu.FixParameter(6);
      
      for(int iens=0;iens<nens;iens++)
        {
          Y_fit[iens]=Y.data[iens].data[iboot];
          err_Y_fit[iens]=Y.data[iens].err();
        }
      
      //minimize
      minu.Migrad();
      
      //get back parameters
      double dum;
      minu.GetParameter(0,A.data[iboot],dum);
      minu.GetParameter(1,B.data[iboot],dum);
      minu.GetParameter(2,C.data[iboot],dum);
      
      double lat_med[4]={lat[0].med(),lat[1].med(),lat[2].med(),lat[3].med()};
      if(iboot==0) C2=chi2(A.data[iboot],B[iboot],C[iboot],lat_med);
    }
  
  //calculate the chi2
  cout<<"A = ("<<A<<"), B=("<<B<<"), C=("<<C<<")"<<endl;
  cout<<"Chi2 = "<<C2<<" / "<<nens-3<<" = "<<C2/(nens-3)<<endl;
  
  delete[] X_fit;
  delete[] Y_fit;
  delete[] err_Y_fit;
}

// --------------------------------- //
void make_lin_fit(double & slope, double & d_slope, double & offset, double & d_offset){
    TMinuit min;
    min.SetPrintLevel(-1);
    //min.SetPrintLevel(0);
    int err = min.DefineParameter(0, "slope", slope, d_slope, 0.05, 1.0);
    assert(err==0);
    err = min.DefineParameter(1, "offset", offset, d_offset, 0.001, 0.2);
    assert(err==0);
    min.SetFCN(chi2_linear);
    min.mnmigr();
    min.GetParameter(0, slope, d_slope);
    min.GetParameter(1, offset, d_offset);
}

void fit(boot &A,boot &B,boot &C,boot &D)
{
  //copy ml
  ml_fit=new double[nens];
  for(int iens=0;iens<nens;iens++) ml_fit[iens]=ml[iens].med();
  
  //alloc dM2Pi and dM2K
  dM2Pi_fit=new double[nens];
  err_dM2Pi_fit=new double[nens];
  dM2K_fit=new double[nens];
  err_dM2K_fit=new double[nens];
  
  TMinuit minu;
  minu.SetPrintLevel(-1);
  
  int npars=4;
  minu.DefineParameter(0,"A",0.0,0.0001,0,0);
  minu.DefineParameter(1,"B",0.0,0.0001,0,0);
  minu.DefineParameter(2,"C",0.0,0.0001,0,0);
  minu.DefineParameter(3,"D",0.0,0.0001,0,0);
  minu.SetFCN(chi2_wr);
  
  double C2;
  for(int iboot=0;iboot<nboot+1;iboot++)
    {
      if(iboot>0) minu.SetPrintLevel(-1);
      
      minu.DefineParameter(4,"a380",lat[0][iboot],0.0001,0,0);
      minu.DefineParameter(5,"a390",lat[1][iboot],0.0001,0,0);
      minu.DefineParameter(6,"a405",lat[2][iboot],0.0001,0,0);
      minu.DefineParameter(7,"a420",lat[3][iboot],0.0001,0,0);
      minu.FixParameter(4);
      minu.FixParameter(5);
      minu.FixParameter(6);
      minu.FixParameter(7);
      
      for(int iens=0;iens<nens;iens++)
        {
          dM2Pi_fit[iens]=dM2Pi.data[iens].data[iboot];
          err_dM2Pi_fit[iens]=dM2Pi.data[iens].err();
          dM2K_fit[iens]=dM2K.data[iens].data[iboot];
          err_dM2K_fit[iens]=dM2K.data[iens].err();
        }
      
      //minimize
      minu.Migrad();
      
      //get back parameters
      double dum;
      minu.GetParameter(0,A.data[iboot],dum);
      minu.GetParameter(1,B.data[iboot],dum);
      minu.GetParameter(2,C.data[iboot],dum);
      minu.GetParameter(3,D.data[iboot],dum);
      
      double lat_med[4]={lat[0].med(),lat[1].med(),lat[2].med(),lat[3].med()};
      if(iboot==nboot) C2=chi2(A.data[iboot],B[iboot],C[iboot],D[iboot],lat_med,true);
    }
  
  //calculate the chi2
  cout<<"A=("<<A<<"), B=("<<B<<"), C=("<<C<<"), D=("<<D<<")"<<endl;
  cout<<"Chi2 = "<<C2<<" / "<<2*nens-npars<<" = "<<C2/(2*nens-npars)<<endl;
  
  delete[] ml_fit;
  delete[] dM2Pi_fit;
  delete[] err_dM2Pi_fit;
  delete[] dM2K_fit;
  delete[] err_dM2K_fit;
}

void parab_fit(double ch2,jack &a,jack &b,jack &c,double *X,jvec Y,const char *path,bool fl_parab=1)
{
  ofstream out(path);
  out<<"@type xydy\n";
  
  double MX=0;
  //copy X
  nc_fit=Y.nel;
  X_fit=new double[nc_fit];
  for(int iel=0;iel<nc_fit;iel++)
    if(use[iel])
    {
      X_fit[iel]=X[iel];
      out<<X[iel]<<" "<<Y[iel]<<endl;
      if(X[iel]>MX) MX=X[iel];
    }
  Y_fit=new double[nc_fit];
  err_Y_fit=new double[nc_fit];
  
  TMinuit minu;
  minu.SetPrintLevel(0);
  
  int npars=3;
  minu.DefineParameter(0,"A",0.0,0.0001,0,0);
  minu.DefineParameter(1,"B",0.0,0.0001,0,0);
  minu.DefineParameter(2,"C",0.0,0.0001,0,0);
  if(!fl_parab)
    {
      minu.FixParameter(0);
      npars--;
    }
  minu.SetFCN(chi2_wr);
  
  for(int ijack=0;ijack<=njack;ijack++)
    {
      minu.DefineParameter(2,"C",fixed_A[ijack],0.0001,0,0);
      minu.FixParameter(2);
      if(ijack>0)
        minu.SetPrintLevel(-1);
      
      for(int iel=0;iel<nc_fit;iel++)
        {
          Y_fit[iel]=Y.data[iel].data[ijack];
          err_Y_fit[iel]=Y[iel].err();
        }
      
      //minimize
      minu.Migrad();
      
      //get back parameters
      double dum;
      minu.GetParameter(0,a.data[ijack],dum);
      minu.GetParameter(1,b.data[ijack],dum);
      minu.GetParameter(2,c.data[ijack],dum);
    }
  
  out<<"&\[email protected] xy\n";
  for(double t=0;t<MX;t+=MX/100)
    out<<t<<" "<<parab(a[njack],b[njack],c[njack],t)<<endl;
  
  out.close();
}

void runAngleOpt() {

  gSystem->AddIncludePath("-I${ANITA_UTIL_INSTALL_DIR}/include");

  double startVal=0;               
  double stepSize=0.1;       
  double minVal=-2;           
  double maxVal=2;            
  Double_t p0 = 0;

  //Load libraries. Need to have ANITA_UTIL_INSTALL_DIR/lib and ROOTSYS/lib in the LD_LIBRARY_PATH               
  gSystem->Load("libfftw3.so");
  gSystem->Load("libMathMore.so");
  gSystem->Load("libPhysics.so");  
  gSystem->Load("libGeom.so");  
  gSystem->Load("libMinuit.so");  
  gSystem->Load("libRootFftwWrapper.so");         
  gSystem->Load("libAnitaEvent.so");      
  gSystem->Load("libAnitaCorrelator.so");

  AnitaGeomTool *fGeomTool = AnitaGeomTool::Instance();
  gSystem->CompileMacro("anglePlotterOpt.C","k");
  
  Double_t relDeltaOut=0;

  TMinuit *myMin = new TMinuit(1);
  myMin->SetObjectFit(anglePlotterOpt);
  myMin->SetFCN(iHateRoot);
  //setArray();

  for(int u = 0; u < 1; u++){

    int middle = 16;

      for(int y = 16; y <32; y++){

	int leftOpt, rightOpt;
	fGeomTool->getThetaPartners(middle,leftOpt,rightOpt); 


	myMin->DefineParameter(0, "antNum", middle, stepSize, minVal, maxVal);
	myMin->FixParameter(0);  
    
	myMin->DefineParameter(1, "deltaT", startVal, stepSize, minVal, maxVal);

  
	Double_t deltaT,deltaTErr;
  
	//*********MINUIT METHOD*******************
	myMin->SetPrintLevel(-1);
	myMin->Migrad();   
	myMin->GetParameter(1,deltaT,deltaTErr);

  

	setValue(rightOpt,deltaT);

	//	printArray();

	//	cout << middle << "  " << rightOpt << "  " << deltaT << endl;
  
	cout << "deltaTArrayMod[" << rightOpt << "] = " << deltaT << ";" << endl;

	middle = rightOpt;

      }

  }

  
  //   myMin->DeleteArrays();

  //   myMin->DeleteArrays();

  
}

void runTopRingOpt() {

  gSystem->AddIncludePath("-I${ANITA_UTIL_INSTALL_DIR}/include");

  double startVal=0;               
  double stepSize=0.1;       
  double minVal=-0.5;           
  double maxVal=0.5;            
  Double_t p0 = 0;

  //Load libraries. Need to have ANITA_UTIL_INSTALL_DIR/lib and ROOTSYS/lib in the LD_LIBRARY_PATH               
  gSystem->Load("libfftw3.so");
  gSystem->Load("libMathMore.so");
  gSystem->Load("libPhysics.so");  
  gSystem->Load("libGeom.so");  
  gSystem->Load("libMinuit.so");  
  gSystem->Load("libRootFftwWrapper.so");         
  gSystem->Load("libAnitaEvent.so");      
  gSystem->Load("libAnitaCorrelator.so");

  AnitaGeomTool *fGeomTool = AnitaGeomTool::Instance();
  gSystem->CompileMacro("topRingOpt.C","k");
  
  Double_t relDeltaOut=0;

  TMinuit *myMin = new TMinuit(150);
  myMin->SetObjectFit(topRingOpt);
  myMin->SetFCN(iHateRoot);
  //setArray();

  double startValDeltaT[16] ={0};
  double startValR[16] ={0};
  double startValPhi[16] ={0};

 //    startValDeltaT[0] =  -0.0519515;     startValR[0] =   -0.0101463;         startValPhi[0] =      -0.00473836 ;   
//     startValDeltaT[1] =     -0.0597062;   startValR[1] =      -0.02577;       startValPhi[1] =   0.00864501 ; 
//     startValDeltaT[2] =     -0.081435;     startValR[2] =     -0.000224044;    startValPhi[2] =   -0.000630649;
//     startValDeltaT[3] =     0.0118873;     startValR[3] =      0.019945;       startValPhi[3] =    0.014016;   
//     startValDeltaT[4] =      0.017917;      startValR[4] =     -0.00297559;    startValPhi[4] =    0.0224936 ;
//     startValDeltaT[5] =      0.0377119;     startValR[5] =     -0.014872;       startValPhi[5] =   0.0163349; 
//     startValDeltaT[6] =     -0.0426158;     startValR[6] =    -0.0562555;      startValPhi[6] =    0.0220065 ; 
//     startValDeltaT[7] =     -0.0221673;     startValR[7] =    -0.034104 ;      startValPhi[7] =    0.0158545 ; 
//     startValDeltaT[8] =      0.0263739;     startValR[8] =     0.00248804;      startValPhi[8] =    0.013246 ; 
//     startValDeltaT[9] =      -0.0938419;    startValR[9] =     -0.00344703;     startValPhi[9] =    -0.00718616; 
//     startValDeltaT[10] =      0.145264;      startValR[10] =     -0.0121874 ;     startValPhi[10] =     0.0156988 ;  
//     startValDeltaT[11] =      0.118105;      startValR[11] =     -0.0337033 ;      startValPhi[11] =   -0.00324182 ;
//     startValDeltaT[12] =      0.321805;      startValR[12] =     0.0134362  ;      startValPhi[12] =   -0.00190277 ;
//     startValDeltaT[13] =      0.0197693;     startValR[13] =     -0.000656063;      startValPhi[13] =   -0.0162318  ;
//     startValDeltaT[14] =      -0.115263;     startValR[14] =     0.0495637   ;      startValPhi[14] =   -0.0198119 ; 
//     startValDeltaT[15] =      -0.255707;   startValR[15] =     0.00189892  ;     startValPhi[15] =     0.0383932  ;
    


 for(int y = 0; y <16; y++){

   ofstream newfile("newSimonNumbers.txt");

   char name[30];
   sprintf(name,"r%d",y);
	myMin->DefineParameter(y, name, startValR[y], stepSize, minVal, maxVal);
   sprintf(name,"z%d",y);
	myMin->DefineParameter(y+16, name, startValDeltaT[y], stepSize, minVal, maxVal);
	   sprintf(name,"phi%d",y);
	myMin->DefineParameter(y+32, name, startValPhi[y], stepSize, minVal, maxVal);
 }


  
	Double_t deltaR[32],deltaRErr[32];
	Double_t deltaZ[32],deltaZErr[32];
	Double_t deltaPhi[32],deltaPhiErr[32];
  
	//*********MINUIT METHOD*******************
	myMin->SetPrintLevel(-1);
	myMin->Migrad();   

	for(int u = 0; u <16; u++){
	  myMin->GetParameter(u,deltaR[u],deltaRErr[u]);
	  //cout << "deltaR[" << u << "] = " << deltaR[u] ;
	
	
	  myMin->GetParameter(u+16,deltaZ[u],deltaZErr[u]);
	  //cout << " deltaZ[" << u << "] = " << deltaZ[u] ;
	
	  myMin->GetParameter(u+32,deltaPhi[u],deltaPhiErr[u]);
	  //cout << " deltaPhi[" << u << "] = " << deltaPhi[u] << ";" << endl;

	  newfile << u << "  " << deltaZ[u]<< "  " << deltaR[u]<< "  " << deltaPhi[u]<< "  " << 0 << endl;


	}


  
}

void AverageMW()
{
   // include NuTeV?
   if( inclNuTeV ) nVals = 5;
   else            nVals = 4;

   // number of parameters
   const int npar = 1;
   
   // create minuit pointer
   TMinuit *minuit = new TMinuit( npar );
   
   // set FCN function
   minuit->SetFCN( &fcn );
 
   // -1=quiet, 0=normal, 1=verbose
   minuit->SetPrintLevel( 0 );

   // start values
   Double_t startVal = 80;
   Double_t fitStep  = 0.1;
   Double_t limitMin = 70;
   Double_t limitMax = 90;
   minuit->DefineParameter( 0, "MW", startVal, fitStep, limitMin, limitMax);

   // minimize with MIGRAD
   Int_t ierr = 0;
   Double_t args[2];
   args[0] = 0; // maximum function calls (0=default)
   args[1] = 0; // tolerance at minimum (0=default)
   minuit->mnexcm( "MIGrad", args, 2, ierr );
   if ( ierr != 0 )
      cout << "Warning: Maybe fit didn't converge!" << endl;

   // fit results
   Double_t fitval_MW[1], fiterr_MW[1]; 
   minuit->GetParameter( 0, fitval_MW[0], fiterr_MW[0] );

   cout << "\n\n*************************************************" << endl;
   cout << "          chi2Min = " << chi2Min  << endl;
   cout << "          n_dof   = " << nVals-1  << endl;
   cout << "          p-value = " << TMath::Prob(chi2Min, nVals-1)  << endl;
   cout << "          MW      = " << fitval_MW[0] << " +- " << fiterr_MW[0] << endl;
   cout << "*************************************************" << endl;

   // make plot
   gStyle->SetOptTitle(0);
   gStyle->SetOptStat(0);
   gStyle->SetOptFit(0);
   TCanvas* c = new TCanvas( "c", "Mean of W boson mass", 0, 0, 850, 500 );
   c->SetGridy();
   TH2F *frame = new TH2F("frame", "Mean of W boson mass", 1, 80, 80.6, 5, 0, 5. );
   frame->SetLineColor(0);
   frame->SetTickLength(0,"Y");
   frame->SetXTitle("M_{W} [GeV]");  
   frame->GetYaxis()->CenterLabels( 1 );
   frame->GetYaxis()->SetNdivisions( frame->GetNbinsY()+10, 1 );
   frame->GetYaxis()->SetBinLabel( 1, "" ); // no labels
   frame->Draw();
   
   // make labels
   MakeLabels(frame);
   
   // draw fit value
   DrawFitValue( frame, fitval_MW[0], fiterr_MW[0] );
      
   // draw measurements
   for( Int_t i=0; i<5; i++ ){
      DrawMeasurement( frame, 4-i, val_MW[i], err_MW[i] ); 
   }

   // Redraw axis
   frame->Draw("sameaxis");
   
}

void runQuickOptAllAntFixedZ() {
  
  gSystem->AddIncludePath("-I${ANITA_UTIL_INSTALL_DIR}/include");
  
  double startVal=0;               
  double stepSize=0.05;       
  double minVal=-0.5;           
  double maxVal=0.5;            
  Double_t p0 = 0;

  //Load libraries. Need to have ANITA_UTIL_INSTALL_DIR/lib and ROOTSYS/lib in the LD_LIBRARY_PATH               
  gSystem->Load("libfftw3.so");
  gSystem->Load("libMathMore.so");
  gSystem->Load("libPhysics.so");  
  gSystem->Load("libGeom.so");  
  gSystem->Load("libMinuit.so");  
  gSystem->Load("libRootFftwWrapper.so");         
  gSystem->Load("libAnitaEvent.so");      
  gSystem->Load("libAnitaCorrelator.so");

  AnitaGeomTool *fGeomTool = AnitaGeomTool::Instance();
  gSystem->CompileMacro("quickOptAllAntFixedZ.C","k");

  fillArrays(eventNumberIndex, thetaWaveIndex, phiWaveIndex, antIndex1, antIndex2, maxCorrTimeIndex, adjacent);
  
  Double_t relDeltaOut=0;

  TMinuit *myMin = new TMinuit(192);
  myMin->SetObjectFit(quickOptAllAntFixedZ);
  myMin->SetFCN(iHateRoot);
//   myMin->SetMaxIterations(2);
  //myMin->SetErrorDef(1000);

//   int ierflg;
//   double eps[1] = {2.};
//   myMin->mnexcm("SET EPS", eps, 1, ierflg);
  //setArray();

  double startValR[MAX_ANTENNAS] ={0};
  double startValZ[MAX_ANTENNAS] ={0};
  double startValPhi[MAX_ANTENNAS] ={0};
  double startValCableDelays[MAX_ANTENNAS] ={0};
  
  
  for(int y = 0; y <MAX_ANTENNAS; y++){
    
    char name[30];
    sprintf(name,"r%d",y);
    myMin->DefineParameter(y, name, startValR[y], stepSize, -0.1, 0.1);
    sprintf(name,"z%d",y);
    myMin->DefineParameter(y+MAX_ANTENNAS, name, startValZ[y], stepSize, -0.1, 0.1);
    myMin->FixParameter(y+MAX_ANTENNAS);
    sprintf(name,"phi%d",y);
    myMin->DefineParameter(y+MAX_ANTENNAS*2, name, startValPhi[y], stepSize, -0.1, 0.1);
    sprintf(name,"cable%d",y);
    myMin->DefineParameter(y+MAX_ANTENNAS*3, name, startValCableDelays[y], stepSize, -0.5, 0.5);
  }
  
  
  
  Double_t deltaR[MAX_ANTENNAS],deltaRErr[MAX_ANTENNAS];
  Double_t deltaZ[MAX_ANTENNAS],deltaZErr[MAX_ANTENNAS];
  Double_t deltaPhi[MAX_ANTENNAS],deltaPhiErr[MAX_ANTENNAS];
  Double_t deltaCableDelays[MAX_ANTENNAS],deltaCableDelaysErr[MAX_ANTENNAS];
  
  //*********MINUIT METHOD*******************
//   myMin->SetPrintLevel(-1);
  myMin->Migrad();
  //  int error_flag;
  //  myMin->mnexcm("MINOS",0,0,error_flag);  

  ofstream newfile("newLindaNumbersCableDelays_ALLANTENNAS_RMS_GRAD_FixedZ.txt");
  for(int u = 0; u <MAX_ANTENNAS; u++){
    myMin->GetParameter(u,deltaR[u],deltaRErr[u]);
    cout << "deltaR[" << u << "] = " << deltaR[u] << " +/- " << deltaRErr[u] << endl;
	
    
    myMin->GetParameter(u+MAX_ANTENNAS,deltaZ[u],deltaZErr[u]);
    cout << " deltaZ[" << u << "] = " << deltaZ[u] << " +/- " << deltaZErr[u] << endl;
    
    myMin->GetParameter(u+MAX_ANTENNAS*2,deltaPhi[u],deltaPhiErr[u]);
    cout << " deltaPhi[" << u << "] = " << deltaPhi[u] << " +/- " << deltaPhiErr[u] << endl;

    myMin->GetParameter(u+MAX_ANTENNAS*3,deltaCableDelays[u],deltaCableDelaysErr[u]);
    cout << " deltaCableDelays[" << u << "] = " << deltaCableDelays[u]  << " +/- " << deltaCableDelaysErr[u] << endl;
    
    newfile << u << "  " << deltaR[u]<< "  " << deltaZ[u]<< "  " << deltaPhi[u]<< "  " << deltaCableDelays[u] << endl;


  }
  
  
  cout << "Easy table" << endl;
  for(int u = 0; u <MAX_ANTENNAS; u++)  cout << u << " & " << deltaR[u]<< " & " << deltaZ[u]<< " & " << deltaPhi[u]<< " & " << deltaCableDelays[u] << endl;

  
}

/*DoMunuit{{{*/
void DoMinuit(double *par)
{
	cout << "      *************************************************" << endl; 
	cout << "      *          Minimization Fitting for             *" << endl;
	cout << "      *              F(y) Function                    *" << endl;
	cout << "      *             Z. Ye 09/20/2010                  *" << endl;
	cout << "      *************************************************" << endl; 
	cout << endl;
	gSystem->Load("libMinuit");
	static const Int_t iParaNum=5;

	TMinuit *pMinuit = new TMinuit(iParaNum);  //initialize TMinuit with a maximum of iParaNum params
	pMinuit->SetFCN(myFcn);

	Double_t arglist[10];
	Int_t ierflg = 0;

	/// set print level
	/*  SET PRIntout  <level>
		Sets the print level, determining how much output will be
		produced. Allowed values and their meanings are displayed
		after a SHOw PRInt command, and are currently <level>=:
		[-1]  no output except from SHOW commands
		[0]  minimum output
		[1]  default value, normal output
		[2]  additional output giving intermediate results.
		[3]  maximum output, showing progress of minimizations. */
	arglist[0] = 1;
	pMinuit->mnexcm("SET PRIntout",arglist,1,ierflg);

	/*
	   SET NOWarnings
	   Supresses Minuit warning messages.
	   SET WARnings
	   Instructs Minuit to output warning messages when suspicious
	   conditions arise which may indicate unreliable results.
	   This is the default.
	   */
	arglist[0] = 1;
	pMinuit->mnexcm("SET NOWarnings",arglist,0,ierflg);

	//Set Name of Parameters
	TString namePar[iParaNum]={"f0","bigB","alpha","a","b"};

	// Double_t inPar[iParaNum]={Target->f0,   
	// 			    Target->B,   
	// 			    Target->alpha,
	// 			    Target->a,
	// 			    Target->b};	
	Double_t inPar[iParaNum]={par[0],par[1],par[2],par[3],par[4]};   

	//Set Stepsize of Parameters
	Double_t step[iParaNum]={ 0.00000001,0.00000001,0.00001,0.00000001,0.00000001};

	//Set Min of Parameters, value==0 for No-Bound
	// Double_t minVal[iParaNum]={ 0.0001,0.0001,100.0,0.0001,0.0001};
	Double_t minVal[iParaNum]={ 0.0,0.0,0.0,0.0,0.0};

	//Set Max of Parameters, value==0 for No-Bound
	// Double_t maxVal[iParaNum]={ 0.01,0.01,300.0,0.01,0.01};
	Double_t maxVal[iParaNum]={ 0.0,0.0,0.0,0.0,0.0};

	//Initializing Parameters
	for(int ii=0;ii<iParaNum;ii++)
	{
		pMinuit->DefineParameter(ii,namePar[ii], 
				inPar[ii], step[ii], minVal[ii], maxVal[ii]);
	}

	//Fix parameters
	//arglist[0] = 3;   pMinuit->mnexcm("FIX ", arglist ,1,ierflg);

	arglist[0] = 1;
	pMinuit->mnexcm("SET ERR", arglist,1,ierflg);

	/*
	//SET LIMits  0->[parno]  1->[lolim]  2->[uplim]
	//	  arglist[0] = 0;	   //this means set all parameters with the same limit
	//	  arglist[0] = 1;	   //this means set 1st parameters with the specified limit
	//	  arglist[0] = 2;	   //this means set 2nd parameters with the sepecified limit
	arglist[0] = 0;
	arglist[1] = 0.;
	arglist[2] = 0.;
	pMinuit->mnexcm("SET LIMits",arglist,3,ierflg);
	*/

	//Set Iteration Times and Use Migrad to start minimization
	arglist[0] = 500;
	arglist[1] = 0.01;
	pMinuit->mnexcm("MIGRAD", arglist ,2,ierflg);

	//Get Results

	for(int ii = 0; ii < iParaNum; ii++)
	{
		pMinuit->GetParameter(ii,outPar[ii],err[ii]);    
		par[ii] = outPar[ii];
	}
	//Put the results into a file
//.........这里部分代码省略.........

void runQuickOptAllAntStepsLDBHPOL() {
  
  gSystem->AddIncludePath("-I${ANITA_UTIL_INSTALL_DIR}/include");
  
  double startVal=0;               
  double stepSize=0.01;       
  double minVal=-0.5;           
  double maxVal=0.5;            
  Double_t p0 = 0;

  //Load libraries. Need to have ANITA_UTIL_INSTALL_DIR/lib and ROOTSYS/lib in the LD_LIBRARY_PATH               
  gSystem->Load("libfftw3.so");
  gSystem->Load("libMathMore.so");
  gSystem->Load("libPhysics.so");  
  gSystem->Load("libGeom.so");  
  gSystem->Load("libMinuit.so");  
  gSystem->Load("libRootFftwWrapper.so");         
  gSystem->Load("libAnitaEvent.so");      
  gSystem->Load("libAnitaCorrelator.so");

  AnitaGeomTool *fGeomTool = AnitaGeomTool::Instance();
  fGeomTool->useKurtAnita3Numbers(1);
  gSystem->CompileMacro("quickOptAllAntStepsLDBHPOL.C","k");

  fillArrays(eventNumberIndex, thetaWaveIndex, phiWaveIndex, antIndex1, antIndex2, maxCorrTimeIndex, adjacent);
  
  Double_t relDeltaOut=0;

  TMinuit *myMin = new TMinuit(192);
  myMin->SetObjectFit(quickOptAllAntStepsLDBHPOL);
  myMin->SetFCN(iHateRoot);
 //   myMin->SetMaxIterations(2);
  //myMin->SetErrorDef(1000);

//   int ierflg;
//   double eps[1] = {2.};
//   myMin->mnexcm("SET EPS", eps, 1, ierflg);
  //setArray();

  double minValCableDelays[MAX_ANTENNAS] ={0};
  double maxValCableDelays[MAX_ANTENNAS] ={0};

  Double_t deltaR[MAX_ANTENNAS] = {0};
  Double_t deltaRErr[MAX_ANTENNAS] = {0};
  Double_t deltaZ[MAX_ANTENNAS] = {0};
  Double_t deltaZErr[MAX_ANTENNAS] = {0};
  Double_t deltaPhi[MAX_ANTENNAS] = {0};
  Double_t deltaPhiErr[MAX_ANTENNAS] = {0};
  Double_t deltaCableDelays[MAX_ANTENNAS] = {0};
  Double_t deltaCableDelaysErr[MAX_ANTENNAS] = {0};

  //  Double_t globalDeltaz = 0.242575;
  Double_t globalDeltaz = 0.03;
  
  //Min global   deltaZ[] = -0.242575 +/- 0.00198264

  for(int y = 0; y <MAX_ANTENNAS; y++){
 
    deltaR[y] = 0;//-0.05;
    deltaCableDelays[y] = 0. - (y==47)*0.35;
    minValCableDelays[y] = deltaCableDelays[y] - 0.15;
    maxValCableDelays[y] = deltaCableDelays[y] + 0.15;
    deltaZ[y] = 0;//(y<16)*(-globalDeltaz) + (y>15)*globalDeltaz;

    char name[30];
    sprintf(name,"r%d",y);
    myMin->DefineParameter(y, name, deltaR[y], stepSize, -0.3, 0.3);
    sprintf(name,"z%d",y);
    myMin->DefineParameter(y+MAX_ANTENNAS, name, deltaZ[y], stepSize, -0.3, 0.3);
    sprintf(name,"phi%d",y);
    myMin->DefineParameter(y+MAX_ANTENNAS*2, name, deltaPhi[y], stepSize, -0.3, 0.3);
    sprintf(name,"cable%d",y);
    myMin->DefineParameter(y+MAX_ANTENNAS*3, name, deltaCableDelays[y], stepSize, minValCableDelays[y], maxValCableDelays[y]);
  }



  for(int y = 0; y <MAX_ANTENNAS; y++){
    myMin->FixParameter(y); // fixed R 
    myMin->FixParameter(y+MAX_ANTENNAS); // fixed Z
    myMin->FixParameter(y+MAX_ANTENNAS*2); // fixed phi
    //    myMin->FixParameter(y+MAX_ANTENNAS*3); // fixed t

  }

  myMin->FixParameter(MAX_ANTENNAS*3); // fixed t0

  //*********MINUIT METHOD*******************
//   myMin->SetPrintLevel(-1);
  myMin->Migrad();
  //  int error_flag;
  //  myMin->mnexcm("MINOS",0,0,error_flag);  



  for(int u = 0; u <MAX_ANTENNAS; u++){
    myMin->GetParameter(u+MAX_ANTENNAS*3,deltaCableDelays[u],deltaCableDelaysErr[u]);
    cout << " deltaCableDelays[" << u << "] = " << deltaCableDelays[u]  << " +/- " << deltaCableDelaysErr[u] << endl;

    //    myMin->GetParameter(u+MAX_ANTENNAS,deltaZ[u],deltaZErr[u]);
//.........这里部分代码省略.........

void fit(boot &A,boot &B,boot &C,boot &D,bvec &X,bvec &Y)
{
  //copy X
  X_fit=new double[nens];
  for(int iens=0;iens<nens;iens++) X_fit[iens]=X[iens].med();
  Y_fit=new double[nens];
  err_Y_fit=new double[nens];
  
  TMinuit minu;
  minu.SetPrintLevel(-1);
  
  int npars=4;
  minu.DefineParameter(0,"A",0.0,0.0001,0,0);
  minu.DefineParameter(1,"B",0.0,0.0001,0,0);
  minu.DefineParameter(2,"C",0.0,0.0001,0,0);
  minu.DefineParameter(3,"D",0.0,0.0001,0,0);
  if(!include_a4)
    {
      minu.FixParameter(3);
      npars--;
    }
  if(!include_ml_term)
    {
      minu.FixParameter(1);
      npars--;
    }
  minu.SetFCN(chi2_wr);
  
  double C2;
  for(int iboot=0;iboot<nboot+1;iboot++)
    {
      if(iboot>0)
        minu.SetPrintLevel(-1);
      
      minu.DefineParameter(4,"a380",lat[0][iboot],0.0001,0,0);
      minu.DefineParameter(5,"a390",lat[1][iboot],0.0001,0,0);
      minu.DefineParameter(6,"a405",lat[2][iboot],0.0001,0,0);
      minu.DefineParameter(7,"a420",lat[3][iboot],0.0001,0,0);
      minu.FixParameter(4);
      minu.FixParameter(5);
      minu.FixParameter(6);
      minu.FixParameter(7);
      
      for(int iens=0;iens<nens;iens++)
        {
          Y_fit[iens]=Y.data[iens].data[iboot];
          err_Y_fit[iens]=Y.data[iens].err();
        }
      
      //minimize
      minu.Migrad();
      
      //get back parameters
      double dum;
      minu.GetParameter(0,A.data[iboot],dum);
      minu.GetParameter(1,B.data[iboot],dum);
      minu.GetParameter(2,C.data[iboot],dum);
      minu.GetParameter(3,D.data[iboot],dum);
      
      double lat_med[4]={lat[0].med(),lat[1].med(),lat[2].med(),lat[3].med()};
      if(iboot==nboot)
        {
          contr_flag=1;
          C2=chi2(A.data[iboot],B[iboot],C[iboot],D[iboot],lat_med);
          contr_flag=0;
        }
    }
  
  int ninc_ens=0;
  for(int iens=0;iens<nens;iens++)
    if(ibeta[iens]!=0 || include_380) ninc_ens++;
  
  //calculate the chi2
  cout<<"A=("<<A<<"), B=("<<B<<"), C=("<<C<<"), D=("<<D<<")"<<endl;
  cout<<"Chi2 = "<<C2<<" / "<<ninc_ens-npars<<" = "<<C2/(ninc_ens-npars)<<endl;
  
  delete[] X_fit;
  delete[] Y_fit;
  delete[] err_Y_fit;
}