C++ TH1I::Fit方法代码示例

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;

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

//.........这里部分代码省略.........
         double delta = selectedTAGHOffset->GetBinContent(i);
         if (useRF) delta *= RF_Period;
         if (ccdb_sum > 0.0 && fabs(delta - mpDelta[tdc_slot-1]) > limit) {
            delta = mpDelta[tdc_slot-1];
         }
         double ccdb = tagh_fadc_time_offsets[i-1];
         double offset = ccdb + delta;
         if (i == 1) c1_adcOffset = offset;
         offset -= c1_adcOffset;
         outFile << i << " " << offset << endl;
         if (verbose) printf("ADC\t%i\t%.3f\t\t%.3f\t\t%.3f\t\t%.3f\n", i, delta, ccdb, mpDelta[tdc_slot-1], offset);
      }
      outFile.close();

      outFile.open(prefix + "tagh_base_time.txt");
      outFile << tagh_t_base_fadc - c1_adcOffset << " " << tagh_t_base_tdc - c1_tdcOffset << endl;
      if (verbose) {
         printf("TAGH ADC Base = %f - (%f) = %f\n", tagh_t_base_fadc, c1_adcOffset, tagh_t_base_fadc - c1_adcOffset);
         printf("TAGH TDC Base = %f - (%f) = %f\n", tagh_t_base_tdc, c1_tdcOffset, tagh_t_base_tdc - c1_tdcOffset);
      }
      outFile.close();
   }

   // We can use the RF time to calibrate the SC time (Experimental for now)
   double meanSCOffset = 0.0; // In case we change the time of the SC, we need this in this scope
   if(useRF){
      TH1F * selectedSCSectorOffset = new TH1F("selectedSCSectorOffset", "Selected TDC-RF offset;Sector; Time", 30, 0.5, 30.5);
      TH1F * selectedSCSectorOffsetDistribution = new TH1F("selectedSCSectorOffsetDistribution", "Selected TDC-RF offset;Time;Entries", 100, -3.0, 3.0);
      TF1* f = new TF1("f","pol0(0)+gaus(1)", -3.0, 3.0);
      for (int sector = 1; sector <= 30; sector++){
         TH1I *scRFHist = ExtractTrackBasedTimingNS::Get1DHistogram("HLDetectorTiming", "SC_Target_RF_Compare", Form("Sector %.2i", sector));
         if (scRFHist == NULL) continue;
         //Do the fit
         TFitResultPtr fr = scRFHist->Fit("pol0", "SQ", "", -2, 2);
         double p0 = fr->Parameter(0);

         f->FixParameter(0,p0);
         f->SetParLimits(2, -2, 2);
         f->SetParLimits(3, 0, 2);
         f->SetParameter(1, 10);
         f->SetParameter(2, scRFHist->GetBinCenter(scRFHist->GetMaximumBin()));
         f->SetParameter(3, 0);

         fr = scRFHist->Fit(f, "SQ", "", -2, 2);
         double SCOffset = fr->Parameter(2);
         selectedSCSectorOffset->SetBinContent(sector, SCOffset);
         selectedSCSectorOffsetDistribution->Fill(SCOffset);
      }
      // Now write out the offsets
      meanSCOffset = selectedSCSectorOffsetDistribution->GetMean();
      if (verbose){
         cout << "Dumping SC results...\n=======================================" << endl;
         cout << "SC mean Offset = " << meanSCOffset << endl;
         cout << "TDC Offsets" << endl;
         cout << "Sector\toldValue\tValueToUse\tmeanOffset\tTotal" << endl;
      }
      outFile.open(prefix + "sc_tdc_timing_offsets.txt");
      for (int sector = 1; sector <= 30; sector++){
         outFile << sc_tdc_time_offsets[sector-1] + selectedSCSectorOffset->GetBinContent(sector) - meanSCOffset << endl;
         if (verbose) printf("%i\t%.3f\t\t%.3f\t\t%.3f\t\t%.3f\n",sector, sc_tdc_time_offsets[sector-1], selectedSCSectorOffset->GetBinContent(sector), meanSCOffset,
               sc_tdc_time_offsets[sector-1] + selectedSCSectorOffset->GetBinContent(sector) - meanSCOffset);
      }
      outFile.close();
      if (verbose){
         cout << "ADC Offsets" << endl;
         cout << "Sector\tvalueToUse\toldValue\tmeanOffset\tTotal" << endl;

void CsIProj()
{
  TFile *file = new TFile("../root/NZ_55_New.root");
  TFile *gates = new TFile("../gates/zlines.root");
  TFile *gates2 = new TFile("../gates/zlines_new.root");

  ofstream ofile("CsI_55A_New.dat");
  
  TCanvas *mycan = (TCanvas*)gROOT->FindObjectAny("mycan");
  if(!mycan)
    {
      mycan = new TCanvas("mycan","mycan");
      mycan->Divide(1,2);
    }
  
  
  ostringstream outstring;
  string name;
  int p1= 30, p2=50; //+- fit limits up to 2 peaks. May be different.
  int const num_par = 5; //number of peaks times 2(pol1)+3(gaus).
  
  for(int ic =0;ic<56;ic++)
    {
      
      outstring.str("");
      outstring << "dEE/dEE_" << ic;
      name = outstring.str();
      
      mycan->cd(1);
      TH2I *hist = (TH2I*)file->Get(name.c_str());
      hist->Draw("col");  
      hist->GetXaxis()->SetRangeUser(200.,1800.); 
      hist->GetYaxis()->SetRangeUser(5.,50.); 
   
      
      if(ic <16 || ic > 31)
	TCutG *mycut = (TCutG*)gates->Get(Form("Zline_%i_2_4",ic));
      else 
	TCutG *mycut = (TCutG*)gates2->Get(Form("Zline_%i_2_4",ic));
      mycut->Draw();

      file->cd();
      outstring.str("");
      outstring << "CsI/CsIGate/ECsI_" << ic << "_Gate";
      name = outstring.str();
      gPad->SetLogz();

      mycan->cd(2);
      TH1I * proj = (TH1I*)file->Get(name.c_str());
      proj->Draw();
      proj->Rebin(4);
      proj->GetXaxis()->SetRangeUser(700.,1800.);

      mycan->Modified();
      mycan->Update();

      TMarker * mark;
      mark=(TMarker*)mycan->WaitPrimitive("TMarker"); //Get the Background limits
      int bkg_lo = mark->GetX();
      delete mark;  
      mark=(TMarker*)mycan->WaitPrimitive("TMarker");
      int bkg_hi = mark->GetX();
      delete mark;
      mark=(TMarker*)mycan->WaitPrimitive("TMarker"); // Get the 1st peak initial guess
      int peak1 = mark->GetX();
      delete mark;
      
      
      double par[num_par] = {0.};
      double out[num_par] = {0.}; 
      int peak1_lo = peak1 - p1, peak1_hi = peak1 + p1; // Peak center and limits
      
      
      TF1 *l1 = new TF1("l1", "pol1", bkg_lo, bkg_hi);
      TF1 *g1 = new TF1("g1", "gaus", peak1_lo,peak1_hi);
      
      TF1 *total = new TF1("total", "pol1(0)+gaus(2)", bkg_lo,bkg_hi);
      
      proj->Fit(l1,"R");
      proj->Fit(g1,"R+");
      
      l1->GetParameters(&par[0]);
      g1->GetParameters(&par[2]);
      
      total->SetParameters(par);
      proj->Fit(total,"R");
      total->GetParameters(out);
      
      
      ofile << ic << " " << out[3] << endl;
      
      outstring.str("");
      outstring << "55A_" << ic;
      name = outstring.str();
      total->SetName(name.c_str());
      total->Draw("same");
      mycan->Modified();
      mycan->Update();
      
      bool IsGood = 0;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
            float timeWindow = 3; //ns (Full Width)
            int binWindow = int(timeWindow / nsPerBin);
            double maxEntries = 0;
            double maxMean = 0;
            for (int j = 1 ; j <= projY->GetNbinsX();j++){
                int minBin = j;
                int maxBin = (j + binWindow) <= projY->GetNbinsX() ? (j + binWindow) : projY->GetNbinsX();
                double sum = 0, nEntries = 0;
                for (int bin = minBin; bin <= maxBin; bin++){
                    sum += projY->GetBinContent(bin) * projY->GetBinCenter(bin);
                    nEntries += projY->GetBinContent(bin);
                    if (bin == maxBin){
                        if (nEntries > maxEntries) {
                            maxMean = sum / nEntries;
                            maxEntries = nEntries;
                        }
                    } 
                }
            }
            //In the case there is RF, our job is to pick just the number of the correct beam bunch, so that's really all we need.
            if(useRF) {
                int beamBucket = int((maxMean / RF_Period) + 0.5); // +0.5 to handle rounding correctly
                selectedTAGMOffset->SetBinContent(i, beamBucket);
                TAGMOffsetDistribution->Fill(beamBucket);
            }
            else{
                selectedTAGMOffset->SetBinContent(i, maxMean);
                TAGMOffsetDistribution->Fill(maxMean);
            }
        }
        /*
        if (!useRF){
            //TFitResultPtr fr1 = selectedTAGMOffset->Fit("pol1", "SQ", "", 0.5, nBinsX + 0.5);
            TFitResultPtr fr1 = selectedTAGMOffset->Fit("pol1", "SQ", "", 5, 50);

            for (int i = 1 ; i <= nBinsX; i++){
                double x0 = fr1->Parameter(0);
                double x1 = fr1->Parameter(1);
                //double x2 = fr1->Parameter(2);
                //double fitResult = x0 + i*x1 + i*i*x2;
                double fitResult = x0 + i*x1;

                double outlierCut = 20;
                double valueToUse = selectedTAGMOffset->GetBinContent(i);
                if (fabs(selectedTAGMOffset->GetBinContent(i) - fitResult) > outlierCut && valueToUse != 0.0){
                    valueToUse = fitResult;
                }

                selectedTAGMOffset->SetBinContent(i, valueToUse);
                if (valueToUse != 0 ) TAGMOffsetDistribution->Fill(valueToUse);
            }
        }
*/
        double meanOffset = TAGMOffsetDistribution->GetMean();
        // This might be in units of beam bunches, so we need to convert
        if (useRF) meanOffset *= RF_Period;
        /*
           for (int i = 1 ; i <= nBinsX; i++){
           double valueToUse = selectedTAGMOffset->GetBinContent(i);
           if (useRF) valueToUse *= RF_Period;
           if (valueToUse == 0) valueToUse = meanOffset;
           outFile.open(prefix + "tagm_tdc_timing_offsets.txt", ios::out | ios::app);
           outFile << "0 " << i << " " << valueToUse + tagm_tdc_time_offsets[i-1] - meanOffset<< endl;
           if (i == 7 || i == 25 || i == 79 || i == 97){
           for(int j = 1; j <= 5; j++){
           outFile << j << " " << i << " " << valueToUse + tagm_tdc_time_offsets[i-1] - meanOffset<< endl;