C++ TH1D::GetMinimum方法代码示例

TH2D* GetJetCorrFunc2D_doublegaussian(int itrg, int jass)
{
                    TH2D* hcorr = (TH2D*)GetRawCorrFunc2D_ratio(itrg,jass);
                    TH2D* hcorr_clone = (TH2D*)hcorr->Clone(Form("corr_clone_itrg%d_jass%d",itrg,jass));
                    hcorr_clone->Reset();
		    for(int ietabin=1;ietabin<=hcorr->GetNbinsX();ietabin++)
		    {
                      TH1D* hcorrphi = (TH1D*)hcorr->ProjectionY(Form("corrphi_%d",ietabin),ietabin,ietabin,"e");  
                      float min = hcorrphi->GetMinimum();
                      hcorrphi->SetAxisRange(-1,1,"X");
                      float nearmax = hcorrphi->GetMaximum();
                      hcorrphi->SetAxisRange(PI-1,PI+1,"X");
                      float awaymax = hcorrphi->GetMaximum();

		      TF1* fitfunc = new TF1("fitfunc",doubleGaussian,-PI/2.,3.*PI/2.,5);		    
                      fitfunc->SetParameters(min,nearmax-min,0.3,awaymax-min,0.5);  
                      fitfunc->SetParLimits(0,0,100000);
                      fitfunc->SetParLimits(1,0,100000);
                      fitfunc->SetParLimits(2,0,100000);
                      fitfunc->SetParLimits(3,0,100000);
                      fitfunc->SetParLimits(4,0,100000);		    
		      for(int ifit=0;ifit<3;ifit++) hcorrphi->Fit(Form("fitfunc_%d",ietabin),"RNO");                     
		      float level = fitfunc->GetParameter(0);

		      for(int iphibin=1;iphibin<=hcorr->GetNbinsY();iphibin++) 
                        hcorr_clone->SetBinContent(ietabin,iphibin,hcorr->GetBinContent(ietabin,iphibin)-level);                      
                      delete fitfunc;
		    }

                    float max = hcorr_clone->GetBinContent(hcorr_clone->GetMaximumBin());
                    hcorr_clone->SetAxisRange(ymin,max*1.3,"Z");
                    return hcorr_clone;
}

void plotgvsr_MC(){
gStyle->SetOptStat(kFALSE);
TCanvas *c1 = new TCanvas();
c1->SetLogy();
TH1D* hFrame = new TH1D("","",1000,0,1);
hFrame->GetYaxis()->SetRangeUser(5e-9,1);
hFrame->GetXaxis()->SetRangeUser(0.15,0.50);
hFrame->GetXaxis()->SetTitle("r");
hFrame->GetYaxis()->SetTitle("|G^{#theta}(ir)|^{2}");
hFrame->SetTitle("");
hFrame->Draw();
TGraph *gr[nnu];
double r0[nnu];
double G2[nnu];
TLine *l[nnu];
for(int i=0;i<nnu;i++){
gr[i]=plotGF(i,xtheta,r0+i,G2+i,marker[i],color[i]);
gr[i]->Draw("Psame");
}
gcl->Draw("same");
TLegend *tg = new TLegend(0.75,0.70-0.10*nnu,0.90,0.70);
tg->SetFillColor(0);
tg->SetBorderSize(0);
tg->SetTextSize(0.04);
for(int i=0;i<nnu;i++)
if(i==0)
tg->AddEntry(gr[i],Form("Prod"),"lp");
else
tg->AddEntry(gr[i],Form("Sum"),"lp");
tg->AddEntry(gcl,Form("Theory"),"lp");
tg->Draw("same");
TLatex *t= new TLatex();
t->SetNDC();
t->SetTextSize(0.05);
t->SetTextFont(42);
if(xtheta==0)
    t->DrawLatex(0.6,0.2,Form("mult = %d, theta = %d", (trkbin[xbin+1]+trkbin[xbin])/2,xtheta));
    else
    t->DrawLatex(0.6,0.2,Form("mult = %d, theta = #frac{%d}{%d}#pi", (trkbin[xbin+1]+trkbin[xbin])/2,xtheta,ntheta*nn));

for(int i=0;i<nnu;i++){
l[i] = new TLine(r0[i],0,r0[i],G2[i]);
l[i]->SetLineStyle(2);
l[i]->SetLineColor(color[i]);
l[i]->Draw("same");
}
linv->SetLineStyle(2);
linv->SetLineColor(2);
linv->SetY1(hFrame->GetMinimum());
linv->Draw("same");
c1->Print("gvsr_thetas_MC.png");


}

//.........这里部分代码省略.........
	// Create some histograms that monitor the progression of the template
	if (fErrorVsPulseAddedHistograms.find(detname) == fErrorVsPulseAddedHistograms.end()) {
	  std::string error_histname = "hErrorVsPulseAdded_" + detname;
	  std::string error_histtitle = "Plot of the Error as each new Pulse is added to the template for the " + detname + " channel";
	  int n_bins = 10000;
	  TH1D* error_hist = new TH1D(error_histname.c_str(), error_histtitle.c_str(), n_bins,0,n_bins);
	  fErrorVsPulseAddedHistograms[detname] = error_hist;
	}

	// all the other pulses will be fitted to the template and then added to it
	// Get some initial estimates for the fitter
	double template_pedestal = hTemplate->GetBinContent(1);
	double template_amplitude;
	double template_time;

	double pulse_pedestal = hPulseToFit->GetBinContent(1);
	double pulse_amplitude;
	double pulse_time;

	double pedestal_offset_estimate = pulse_pedestal; // now we're dealing with actual pulses since we subtract the template_pedestal in the transformation
	double amplitude_scale_factor_estimate;
	double time_offset_estimate;
	if (TSetupData::Instance()->GetTriggerPolarity(bankname) == 1) { 
	  template_amplitude = (hTemplate->GetMaximum() - template_pedestal);
	  template_time = hTemplate->GetMaximumBin() - 1; // go from bin numbering (1, n_samples) to clock ticks (0, n_samples-1)

	  pulse_amplitude = (hPulseToFit->GetMaximum() - pulse_pedestal);
	  pulse_time = hPulseToFit->GetMaximumBin() - 1;

	  amplitude_scale_factor_estimate = pulse_amplitude / template_amplitude;  // estimated scale factor
	  time_offset_estimate = pulse_time - template_time;
	}
	else if (TSetupData::Instance()->GetTriggerPolarity(bankname) == -1) {
	  template_amplitude = (template_pedestal - hTemplate->GetMinimum());
	  template_time = hTemplate->GetMinimumBin() - 1; // go from bin numbering (1, n_samples) to clock ticks (0, n_samples-1)

	  pulse_amplitude = (pulse_pedestal - hPulseToFit->GetMinimum());
	  pulse_time = hPulseToFit->GetMinimumBin() - 1; // go from bin numbering (1, n_samples) to clock ticks (0, n_samples-1)

	  amplitude_scale_factor_estimate = pulse_amplitude / template_amplitude;  // estimated scale factor
	  time_offset_estimate = pulse_time - template_time;
	}

	fTemplateFitter->SetInitialParameterEstimates(pedestal_offset_estimate, amplitude_scale_factor_estimate, time_offset_estimate);
	
	if (Debug()) {
	  std::cout << "TemplateCreator: " << detname << "(" << bankname << "): Pulse #" << pulseIter - thePulseIslands.begin() << ": " << std::endl
		    << "TemplateCreator: Template: pedestal = " << template_pedestal << ", amplitude = " << template_amplitude << ", time = " << template_time << std::endl
		    << "TemplateCreator: Pulse: pedestal = " << pulse_pedestal << ", amplitude = " << pulse_amplitude << ", time = " << pulse_time << std::endl
		    << "TemplateCreator: Initial Estimates: pedestal = " << pedestal_offset_estimate << ", amplitude = " << amplitude_scale_factor_estimate 
		    << ", time = " << time_offset_estimate << std::endl;
	}

	int fit_status = fTemplateFitter->FitPulseToTemplate(hTemplate, hPulseToFit, bankname);
	++n_fit_attempts;
	if (fit_status != 0) {
	  if (Debug()) {
	    std::cout << "TemplateCreator: Problem with fit (status = " << fit_status << ")" << std::endl;
	  }
	  delete hPulseToFit; // delete this here since it is no longer needed
	  continue;
	}
	++n_successful_fits;

	if (Debug()) {
	  std::cout << "Template Creator: Fitted Parameters: PedOffset = " << fTemplateFitter->GetPedestalOffset() << ", AmpScaleFactor = " << fTemplateFitter->GetAmplitudeScaleFactor()

TH1D* readData(TString dataset, TString name, TString wavel, Double_t peak1, Double_t peak2, Double_t intTime, Double_t shift){

    TString fileName = "./Data/" + dataset + ".txt";

    Double_t histoMin = 99999, histoMax = 0;

    ifstream ifs(fileName); if(!ifs.is_open()){cout << "Error. File " << fileName << " not found. Exiting...\n"; return NULL;}

    Double_t entry, wavelength, photons, stepsize;
    
    while(ifs >> entry >> wavelength >> photons){

        if(entry == 1) stepsize = wavelength;
        else if(entry == 2) stepsize = wavelength - stepsize;
        
        if(wavelength < histoMin) histoMin = wavelength;
        if(wavelength > histoMax) histoMax = wavelength;

    }

    ifs.close();

    TH1D *histo = new TH1D(name + wavel + TString::Format(" %1.2f", shift), name + wavel + TString::Format(" %1.2f", shift), Nint(entry), histoMin + shift - stepsize/2.0, histoMax + shift + stepsize/2.0);


    ifs.open(fileName); if(!ifs.is_open()){cout << "Error. File " << fileName << " not found. Exiting...\n"; return NULL;}
    
    while(ifs >> entry >> wavelength >> photons){
        
        for(Int_t x = 0; x< photons/(1000.0/intTime); x++){
            histo->Fill(wavelength + shift);
        }

    }

    ifs.close();

    TF1 *f;
    if(peak2 == 0.){
        f = new TF1("fit", "[0]+[1]*TMath::Voigt(x-[2],[3],[4])");
        f->SetParameter(0, histo->GetMinimum());
        f->SetParameter(1, histo->GetMaximum()/25.);
        f->SetParameter(2, peak1 + shift);
        f->SetParameter(3, .02);
        f->SetParameter(4, .02);
        f->SetParName(0, "Offset");
        f->SetParName(1, "Normalization");
        f->SetParName(2, "Mean");
        f->SetParName(3, "Sigma");
        f->SetParName(4, "Gamma");
    }
    else{
        f = new TF1("fit", "[0]+[1]*TMath::Voigt(x-[2],[3],[4])+[5]*TMath::Voigt(x-[6],[7],[8])");
        f->SetParameter(0, histo->GetMinimum());
        f->SetParameter(1, histo->GetMaximum()/20.);
        f->SetParameter(2, peak1+shift-0.1);
        f->SetParameter(3, .02);
        f->SetParameter(4, .02);
        f->SetParameter(5, histo->GetMaximum()/20.);
        f->SetParameter(6, peak2+shift+0.1);
        f->SetParameter(7, .02);
        f->SetParameter(8, .02);
        f->SetParName(0, "Offset");
        f->SetParName(1, "Normalization 1");
        f->SetParName(2, "Mean 1");
        f->SetParName(3, "Sigma 1");
        f->SetParName(4, "Gamma 1");
        f->SetParName(5, "Normalization 2");
        f->SetParName(6, "Mean 2");
        f->SetParName(7, "Sigma 2");
        f->SetParName(8, "Gamma 2");
    }

    histo->SetLineWidth(2);
    histo->SetMarkerStyle(20);
    histo->SetMarkerSize(1.5);
    histo->SetMarkerColor(kRed);
    histo->SetLineColor(kBlack);
    histo->GetXaxis()->SetTitle("Wavelength (#AA)");
    histo->GetXaxis()->CenterTitle();
    histo->GetXaxis()->SetTitleSize(0.045);
    histo->GetXaxis()->SetTitleOffset(1.00);
    histo->GetXaxis()->SetLabelOffset(0.010);
    histo->GetXaxis()->SetLabelSize(0.045);
    histo->GetYaxis()->SetTitle("Photon count");
    histo->GetYaxis()->CenterTitle();
    histo->GetYaxis()->SetTitleSize(0.050);
    histo->GetYaxis()->SetTitleOffset(0.95);
    histo->GetYaxis()->SetLabelSize(0.045);
    gStyle->SetTitleSize(0.070, "t");
    if(peak2 == 0.) histo->SetTitle(name + TString::Format(" %4.2f #AA", f->GetParameter(2)));
    else histo->SetTitle(name + TString::Format(" %4.2f / %4.2f #AA", f->GetParameter(2), f->GetParameter(6)));
    //histo->SetTitle("Sodium 5890.522 / 5890.567 #AA");
    cout << shift << endl;

    f->SetLineColor(kBlue);
    f->SetLineWidth(4);
    cout << "Fitting " << name << endl;
    histo->Fit(f, "ME");

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

//.........这里部分代码省略.........
    z_inv->SetLineColor(kBlack);
    z_inv->SetLineStyle(1);
    z_inv->SetLineWidth(1);
    z_jets->SetLineColor(kBlack);
    z_jets->SetLineStyle(3);
    z_jets->SetLineWidth(1);
  }
  lm0->SetLineColor(kRed);
  lm0->SetLineStyle(1);
  lm0->SetLineWidth(2);
  lm1->SetLineColor(kRed);
  lm1->SetLineStyle(3);
  lm1->SetLineWidth(2);

  // Populate legend
  legend->AddEntry( qcd, " QCD", "f" );
  legend->AddEntry( lm0, " SUSY LM0", "L" );
  legend->AddEntry( lm1, " SUSY LM1", "L" );
  legend->AddEntry( tt_jets, " t#bar{t}+jets", "L" );
  legend->AddEntry( w_jets, " W+jets", "L" );
  if ( combine ) {
    legend->AddEntry( z_all, " Z", "L" );
  } else {
    legend->AddEntry( z_jets, " Z+jets", "L" );
    legend->AddEntry( z_inv, " Z#rightarrow#nu#nu", "L" );
  }

  // Calc maximum number of entries
  double aMax = 0.;
  if ( qcd->GetMaximum()     > aMax ) { aMax = qcd->GetMaximum(); }
  if ( lm0->GetMaximum()     > aMax ) { aMax = lm0->GetMaximum(); }
  if ( lm1->GetMaximum()     > aMax ) { aMax = lm1->GetMaximum(); }
  if ( tt_jets->GetMaximum() > aMax ) { aMax = tt_jets->GetMaximum(); }  
  if ( w_jets->GetMaximum()  > aMax ) { aMax = w_jets->GetMaximum(); }  
  if ( combine ) {
    if ( z_all->GetMaximum()  > aMax ) { aMax = z_all->GetMaximum(); }  
  } else {
    if ( z_inv->GetMaximum()   > aMax ) { aMax = z_inv->GetMaximum(); }  
    if ( z_jets->GetMaximum()  > aMax ) { aMax = z_jets->GetMaximum(); }  
  }

  // Calc minimum number of entries
  double aMin = 1.e12;
  if ( qcd->GetMinimum(1.e-12)     < aMin ) { aMin = qcd->GetMinimum(1.e-12); }
  if ( lm0->GetMinimum(1.e-12)     < aMin ) { aMin = lm0->GetMinimum(1.e-12); }
  if ( lm1->GetMinimum(1.e-12)     < aMin ) { aMin = lm1->GetMinimum(1.e-12); }
  if ( tt_jets->GetMinimum(1.e-12) < aMin ) { aMin = tt_jets->GetMinimum(1.e-12); }  
  if ( w_jets->GetMinimum(1.e-12)  < aMin ) { aMin = w_jets->GetMinimum(1.e-12); }  
  if ( combine ) {
    if ( z_all->GetMinimum(1.e-12)   < aMin ) { aMin = z_all->GetMinimum(1.e-12); }  
  } else {
    if ( z_inv->GetMinimum(1.e-12)   < aMin ) { aMin = z_inv->GetMinimum(1.e-12); }  
    if ( z_jets->GetMinimum(1.e-12)  < aMin ) { aMin = z_jets->GetMinimum(1.e-12); }  
  }

  if ( qcd ) qcd->GetYaxis()->SetTitleOffset(1.43);
  if ( qcd ) qcd->GetYaxis()->SetTitleSize(0.06);
  if ( qcd ) qcd->GetXaxis()->SetTitleSize(0.06);
  if ( qcd ) qcd->GetXaxis()->SetTitleOffset(0.9);

  if ( log ) {
    if ( qcd ) qcd->SetMaximum( aMax * 10. );
    if ( qcd ) qcd->SetMinimum( aMin * 0.1 );
  } else {
    if ( qcd ) qcd->SetMaximum( aMax * 1.1 );
    if ( qcd ) qcd->SetMinimum( aMin * 0.9 );
  }

  if ( norm ) {
    if ( qcd ) qcd->DrawNormalized("Ehist");
    if ( lm0->GetEntries() > 0. )     { lm0->DrawNormalized("hsame"); }
    if ( lm1->GetEntries() > 0. )     { lm1->DrawNormalized("hsame"); }
    if ( tt_jets->GetEntries() > 0. ) { tt_jets->DrawNormalized("hsame"); }
    if ( w_jets->GetEntries() > 0. )  { w_jets->DrawNormalized("hsame"); }
    if ( combine ) {
      if ( z_all->GetEntries() > 0. )   { z_all->DrawNormalized("hsame"); }
    } else {
      if ( z_inv->GetEntries() > 0. )   { z_inv->DrawNormalized("hsame"); }
      if ( z_jets->GetEntries() > 0. )  { z_jets->DrawNormalized("hsame"); }
    }
  } else {
    if ( qcd ) qcd->Draw("h");
    lm0->Draw("sameH");
    lm1->Draw("sameH");
    if ( tt_jets ) tt_jets->Draw("sameh");
    w_jets->Draw("sameH");
    if ( combine ) {
      z_all->Draw("sameH");
    } else {
      z_inv->Draw("sameH");
      z_jets->Draw("sameH");
    }
  }
  
  file->cd();
  legend->Draw("same");
  aCanvas->Write();
  return aCanvas;

}

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

//FIXME need ratio plots
    TCanvas *c1 = new TCanvas("CrossSect1D","CrossSect1D",600,600);
    c1->Draw();
    c1->cd();
    TPad *p1 = new TPad("p1", "",0.0,0.25,1.0,0.98,0,0,0);
    p1->Draw();
    p1->cd();
    //p1->SetTickx(kFALSE);

    //TAxis *xaxis = h_theory4->GetXaxis();
    //xaxis->SetMoreLogLabels();
    //xaxis->SetNoExponent();

    TAxis *axis = h_theory4->GetYaxis();
    axis->SetTickLength(axis->GetTickLength()/1.35);
    //axis->SetNdivisions(506);


    p1->SetRightMargin(0.055);
    p1->SetBottomMargin(0.01);
    p1->SetTopMargin(0.1);
    p1->SetLogy();
    p1->SetLogx();
    h_theory4->GetXaxis()->SetMoreLogLabels();
    if (IS_RSHAPE) h_theory4->GetYaxis()->SetTitle("1/#sigma_{Z}d#sigma/dm [GeV^{-1}]");
    else h_theory4->GetYaxis()->SetTitle("d#sigma/dm [pb/GeV]");
    h_theory4->GetXaxis()->SetTitle("m_{#mu#mu} [GeV]");
    h_theory4->SetLineColor(kBlue);
    h_theory4->SetLineWidth(0.1);
    h_theory4->SetMarkerSize(0);
    //h_theory4->SetMaximum(h_theory4->GetMaximum()*100.);
    h_theory4->GetYaxis()->SetRangeUser(h_theory4->GetMinimum()/10.,h_theory4->GetMaximum()*10.);
    h_theory4->SetMinimum(1.5*10E-9);
    h_theory4->Draw("L");

    TH1D* h_theory4_clone = (TH1D*)h_theory4->Clone();
    h_theory4_clone->SetLineWidth(1.5);
    //h_theory4_clone->GetXaxis()->SetRangeUser(40,180);
    h_theory4_clone->Draw("Lsame");
    h_theory4_band->SetLineColor(kBlue);
    h_theory4_band->SetLineWidth(3);
    h_theory4_band->SetMarkerSize(0);
    h_theory4_band->SetFillColor(kBlue);
    h_theory4_band->Draw("E3same");
    //data->Draw("Psame");

    gdata->SetMarkerStyle(20);
    gdata->SetMarkerSize(0.7);
    gdata->Draw("Psame");

    double yminl_ = h_theory4->GetMinimum();
    double ymaxl_ = h_theory4->GetMaximum();
    //draw_bin_grid2( yminl_, 1.7*yminl_ );
    //draw_bin_grid( yminl_, 1.2*yminl_ );
    //draw_bin_grid2( ymaxl_/1.7, ymaxl_ );
    //draw_bin_grid( ymaxl_/1.2, ymaxl_ );

    // text2
    size_t ntxt = 3;
    TString txt[3];
    float txtSize[3];
    float txtX[3];
    float txtY[3];
    int   txtAlign[3];

//.........这里部分代码省略.........
	      break;
	    systTmpMmuo.push_back (dat);
	  }
	inM2.close ();
	// ------------------------------------------------------------

        TH1D *leadingRatioSystematics;
	leadingRatioSystematics = (TH1D *) leading->Clone ("leading");

	if (systTmpM.size () != leadingRatioSystematics->GetNbinsX ())
	  cout << "TE SON MONA! WRONG NUMBER OF BINS (# syst from file->" <<systTmpM.size()<<" - # bins->"<<leadingRatioSystematics->GetNbinsX()<<")"<<endl;
	for (int i = 0; i < leadingRatioSystematics->GetNbinsX (); i++)
	  {
	    //leadingRatioSystematics->SetBinContent(i + 1, 1.0);
	    //leadingRatioSystematics->SetBinError (i + 1,sqrt(systTmpM[i]*systTmpM[i] + systTmpMmuo[i]*systTmpMmuo[i]));
	    //	    leadingRatioSystematics->SetBinContent(i + 1, 1.0);
	    //	    leadingRatioSystematics->SetBinError (i + 1,max(systTmpM[i]/100.0,systTmpMmuo[i]/100.0));
	  }

	plots->cd ();
	TPad *pad1 = new TPad("pad1","pad1",0.01,0.3,0.99,0.99);

	pad1->Draw();
	pad1->cd();
	pad1->SetTopMargin(0.1);
	pad1->SetBottomMargin(0.0);
	pad1->SetLeftMargin(0.2);
	pad1->SetRightMargin(0.0);
	pad1->SetFillStyle(0);
	if (use_case !=3) pad1->SetLogy(1); 
	else pad1->SetLogy(0);

	if (use_case ==3){
	  leading->SetMinimum((0.5-0.05*(whichjet-1))*leading->GetMinimum());
	  leading->SetMaximum((1.25+0.35*(whichjet-1))*leading->GetMaximum());
	}
	leading->SetLineColor (kRed+1);
	leading->SetMarkerStyle (20);
	leading->SetFillColor (kRed+1);
	leading->SetMarkerColor (kRed+1);

	leading->GetXaxis ()->SetTitleOffset (1.1);
	leading->GetXaxis ()->SetTitleSize (0.05);
	leading->GetXaxis ()->SetLabelSize (0.0);
	leading->GetXaxis ()->SetLabelFont (42);
	leading->GetXaxis ()->SetTitleFont (42);

	leading->GetYaxis ()->SetTitleOffset (1.);
	leading->GetYaxis ()->SetTitleSize (0.07);
	leading->GetYaxis ()->SetLabelSize (0.06);
	leading->GetYaxis ()->SetLabelFont (42);
	leading->GetYaxis ()->SetTitleFont (42);

	leading->SetTitle ();
	leading->GetXaxis ()->SetTitle ();
	    

	if (use_case ==1) {
	  if (absoluteNormalization) leading->GetYaxis ()->SetTitle ("d#sigma/dN [pb]");
	  else {
	    if (lepton == 1) leading->GetYaxis ()->SetTitle ("(1/#sigma_{Z #rightarrow e^{+}e^{-}}) d#sigma/dN");
	    if (lepton == 2) leading->GetYaxis ()->SetTitle ("(1/#sigma_{Z #rightarrow #mu^{+}#mu^{-}}) d#sigma/dN");
	    if (lepton == 3) leading->GetYaxis ()->SetTitle ("(1/#sigma_{Z #rightarrow l^{+}l^{-}}) d#sigma/dN");
	  }
	}

//.........这里部分代码省略.........
        c1[iPID]->Update() ;
      else
        c1b[iPID]->Update() ;

      if(i<17)
        c2[iPID]->cd(i) ;
      else
        c2b[iPID]->cd(i-16) ;
      
      fbg1->SetParameters(fit1->GetParameter(3),fit1->GetParameter(4),fit1->GetParameter(5),fit1->GetParameter(6)); 
      fbg2->SetParameters(fit2->GetParameter(5),fit2->GetParameter(6),fit2->GetParameter(7),fit2->GetParameter(8)); 
      
      Double_t intRangeMin = PeakPosition(pt)-3.*PeakWidth(pt) ;
      Double_t intRangeMax = PeakPosition(pt)+3.*PeakWidth(pt) ;
      Int_t    intBinMin   = hp->GetXaxis()->FindBin(intRangeMin) ;
      Int_t    intBinMax   = hp->GetXaxis()->FindBin(intRangeMax) ;
      Double_t errStat     = hpm->Integral(intBinMin,intBinMax); 

      hpm ->Multiply(fbg1) ;
      hpm2->Multiply(fbg2) ;
      if(pt<10.){
        hp  ->Add(hpm ,-1.) ;
        hp2 ->Add(hpm2,-1.) ;
      }
      
      Int_t binPi0 = hp->FindBin(kMean);
      fgs->SetParameters(hp->Integral(binPi0-1,binPi0+1)/3.,fit1->GetParameter(1),0.006) ;
      fgs->SetParLimits(0,0.000,10.*hp->GetMaximum()) ;
      fgs->SetParLimits(1,0.120,0.175) ;
      fgs->SetParLimits(2,0.0055,0.010) ;
      hp->Fit(fgs,"QL","",rangeMin,rangeMax) ;   
      hp->Fit(fgs,"QML","",rangeMin,rangeMax) ;   
      hp->SetMaximum(hp->GetMaximum()*1.4) ;
      hp->SetMinimum(hp->GetMinimum()*1.1) ;
      hp->SetMarkerStyle(20) ;
      hp->SetMarkerSize(0.7) ;
      mr1[iPID]->SetBinContent(i,fgs->GetParameter(1)) ;
      mr1[iPID]->SetBinError  (i,fgs->GetParError(1) ) ;
      sr1[iPID]->SetBinContent(i,TMath::Abs(fgs->GetParameter(2))) ;
      sr1[iPID]->SetBinError  (i,fgs->GetParError(2) ) ;

      Double_t y=fgs->GetParameter(0)/hp->GetXaxis()->GetBinWidth(1) ;
      nr1[iPID]->SetBinContent(i,y) ;
      Double_t ey=fgs->GetParError(0)/hp->GetXaxis()->GetBinWidth(1) ;
      nr1[iPID]->SetBinError(i,ey) ;

      Double_t npiInt = hp->Integral(intBinMin,intBinMax)-(intBinMax-intBinMin)*fgs->GetParameter(3) ;
      Double_t norm   = fbg1->GetParameter(0) ;
      Double_t normErr= fbg1->GetParError(0) ;
      if(npiInt>0.){
        nr1int[iPID]->SetBinContent(i,npiInt) ;
        nr1int[iPID]->SetBinError(i,TMath::Sqrt(npiInt + norm*errStat + normErr*normErr*errStat*errStat + norm*norm*errStat)) ;
      }
      
//printf(" Nint1 =%f+-%f \n",npiInt,TMath::Sqrt(npiInt + norm*errStat + normErr*normErr*errStat*errStat + norm*norm*errStat)) ;
      hp2->GetXaxis()->SetRangeUser(rangeMin,rangeMax) ;
      hp2->SetMaximum(hp2->GetMaximum()*1.4) ;
      hp2->SetMinimum(hp2->GetMinimum()*1.1) ;
      hp2->SetMarkerStyle(24) ;
      hp2->SetMarkerSize(0.8) ;

      fcb->SetParameters(hp->Integral(binPi0-1,binPi0+1)/3.,fit2->GetParameter(1),0.006,fit2->GetParameter(3),fit2->GetParameter(4)) ;
      fcb->SetParLimits(0,0.000,10.*hp->GetMaximum()) ;
      fcb->SetParLimits(1,0.120,0.175) ;
      fcb->SetParLimits(2,0.005,0.010) ;
      hp2->Fit(fcb,"QL","",rangeMin,rangeMax) ;

//.........这里部分代码省略.........
      if(type == 2) corr = new TH1D (*(TH1D*)file->Get("JF_DEC"));
      if(type == 3) corr = new TH1D (*(TH1D*)file->Get("JF_DIR"));
    }
    else{
       my_corr = new TH1D (*(TH1D*)f11->Get(cfname.str().c_str()));
      if(type == 0) corr = new TH1D (*(TH1D*)file->Get("hDPHI_INC"));
      if(type == 1) corr = new TH1D (*(TH1D*)file->Get("hDPHI_PI0"));
      if(type == 2) corr = new TH1D (*(TH1D*)file->Get("hDPHI_DEC"));
    }
    
    int nReb;
    if(photon_bin == 5) nReb = 3;
    else if(photon_bin == 7) nReb = 5;
    else nReb = 6;
    corr->Rebin(nReb);
    my_corr->Rebin(nReb);
    
    // if(!isjf){
      double my_R = my_corr->Integral("width")/PI;
      cout<<"my_R = "<<my_R<<endl;
      my_corr->Scale(1./my_R);
      double R = corr->Integral("width")/PI;
      cout<<"R = "<<R<<endl;
      //corr->Sumw2();
      cout<<"my_R/R = "<<my_R/R<<endl;
      corr->Scale(1./R);
      //   }
    
    corr->SetMarkerStyle(20);
    corr->SetMarkerSize(0.8);
    my_corr->SetMarkerColor(1);
    my_corr->SetLineColor(1);
    if(isjf){
      double ymin = corr->GetMinimum();
      double ymax = corr->GetMaximum();
      if(my_corr->GetMaximum() > ymax) ymax = my_corr->GetMaximum();
      if(my_corr->GetMinimum() < ymin) ymin = my_corr->GetMinimum();
      corr->SetAxisRange(ymin-0.01, 1.2*ymax,"Y");
    }
    corr->GetYaxis()->SetTitleOffset(1.9);
    corr->Draw();
    my_corr->Draw("same");

    TLegend *l1 = new TLegend(0.7,0.7,0.9,0.9,"","brNDC");
    l1->AddEntry(my_corr,"Run11","lpf");
    l1->AddEntry(corr,"Run10","lpf");
    if(hadron_bin==0) l1->Draw("same");

    TLatex *la = new TLatex(0.3, 0.75, latex_name.str().c_str());
    la->SetNDC();
    la->Draw("same"); 
      
    string pbin_string = cfname.str().substr(7,1);
    string hbin_string = cfname.str().substr(10,1);
    if(pbin_string.compare("0") == 0){
      if(hbin_string.compare("0") == 0) can_ratio[0]->cd(1);
      if(hbin_string.compare("1") == 0) can_ratio[0]->cd(2);
      if(hbin_string.compare("2") == 0) can_ratio[0]->cd(3);
      if(hbin_string.compare("3") == 0) can_ratio[0]->cd(4);
      if(hbin_string.compare("4") == 0) can_ratio[0]->cd(5);
    }
    if(pbin_string.compare("1") == 0){
      if(hbin_string.compare("0") == 0) can_ratio[1]->cd(1);
      if(hbin_string.compare("1") == 0) can_ratio[1]->cd(2);
      if(hbin_string.compare("2") == 0) can_ratio[1]->cd(3);
      if(hbin_string.compare("3") == 0) can_ratio[1]->cd(4);