C++ TF1::GetChisquare方法代码示例

void LongiProfile(int E, std::string runPeriod)
{
  char inputFile[200];
  sprintf(inputFile,"%s%d%s","/lyoserv/home/ilc/steen/resultRootFile/tb_data/DHCAL_",Run(E,runPeriod),".root");
  std::cout << "inputFile = " << inputFile << std::endl;
  PionTimeCalibration* timeCalib=new PionTimeCalibration(inputFile);
  timeCalib->ShowerBarycenterCut(runPeriod);
  timeCalib->Loop();
  
  std::cout<< timeCalib->longitudinalVec.size()<<std::endl;
  TF1 *func;
  char outputFile[200];
  sprintf(outputFile,"%s%d%s","/home/steen/timeCalib/LongiProfile/calib_",Run(E,runPeriod),".txt");
  fstream out;
  out.open(outputFile,ios::out);
  for(unsigned int k=0; k<48; k++){
    TH2D* hRadPro=new TH2D("hName","",100,0,10,4000,0,200);
    func=new TF1("func","pol2",0,10);
    for(unsigned int i=0; i<timeCalib->nhit1.size(); i++){
      hRadPro->Fill(timeCalib->time.at(i),timeCalib->longitudinalVec.at(i).layer[k]);
    }
    hRadPro->ProfileX()->Fit(func,"NQ");
    out << k 
	<< " " << func->GetParameter(0) << " " << func->GetParError(0) 
	<< " " << func->GetParameter(1) << " " << func->GetParError(1)
	<< " " << func->GetParameter(2) << " " << func->GetParError(2)
	<< " " << func->GetChisquare()/func->GetNDF() << std::endl;
    delete hRadPro;
  }
  out.close();
  delete func;
  delete timeCalib;
}

void Scan( TString tgt = "full", Int_t chan = 300, Bool_t plot = kTRUE,
		UInt_t rebin = 1)
{
	UInt_t bin;
	TString name;

	TH1D* h1;
	TH2D* h2;

	if ( tgt == "full") h2 = h2d_f;
	else if ( tgt == "empty") h2 = h2d_e;

	if ( ( chan >= 0) && ( chan <= 351)) {
		bin = h2->GetXaxis()->FindBin( chan);
		name = Form( "Channel %d", chan);
		h1 = h2->ProjectionY( name, bin, bin);
		name = Form( "E_{#gamma} = %5.1f #pm %3.1f MeV", tcd[chan].energy,
				tcd[chan].denergy);
	}
	else if ( chan == -1) {
		name = "All Channels";
		h1 = h2->ProjectionY( name, 0, 351);
	}
	else {
		cout << "Error: channels must be 0-351 or -1" << endl;
		break;
	}

	h1->SetTitle( name);

	TString opt;
	opt = "NQR+";
	if ( plot == kTRUE) {
		h1->Draw();
		opt = "QR+";
	}

	TF1 *f1 = new TF1( "f1", "pol1", 20, 80);
	TF1 *f2 = new TF1( "f2", "gaus", 95, 105);
	h1->Fit( "f1", "NQR");
	f1->GetParameters( &par[0]);
	h1->Fit( "f2", "NQR+");
	f2->GetParameters( &par[2]);
	TF1 *fall = new TF1( "fall", "pol1(0) + gaus(2)", 20, 180);
	fall->SetParameters( par);
	h1->Fit( "fall", opt);
	Int_t bins = h1->GetXaxis()->FindBin( 180) - h1->GetXaxis()->FindBin( 20);
	red_chisq = fall->GetChisquare()/(bins-5);

//	gStyle->SetOptStat( 0);
	gStyle->SetOptFit( 1);

	cout << chan;
	cout << " " << par[3];
	cout << " " << par[4];
	cout << " " << red_chisq;
	cout << endl;

}

/************************************************************
*  Fit ADC spectrum
************************************************************/
int FitAdcHist(TH1F *h, float data[NDATA], const int print ) {
  TF1 *fadc;
  TPaveText *adcPT;

  // Protect for null histograms
  if( h == 0 || h->GetEntries()<1000. ) return -1;
  printf("Fitting %s\n",h->GetName());

  gStyle->SetOptStat("emr"); 
  adcPT = new TPaveText(0.58,0.315,0.98,0.635,"NDC");

//  Refit ADC
  Double_t n    = h->GetEntries();
  Double_t rms  = h->GetRMS();
  Double_t mean = h->GetMean();
  Double_t ymin = mean-2*rms;
  if( ymin < 60. ) ymin=60.;
  Double_t ymax = mean+1.5*rms;
  fadc = new TF1("adc_fun",skewnormal,ymin,ymax,4);
  fadc->SetParameters(n,mean,rms,1.5);
  fadc->SetParLimits(1,mean-rms,mean+0.5*rms);
  fadc->SetParLimits(2,rms*0.5,rms*1.5);
  fadc->SetParLimits(3,1.,4.);
  h->Fit("adc_fun","r"); 
//  If skew is too high refit with fixed skew
  if( fadc->GetParameter(3) > 2.5 ) {
    printf("REFIT %s\n",h->GetName());
    fadc->SetParameter(3,1.6);
    fadc->SetParLimits(3,1.6,1.6);
    h->Fit("adc_fun","r") ;
  }
// if have enough hits use fit results  
  if( n > MINHITS ) {
    data[13] = fadc->GetParameter(1);
    data[14] = fadc->GetParError(1);
    data[15] = fadc->GetParameter(2);
    data[16] = fadc->GetParError(2);
    data[17] = fadc->GetParameter(3);
    data[18] = fadc->GetParError(3);
    double dof = fadc->GetNDF();
    if( dof > 0. ) data[19] = fadc->GetChisquare()/dof;
    adcPT->AddText(Form("Peak %.1lf #pm %.1lf",data[13],data[14]));
    adcPT->AddText(Form("Width %.1lf #pm %.1lf",data[15],data[16]));
    adcPT->AddText(Form("Skew %.2lf #pm %.2lf",	data[17],data[18]));
    if( dof > 0. ) adcPT->AddText(Form("#Chi^{2}/DoF     %.2lf",data[19]));
    adcPT->SetTextColor(2);
    adcPT->Draw();
// if have too few hits use fit stats instead of fit results  
  } else {
    data[12] = mean;
    data[13] = h->GetMeanError();
    data[14] = rms;
    data[15] = h->GetRMSError();
    data[16] = h->GetSkewness();
  }
  if( print ) c1->Print(Form("%s_%s.png",fname,h->GetName()));
  return 0;
}    //end FitAdcHist

	void FillFromFunction(TF1& infunc){
		norm = normscale * infunc.GetParameter(0);
		norm_err = normscale * infunc.GetParError(0);
		gamma = infunc.GetParameter(2);
		gamma_err = infunc.GetParError(2);
		ang = infunc.GetParameter(1);
		ang_err = infunc.GetParError(1);
		chi2 = infunc.GetChisquare();
		ndf = infunc.GetNDF();
	}

drawGraph( TGraph* gr, const Char_t* ytitle = 0 ) 
{
	gr->Draw( "AP*" );
	TH1* h = gr->GetHistogram( );

	std::cout << gr->GetMean( ) << std::endl;

	gr->Fit( "pol0", "Q", "SAME" );
	TF1* fit = gr->GetFunction( "pol0" );	
	Double_t chi2 	= fit->GetChisquare( );
	Double_t par0 	= fit->GetParameter( 0 );
	Double_t err 	= fit->GetParError( 0 );	

	TLatex* l = new TLatex;
	l->SetNDC( );
	l->DrawLatex( 0.85, 0.80, Form( "#splitline{Entries = %i}{RMS = %2.3f}",
		gr->GetN(), gr->GetMean(2)) );

	TLatex* l2 = new TLatex;
	l2->SetNDC( );
	l2->DrawLatex( 0.85, 0.7, Form("#splitline{fit = %2.3f#pm%2.3f}{#chi^{2} = %2.3f}", 
		par0, chi2, err) );
		
//	TLegend* leg = new TLegend( 0.8, 0.9, 0.95, 0.95 );
//	leg->SetFillColor( kWhite );
//	leg->SetBorderSize( 0 );
//	leg->AddEntry( "gr", "CPU+AFAR+MBED+ATWD+AMPS", "P" );
//	leg->Draw( );

	if( h != 0 ) 
	{
		h->GetXaxis()->SetTimeFormat("%m/%d %H:%M");
		h->GetXaxis()->SetNdivisions( 10, 10, 0 );
		h->GetXaxis()->SetTimeDisplay( 1 );
		h->GetXaxis()->SetTitle("Date / Time (UTC)");
		h->GetXaxis()->SetLabelSize( 0.06 );
		h->GetYaxis()->SetLabelSize(0.07);
		h->GetYaxis()->SetTitleSize(0.06);
		h->GetXaxis()->SetTitleSize(0.055);
		h->GetYaxis()->SetTitleOffset( 0.9 );
		h->GetXaxis()->SetTitleOffset( 0.9 );
		h->SetMarkerSize( 5 );

		if( ytitle != 0 ) 
		{
			h->GetYaxis()->SetTitle(ytitle);
      	}
	}
   
	if( gPad != 0 ) 
	{
		gPad->Update();
	}
}

TF1 *langaufit(TH1F *his, Double_t *fitrange, Double_t *startvalues, Double_t *parlimitslo, Double_t *parlimitshi, Double_t *fitparams, Double_t *fiterrors, Double_t *ChiSqr, Int_t *NDF)
{
   // Once again, here are the Landau * Gaussian parameters:
   //   par[0]=Width (scale) parameter of Landau density
   //   par[1]=Most Probable (MP, location) parameter of Landau density
   //   par[2]=Total area (integral -inf to inf, normalization constant)
   //   par[3]=Width (sigma) of convoluted Gaussian function
   //
   // Variables for langaufit call:
   //   his             histogram to fit
   //   fitrange[2]     lo and hi boundaries of fit range
   //   startvalues[4]  reasonable start values for the fit
   //   parlimitslo[4]  lower parameter limits
   //   parlimitshi[4]  upper parameter limits
   //   fitparams[4]    returns the final fit parameters
   //   fiterrors[4]    returns the final fit errors
   //   ChiSqr          returns the chi square
   //   NDF             returns ndf

   Int_t i;
   Char_t FunName[100];

   sprintf(FunName,"Fitfcn_%s",his->GetName());

   TF1 *ffitold = (TF1*)gROOT->GetListOfFunctions()->FindObject(FunName);
   if (ffitold) delete ffitold;

   TF1 *ffit = new TF1(FunName,langaufun,fitrange[0],fitrange[1],4);
   ffit->SetParameters(startvalues);
   ffit->SetParNames("Width","MP","Area","GSigma");

   for (i=0; i<4; i++) {
      ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]);
   }

   his->Fit(FunName,"RB0");   // fit within specified range, use ParLimits, do not plot

   ffit->GetParameters(fitparams);    // obtain fit parameters
   for (i=0; i<4; i++) {
      fiterrors[i] = ffit->GetParError(i);     // obtain fit parameter errors
   }
   ChiSqr[0] = ffit->GetChisquare();  // obtain chi^2
   NDF[0] = ffit->GetNDF();           // obtain ndf

   return (ffit);              // return fit function

}

void FillGoodRun(int ihar){
    float pi = acos(-1.0);
    TString str;
    TFile *fin;
    int nrun = GetTotalRun();
    if(nrun<0) exit(1);
     for(int icent=0;icent<ncent;icent++){
 //     for(int ihar=0;ihar<nhar;ihar++){
 //         if(ihar!=1) continue;
       for(int isub=0;isub<nsub;isub++){
            str = choosesub(isub);
            if(str=="ABORT") continue;
            for(int irun=0;irun<nrun;irun++){
      //std::cout<<icent<<" "<<ihar<<" "<<isub<<" "<<irun<<std::endl;
        //fin = TFile::Open(Form("/phenix/plhf/xuq/taxi/%s%s/%d/data/%s.root",dataset.Data(),pro.Data(),taxi,GetRun(irun).Data()));
         fin = TFile::Open(Form("/gpfs/mnt/gpfs02/phenix/plhf/plhf1/xuq/phenix/flow/Run16dAu/work/62GeV/output/%s",GetRun(irun).Data()));
        TH1F* hpsi = new TH1F("psi","psi",100,-pi,pi);
        for(int ibbcz=0;ibbcz<nbbcz;ibbcz++){
          TH1F* hpsitemp = (TH1F*)fin->Get(Form("psiFla_0_0_%d_%d_%d_%d",icent,ibbcz,ihar,isub));
          hpsi->Add(hpsitemp);
        }
        TF1 *fun = new TF1("fun","pol0",-pi,pi);
      if(hpsi->GetEntries()>1000){
	//hpsi->SetMarkerStyle(20);
	//hpsi->SetMarkerSize(0.6);
	//hpsi->SetMarkerColor(4);
	hpsi->SetMinimum(10);
	hpsi->Fit("fun","QR0");
	//float par=fun->GetParameter(0);
	//hpsi->SetMaximum(1.5*par);
	//hpsi->Draw();
	GoodRunFit[icent][ihar][isub][irun] = fun->GetChisquare()/fun->GetNDF();
        fin->Close();
      }
      else{
        GoodRunFit[icent][ihar][isub][irun] = -9999;
        fin->Close();
     }
     // GoodRunFit[icent][ihar][isub][irun] = 1.;
    }
       }
   //   }
     }
}

void Chi(int num, double* arr)
{
  //double chi[3] = arr;
  //double* pointer = χ
  //if(num > -1)
    {
      TF1* f = new TF1("f", "5*(TMath::Power(x,-4))",1,100);
      TH1D* h1 = new TH1D("h1","h1", 100, -0.5, 99.5);

      for(int i = 0; i < 1000000; i++)
	{
	  h1->Fill(f->GetRandom(1,100));
	  //cout << f.GetRandom(1,100) << endl;
	}
      
      //cout << __FILE__ << __LINE__<< endl;

      TF1* pwrLw = new TF1("pwrLw", "[0]*TMath::Power(x,[1])",1,1000);//changed to 0 parameter
      pwrLw->SetParameter(0,5);
      pwrLw->SetParameter(1,-4);
      //cout << pwrLw->GetParameter(0) << "\t" << pwrLw->GetParameter(1) << endl;

      //cout << __FILE__ << __LINE__<< endl;
      
      h1->Fit("pwrLw","0","",2,25);
      TF1* func = h1->GetFunction("pwrLw");
      double chi = func->GetChisquare();
      double dof = func->GetNDF();
      arr[0] = chi/dof;
      double par0 = func->GetParameter(0);
      double par1 = func->GetParameter(1);
      //cout << par0 << "\t" << par1 << endl;
      arr[1] = CalcChiSqr(h1,par0,par1,2,25);
      //chi[2] = CalcChiSqr(h1,5,-4,2,25);
      //cout << "Outputing values from Chi()" << endl;
      //cout << arr[0] << "\t" << arr[1] << endl;//"\t" << chi[2] << endl;
      //DoCleanUp(h1, f, pwrLw);
    }
    //arr = chi;
    //return chi;
}

TF1* getLandau(TH1* InputHisto, double* FitResults, double LowRange, double HighRange)
{
   FitResults[0]         = -0.5;  //MPV
   FitResults[1]         =  0;    //MPV error
   FitResults[2]         = -0.5;  //Width
   FitResults[3]         =  0;    //Width error
   FitResults[4]         = -0.5;  //Fit Chi2/NDF

   // perform fit with standard landau
   TF1* MyLandau = new TF1("MyLandau","landau",LowRange, HighRange);
   MyLandau->SetParameter(1,300);
   InputHisto->Fit("MyLandau","0QR WW");

   // MPV is parameter 1 (0=constant, 1=MPV, 2=Sigma)
   FitResults[0]         = MyLandau->GetParameter(1);  //MPV
   FitResults[1]         = MyLandau->GetParError(1);   //MPV error
   FitResults[2]         = MyLandau->GetParameter(2);  //Width
   FitResults[3]         = MyLandau->GetParError(2);   //Width error
   FitResults[4]         = MyLandau->GetChisquare() / MyLandau->GetNDF();  //Fit Chi2/NDF

   return MyLandau;
}

int main(int argc, char *argv[])
{
  int binCount = 100;
  if (argc > 2) {
    std::cout << "too many argument. exit." << std::endl;
    std::exit(1);
  } else if (argc == 2) {
    binCount = std::atoi(argv[1]);
    --argc;
  }
  TApplication app("app", &argc, argv);
  gStyle->SetOptStat(0);

  TH1D *hist = new TH1D("hist", "ns vs count", binCount, minNanoSec, maxNanoSec);

  char filename[256];
  GetFileName(filename, "../data/muon_lifetime/muon_lifetime.dat");
  std::ifstream ifs(filename);
  for (int i = 0; !ifs.eof(); ++i) {
    Double_t data;
    ifs >> data;
    if (data) {
      NanoSecWithError ns;
      ConvertTdcChannelToNanoSec(i%8, data, &ns);
      hist->Fill(ns.time);
    }
  }
  hist->Draw();

  TF1 *fit = new TF1("fit", "[0] * exp(- x / [1]) + [2]", 1000, 20000);
  fit->SetParameters(200, 2100, 40);
  hist->Fit(fit, "R+");
  std::cout << "chi^2/ndf: ";
  std::cout << fit->GetChisquare() / fit->GetNDF() << std::endl;

  app.Run();
  return 0;
}

void FillGoodRun(int icent,int ihar){
    float pi = acos(-1.0);
    TString str;
    TFile *fin;
    int nrun = GetTotalRun();
    if(nrun<0) exit(1);
       for(int isub=0;isub<nsub;isub++){
            str = choosesub(isub);
            if(str=="ABORT") continue;
            for(int irun=0;irun<nrun;irun++){
      std::cout<<icent<<" "<<ihar<<" "<<isub<<" "<<irun<<std::endl;
         fin = TFile::Open(Form("/store/user/qixu/flow/Run16dAu/62GeV/%s",GetRun(irun).Data()));
        TH1F* hpsi = new TH1F("psi","psi",100,-pi,pi);
        for(int ibbcz=0;ibbcz<nbbcz;ibbcz++){
          TH1F* hpsitemp = (TH1F*)fin->Get(Form("psiFla_%d_%d_%d_%d",icent,ibbcz,ihar,isub));
          hpsi->Add(hpsitemp);
        }
        TF1 *fun = new TF1("fun","pol0",-pi,pi);
      if(hpsi->GetEntries()>1000){
	//hpsi->SetMarkerStyle(20);
	//hpsi->SetMarkerSize(0.6);
	//hpsi->SetMarkerColor(4);
	hpsi->SetMinimum(10);
	hpsi->Fit("fun","QR0");
	//float par=fun->GetParameter(0);
	//hpsi->SetMaximum(1.5*par);
	//hpsi->Draw();
	GoodRunFit[icent][ihar][isub][irun] = fun->GetChisquare()/fun->GetNDF();
        fin->Close();
      }
      else{
        GoodRunFit[icent][ihar][isub][irun] = -9999;
        fin->Close();
     }
     // GoodRunFit[icent][ihar][isub][irun] = 1.;
    }
       }
}

float GoodRun(int icent, int ihar, int isub, int irun){
    float pi = acos(-1);
    TF1 *fun = new TF1("fun","pol0",-pi,pi);
    TString str;
    TFile *fin;

     ofstream fout;
        if(isub==1){
         str = "FVTX1S";
        }
        else if(isub==2){
         str = "FVTX2S";
        }
        else return -9999;
         fin = TFile::Open(Form("Run15pAu200MinBias/output_fvtxwithcntEP_%d.root",GetRun(irun)));
        TH1F* hpsi = new TH1F("psi","psi",100,-pi,pi);
        for(int ibbcz=0;ibbcz<nbbcz;ibbcz++){
          hpsitemp = (TH1F*)fin->Get(Form("psi_%d_%d_%d_%d",icent,ibbcz,ihar,isub));
          hpsi->Add(hpsitemp);
        }
      if(hpsi->GetEntries()>10000){
	hpsi->SetMarkerStyle(20);
	hpsi->SetMarkerSize(0.6);
	hpsi->SetMarkerColor(4);
	hpsi->SetMinimum(10);
	hpsi->Fit("fun","QR0");
	float par=fun->GetParameter(0);
	hpsi->SetMaximum(1.5*par);
	//hpsi->Draw();
        fin->Close();
	return fun->GetChisquare()/fun->GetNDF();
      }
      else{
        fin->Close();
        return -9999;
      }
    }

void Nhit5by5Calibration(int E, std::string runPeriod)
{
  char inputFile[200];
  sprintf(inputFile,"%s%d%s","/lyoserv/home/ilc/steen/resultRootFile/tb_data/DHCAL_",Run(E,runPeriod),".root");
  std::cout << "inputFile = " << inputFile << std::endl;
  PionTimeCalibration* timeCalib=new PionTimeCalibration(inputFile);
  timeCalib->ShowerBarycenterCut(runPeriod);
  timeCalib->Loop();
  
  std::cout<< timeCalib->nhit5by5.size()<<std::endl;
  TF1 *func;
  char outputFile[200];
  sprintf(outputFile,"%s%d%s","/home/steen/timeCalib/Nhit5by5/calib_",Run(E,runPeriod),".txt");
  fstream out;
  out.open(outputFile,ios::out);
  TH2D* hRadPro=new TH2D("hName","",100,0,10,600,0,600);
  func=new TF1("func","pol1",0,10);
  for(unsigned int i=0; i<timeCalib->nhit5by5.size(); i++){
    hRadPro->Fill(timeCalib->time.at(i),timeCalib->nhit5by5.at(i));
  }
  //TCanvas *cc=new TCanvas();
  //hRadPro->Draw("colz");
  //cc->Update();
  //cc->WaitPrimitive();
  hRadPro->ProfileX()->Fit(func,"NQ");
  //cc->Update();
  //cc->WaitPrimitive();
  out << func->GetParameter(0) << " " << func->GetParError(0) 
      << " " << func->GetParameter(1) << " " << func->GetParError(1) 
      << " " << func->GetChisquare()/func->GetNDF() << std::endl;
  delete hRadPro;
  out.close();
  delete func;
  delete timeCalib;
  //delete cc;
}

int main()
{

  float const conS =.3/2.35;
  ofstream fout("cal/backN.cal");
  ofstream fwhm("cal/fwhmback.dat");
  TFile f("sort.root");
  TCanvas* canvas[14];
  int Ntele = 14;
  int Nstrip = 32;
  ostringstream outstring;
  TH2I frame("frame","",10,4.5,9.5,10,0,130);
  frame.SetStats(kFALSE);

  double xx[14*32];
  double yy[14*32];

  TF1 *func = new TF1("fit",ThPeaks,3,9,4);
  double para[5];
  ifstream file("cal/back.cal");
  float intercept, slope;
  int i1,i2;
  string name;

  TH1F con("con","",500,0,10);
  for (int itele=0;itele<Ntele;itele++)
    {
      outstring.str("");
      outstring << "B"<<itele;
      name = outstring.str();
      canvas[itele] = new TCanvas(name.c_str());
      canvas[itele]->Divide(6,6);
      for (int istrip =0;istrip<Nstrip;istrip++)
        {
      
          canvas[itele]->cd(istrip+1);
          file >> i1 >> i2  >> slope >> intercept;


          outstring.str("");
          outstring << "back/cal/EBC"<<itele<<"_"<<istrip;
          string name = outstring.str();
          cout <<  name << endl;
          TH1I * hist = (TH1I*) f.Get(name.c_str());
	  frame.Draw();


          hist->SetStats(kFALSE);
          hist->GetXaxis()->SetRangeUser(4.5,9.5);
          con.GetXaxis()->SetRangeUser(4.5,9.5);
	  for (int i=1;i<=500;i++)
	    for (int j=1;j<500;j++)
	    {
              float deltax = hist->GetBinCenter(i)-con.GetBinCenter(j);
	      if (fabs(deltax) > 10.*conS)continue;
              float fact = gauss(deltax,0.,conS);
	      float y = fact*hist->GetBinContent(i)*hist->GetBinWidth(i);
              con.SetBinContent(j,y+con.GetBinContent(j));
	    }



	  for (int i=1;i<=500;i++) 
	    {
	      hist->SetBinContent(i,con.GetBinContent(i));
	      con.SetBinContent(i,0.);
	    }

          hist->Draw("same");



          func->SetParameter(0,0);
          func->SetParameter(1,1.);
	  func->FixParameter(2,conS);
          //func->SetParameter(2,0.1);
          func->SetParameter(3,8.);
          func->SetLineColor(2);
          //func->Draw("same");



	  hist->Fit(func);
          func->GetParameters(para);
          cout << "chisq=" << func->GetChisquare() << endl;
           if (fabs(para[1]-1.) < .2) 
	     { 

              slope *= para[1];
              intercept = intercept*para[1] + para[0];
	     }
            fout << itele << " " << istrip << " " 
                 << slope << " " << intercept << endl;
            fwhm << itele << " " << istrip << " " 
                  << para[2]*2.35 << endl;
            int ii = itele*32+istrip;
            xx[ii] = (float)ii;
            yy[ii] = para[2]*2.35;
            cout << para[0] << " " << para[1] << " " << para[2] << endl;

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

//.........这里部分代码省略.........
	histXZ->SetMarkerStyle(kOpenTriangleDown); //  - open triangle down not found on noise server ? 
	
	double * xc = new double[3];
	double * yc = new double[3];
	double * zc = new double[3];
	
	vector<double> coordinates ;
	double xCoordinate = 0;
	int yCoordinate = 0;
	int zCoorinate = 0;
	
	coordinates = this->getChamberNumber(1)->getXYCoordinatesOfCluster(i+1);
	xCoordinate = coordinates.at(0);
	yCoordinate = coordinates.at(1);
	zCoorinate = 10*vectorOfReferenceChambers[0];
	
	int prevPartition = yCoordinate;
	xc[0] = xCoordinate;
	yc[0] = yCoordinate;
	zc[0] = 1*10;
	
	histXZ->Fill(zc[0],xCoordinate);
	cout << xCoordinate << " " << yCoordinate << endl;
	
	coordinates = this->getChamberNumber(4)->getXYCoordinatesOfCluster(j+1);
	xCoordinate = coordinates.at(0);
	yCoordinate = coordinates.at(1);
	
	xc[1] = xCoordinate;
	yc[1] = yCoordinate;
	zc[1] = 4*10;
	
	histXZ->Fill(zc[1],xCoordinate);
	prevPartition = yCoordinate;
	
	cout << xCoordinate << " " << yCoordinate << endl;	
	coordinates = this->getChamberNumber(6)->getXYCoordinatesOfCluster(k+1);
	xCoordinate = coordinates.at(0);
	yCoordinate = coordinates.at(1);
	
	xc[2] = xCoordinate;
	yc[2] = yCoordinate;
	zc[2] = 6*10;
	
	histXZ->Fill(zc[2],xCoordinate);
	cout << xCoordinate << " " << yCoordinate << endl;
	
	if ( positive && negative ) continue;
	
	TF1 * fitfunc = new TF1("FitTrack","[0]+x*[1]",0,100);
	Double_t * params = new Double_t[2];
	histXZ->Fit(fitfunc);
	fitfunc->GetParameters(params);
	cout << "par1 " << params[0] << " par2 " << params[1] << " chi2 " << fitfunc->GetChisquare() << endl;
	double channelToSearchHitIn ;
	
	for (int jj = 0 ; jj < this->getNumberOfChambers() ; jj++){
	  if (jj+1 != vectorOfReferenceChambers[0] || jj+1 != vectorOfReferenceChambers[1] || jj+1 != vectorOfReferenceChambers[1])
	    // add additional rule that the chamber should exist in the calibration object 
	  {
	    
	    channelToSearchHitIn = fitfunc->Eval((jj+1)*10);
	    cout << "Evaluated for chamber number " << jj+1 << " value : " << channelToSearchHitIn << endl;
	    
	  }
	  
	}
	// now here - what to return, and how to get the hits in the chambers under test from the function
	
	
	
	if (fitfunc->GetChisquare() > 20) continue; // cut the execution 
	if(fitfunc->GetChisquare() < 20){
	  //goodTracks->Write(trackHistoName.c_str());
	  histClustersPartitionDistr->Fill(partitionPenetrated);
	  //histXZ->SaveAs((trackHistoName+".root").c_str());  
	  // here search for hits in the chambers under test  
	  
	}
	
	/*
	else{
	  badTracks->Write(trackHistoName.c_str());
	}	
	*/
	
	//trackHistoName+=".root";
	//histXZ->SaveAs(trackHistoName.c_str());
	
	histXZ->Delete();
	
      }
    }
  }
//   badTracks->Close("R");
//   badTracks->Delete();
//   goodTracks->Close("R");
//   goodTracks->Delete();

}