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

double extractLimitAtQuantile(TString inFileName, TString plotName, double d_quantile ){
	TFile *f = TFile::Open(inFileName);
	TF1 *expoFit = new TF1("expoFit","[0]*exp([1]*(x-[2]))", rMin, rMax);
	TGraphErrors *limitPlot_ =  new TGraphErrors();

/* 	bool done = false; */
	if (_debug > 0) std::cout << "Search for upper limit using pre-computed grid of p-values" << std::endl;

	readAllToysFromFile(limitPlot_, f, d_quantile ); 
	f->Close();
	limitPlot_->Sort();
	double minDist=1e3;
	int n= limitPlot_->GetN();
	cout<<" Number of points in limitPlot_ : "<<n<<endl;
	if(n<=0) return 0;

	clsMin.first=0;
	clsMin.second=0;
	clsMax.first=0;
	clsMax.second=0;

	limit = 0; limitErr = 0;
	for (int i = 0; i < n; ++i) {
		double x = limitPlot_->GetX()[i], y = limitPlot_->GetY()[i]; //, ey = limitPlot_->GetErrorY(i);
		if (fabs(y-clsTarget) < minDist) { limit = x; minDist = fabs(y-clsTarget); }
	}
	int ntmp =0;
	for (int j = 0; j < n; ++j) {
		int i = n-j-1;
		double x = limitPlot_->GetX()[i], y = limitPlot_->GetY()[i], ey = limitPlot_->GetErrorY(i);
		if (y-3*ey >= clsTarget && ntmp<=2) { 
			rMin = x; clsMin = CLs_t(y,ey); 
			ntmp ++ ;
		}
	}
	ntmp =0;
	for (int i = 0; i < n; ++i) {
		double x = limitPlot_->GetX()[i], y = limitPlot_->GetY()[i], ey = limitPlot_->GetErrorY(i);
		if (y+3*ey <= clsTarget && ntmp<=2) {
		       	rMax = x; clsMax = CLs_t(y,ey);
			ntmp ++ ;
		}
	}
	if((clsMin.first==0 and clsMin.second==0) )
	{
		rMin = limitPlot_->GetX()[0]; clsMin=CLs_t(limitPlot_->GetY()[0], limitPlot_->GetErrorY(0)); 
	}
	if((clsMax.first==0 and clsMax.second==0))
	{
		rMax = limitPlot_->GetX()[n-1]; clsMax=CLs_t(limitPlot_->GetY()[n-1], limitPlot_->GetErrorY(n-1));
	}

	if (_debug > 0) std::cout << " after scan x ~ " << limit << ", bounds [ " << rMin << ", " << rMax << "]" << std::endl;
	limitErr = std::max(limit-rMin, rMax-limit);
	expoFit->SetRange(rMin,rMax);
	//expoFit.SetRange(limitPlot_->GetXaxis()->GetXmin(),limitPlot_->GetXaxis()->GetXmax());
	//expoFit->SetRange(1.7,2.25);

	if (limitErr < std::max(rAbsAccuracy_, rRelAccuracy_ * limit)) {
		if (_debug > 1) std::cout << "  reached accuracy " << limitErr << " below " << std::max(rAbsAccuracy_, rRelAccuracy_ * limit) << std::endl;
/* 		done = true;  */
	}

	//if (!done) { // didn't reach accuracy with scan, now do fit
	if (1) { // didn't reach accuracy with scan, now do fit
		if (_debug) {
			std::cout << "\n -- HybridNew, before fit -- \n";
			std::cout << "Limit: r" << " < " << limit << " +/- " << limitErr << " [" << rMin << ", " << rMax << "]\n";

			std::cout<<"rMin="<<rMin<<" clsMin="<<clsMin.first<<",  rMax="<<rMax<<" clsMax="<<clsMax.first<<endl;
		}

		expoFit->FixParameter(0,clsTarget);
		expoFit->SetParameter(1,log(clsMax.first/clsMin.first)/(rMax-rMin));
		expoFit->SetParameter(2,limit);
		double rMinBound, rMaxBound; expoFit->GetRange(rMinBound, rMaxBound);
		limitErr = std::max(fabs(rMinBound-limit), fabs(rMaxBound-limit));
		int npoints = 0; 
		for (int j = 0; j < limitPlot_->GetN(); ++j) { 
			if (limitPlot_->GetX()[j] >= rMinBound && limitPlot_->GetX()[j] <= rMaxBound) npoints++; 
		}
		for (int i = 0, imax = 0; i <= imax; ++i, ++npoints) {
			limitPlot_->Sort();
			limitPlot_->Fit(expoFit,(_debug <= 1 ? "QNR EX0" : "NR EXO"));
			if (_debug) {
				std::cout << "Fit to " << npoints << " points: " << expoFit->GetParameter(2) << " +/- " << expoFit->GetParError(2) << std::endl;
			}

			// only when both  "cls+3e<0.05 and cls-3e>0.05" are satisfied, we require below ...
			//	if ((rMin < expoFit->GetParameter(2))  && (expoFit->GetParameter(2) < rMax) && (expoFit->GetParError(2) < 0.5*(rMaxBound-rMinBound))) { 
			// sanity check fit result
			limit = expoFit->GetParameter(2);
			limitErr = expoFit->GetParError(2);
			if (limitErr < std::max(rAbsAccuracy_, rRelAccuracy_ * limit)) break;
			//	}
		} 
	}

	if (limitPlot_) {
		TCanvas *c1 = new TCanvas("c1","c1");
//.........这里部分代码省略.........

//.........这里部分代码省略.........
         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;
      }
      outFile.open(prefix + "sc_adc_timing_offsets.txt");
      for (int sector = 1; sector <= 30; sector++){
         outFile << sc_fadc_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_fadc_time_offsets[sector-1], selectedSCSectorOffset->GetBinContent(sector), meanSCOffset,
               sc_fadc_time_offsets[sector-1] + selectedSCSectorOffset->GetBinContent(sector) - meanSCOffset);

//.........这里部分代码省略.........
  TH1D *diMuonsPt_GenA1[100];
  TH1D *diMuonsPt_RecA1[100];
  TF1 *backfun_1;
  char namePt_1B[500];//for bkg func
 
  for(Int_t ih = 0; ih < Nptbin; ih++){
    diMuonsInvMass_RecA1[ih] = diMuonsInvMass_RecA[0][ih];
    diMuonsInvMass_GenA1[ih] = diMuonsInvMass_GenA[0][ih];
    diMuonsPt_GenA1[ih] = diMuonsPt_GenA[0][ih];
    diMuonsPt_RecA1[ih] = diMuonsPt_RecA[0][ih];
    
    for (int ifile = 1; ifile <= 5; ifile++) {
      diMuonsInvMass_RecA1[ih]->Add(diMuonsInvMass_RecA[ifile][ih]);
      diMuonsInvMass_GenA1[ih]->Add(diMuonsInvMass_GenA[ifile][ih]);     
      
      diMuonsPt_GenA1[ih]->Add(diMuonsPt_GenA[ifile][ih]); 
      diMuonsPt_RecA1[ih]->Add(diMuonsPt_RecA[ifile][ih]); 
    }
  }
  
  //===========================Fitting===================================================================//
  // Fit ranges
  double mass_low, mass_high;
  double MassUpsilon, WeidthUpsilon;
  
  // Low mass range upsilon width 54 KeV
  MassUpsilon = 9.46; WeidthUpsilon = 0.055;
  //MassUpsilon = 9.46; WeidthUpsilon = 0.068;
  mass_low = 9.0; mass_high = 10.0;  // Fit ranges
  
  // Fit Function crystall ball
  TF1 *GAUSPOL = new TF1("GAUSPOL",CrystalBall,8.0,12.0,6);
  GAUSPOL->SetParNames("Yield (#Upsilon)","BinWidth","Mean","Sigma","#alpha","n");
  GAUSPOL->SetParameter(2, MassUpsilon);
  GAUSPOL->SetParameter(3, WeidthUpsilon);
  //GAUSPOL->SetParLimits(3, 0.1*WeidthUpsilon,2.0*WeidthUpsilon);
  GAUSPOL->SetParameter(4, 1.0);
  GAUSPOL->SetParameter(5, 20.0);
  GAUSPOL->SetLineWidth(2.0);
  GAUSPOL->SetLineColor(2);




  //=====================Loop for eff===========================================================
  double GenNo[100]={0};
  double Eff[100]={0};
  double GenError[100]={0};
  double RecError[100]={0};
  double errEff_cat_S1[100]={0};
  double errEff_cat_S2[100]={0};
  double errEff_cat_S1_1[100]={0},errEff_cat_S1_2[100]={0};
  double errEff_cat_S2_1[100]={0},errEff_cat_S2_2[100]={0};
  char PlotName[500],PlotName1[500], PlotName2[500]; 
  char GPlotName[500],GPlotName1[500], GPlotName2[500];

  for (Int_t ih = 0; ih < Nptbin; ih++) {
   
    cout<<" no of gen dimuons from diMuons Pt histo    "<<diMuonsPt_GenA1[ih]->Integral(1,100)<<endl;
    cout<<" no of gen dimuons from diMuons Mass histo  "<<diMuonsInvMass_GenA1[ih]->Integral(1,100)<<endl;
    

    //from pT histogram
    //gen_pt[ih] =diMuonsPt_GenA1[ih]->IntegralAndError(1,100,genError);
   
    gen_pt[ih] = diMuonsInvMass_GenA1[ih]->IntegralAndError(1,100,genError);

//.........这里部分代码省略.........
					{
						cout<<"\nWARNING !!! Impossible to identify the process of the signal from the sample name: "<<samplename<<" -> assuming full chain llqq with l=e or mu ."<<endl;
						BR_sf=1.0;
					}

					h1b->SetBinContent(binToFill,eff/BR_sf);
					sig_masses.push_back(indM*step+startM);
					sig_effs.push_back(eff/BR_sf);
					indM++;

					cout<<"Array indexes: "<<isig<<"  "<<iCat<<endl;
					SigEff[isig][iCat]=eff/BR_sf;
					SigNgen[isig][iCat]=Ngen;//last entry will be used, should be ok because always the same 
					isig++;
					cout<<"isig="<<isig<<"    "<<"nsig="<<nsig<<endl;
				
					if(isig==nsig)iCat++;
				}//end if sample name is a signal sample name

				h1->Fill(samplename,eff);	
			}//end of sample loop

			c1->SetGridy(1);

			/// Here you choose the functional form of the efficiency.
			/// You have to be careful with the initial values of the parameters...
			bool fitEff=false;
			TF1 *fitFunc = 0;
			if(fitEff)
			{
				fitFunc = new TF1("fitPoly2","tanh([0]*(x-[1]))*([2]+[3]*x)",580,2550);
				fitFunc->SetLineColor(kRed);
				fitFunc->SetLineWidth(2);
				fitFunc->SetParameter(0,-0.002);
				fitFunc->SetParameter(1,300);
				fitFunc->SetParameter(2,-0.3);
				fitFunc->SetParameter(3,3E-5);
				fitFunc->SetLineColor(kRed);
				h1b->Fit(fitFunc,"R");

				cout<<"\n\nEff extrapolated at 2000 GeV: "<<fitFunc->Eval(2000.0)<<endl;
			}

			TString lepStr;
			if(lepCut==1)lepStr="MU";
			if(lepCut==0)lepStr="ELE";
			TString PurStr;
			if(iPur==1)PurStr="HP";
			if(iPur==0)PurStr="LP";
			if(nxjCut==2)PurStr="";

			TString SaveName="_"+lepStr+"_"+PurStr+"_";
			SaveName+=iNJ;
			SaveName+="J";
			//outFile<<SaveName<<endl;

			/// The original plot with efficiencies for both signals and backgrounds
			h1->SetTitle(cut);			
			h1->Draw();
			h1->Draw("TEXT0same");
			c1->SaveAs("SignalEffPlots/"+cut+SaveName+".png");

			/// The plot with signal efficiency and the fit.
			TCanvas *c2=new TCanvas("cFit","cFit",800,800);
			c2->cd();
			h1b->SetMarkerStyle(20);

void kees_gen() {
  gROOT->SetStyle("HALLA");
  TCanvas *cn = new TCanvas("cn");
  cn->Draw();
  cn->UseCurrentStyle();
  TH1F *frm = new TH1F("frm","",100,0.,10.);
  frm->GetXaxis()->SetTitle("Q^{2}  [GeV^{2}]");
  frm->GetYaxis()->SetTitle("G_{E}^{n}");
  frm->SetMinimum(-.02);
  frm->SetMaximum(0.1);
  frm->UseCurrentStyle();
  frm->Draw();
  frm->SetAxisRange(0.,5.,"X");

  TF1 *genf = new TF1("genf",genff,1.,10.,1);
  genf->SetLineColor(2);
  genf->SetLineStyle(2);
  genf->SetParameter(0,1.);
  genf->SetParameter(1,.3);
  genf->SetParameter(0,-0.632);
  // match to Madey point just below 1.5
  //  genf->SetParameter(0,.0411/genf->Eval(1.45));
  
  TMultiGraph* mgrDta = new TMultiGraph("Data","G_{E}^{n}");
  TLegend *legDta = new TLegend(.54,.6,.875,.90,"","brNDC");

  TMultiGraph* wgr = mgrDta;
  TLegend *wlg = legDta;

   // the data
  legDta->SetBorderSize(0); // turn off border
  legDta->SetFillStyle(0);
  
  datafile_t *f = datafiles;
  TGraph* gr=0;
  while ( f && f->filename ) {
    gr=OneGraph(f);
    if (gr) {
      if (f->lnpt) {
	mgrDta->Add(gr,f->lnpt);
	legDta->AddEntry(gr,f->label,f->lnpt);
      }
      else if (gr->GetMarkerStyle()>=20) {
	mgrDta->Add(gr,"p");
	legDta->AddEntry(gr,f->label,"p");
      }	
      else {
	mgrDta->Add(gr,"l");
	legDta->AddEntry(gr,f->label,"l");
      }
    }
    f++;
  }
    

  mgrDta->Draw("p");
  //  legDta->Draw();   don't draw the data legend
  
  TMultiGraph* mgrThry = new TMultiGraph("Theory","G_{E}^{n}");
  TLegend *legThry = new TLegend(.54,.6,.875,.9,"","brNDC");

  wgr = mgrThry;
  wlg = legThry;

  // the theory
  wlg->SetBorderSize(0); // turn off border
  wlg->SetFillStyle(0);
  
  f = theoryfiles1;
  gr=0;
  while ( f && f->filename ) {
    gr=OneGraph(f);
    if (gr) {
      TGraphAsymmErrors *egr = dynamic_cast<TGraphAsymmErrors*>(gr);
      if (egr && egr->GetN()>1 && egr->GetEYhigh() && egr->GetEYhigh()[1]>0) {
	gr = toerror_band(egr);
	gr->SetFillStyle(3000+f->style);
      }
      if (f->lnpt) {
	wgr->Add(gr,f->lnpt);
	wlg->AddEntry(gr,f->label,f->lnpt);
      }
      else if (gr->GetMarkerStyle()>=20) {
	wgr->Add(gr,"p");
	wlg->AddEntry(gr,f->label,"p");
      }	
      else {
	wgr->Add(gr,"l");
	wlg->AddEntry(gr,f->label,"l");
      }
    }
    f++;
  }

  genf->Draw("same");
  mgrThry->Draw("c");
  legThry->AddEntry(genf,"F_{2}/F_{1} #propto ln^{2}(Q^{2}/#Lambda^{2})/Q^{2}","l");
  legThry->Draw();

  // draw a line at 1
//.........这里部分代码省略.........

int compareL2L3Correction(const std::string& kalibri, const std::string& jetMETL2, const std::string& jetMETL3) {
//int compareL2L3Correction(const std::string& jetMETL2, const std::string& jetMETL3) {
    // -- Read calibration constants ---------------
    std::cout << "Reading calibration constants...\n";
    parmap corrKalibri  = readParameters(kalibri);
    parmap corrJetMETL2 = readParameters(jetMETL2);
    parmap corrJetMETL3 = readParameters(jetMETL3);
    if( (1 + corrKalibri.size() ) != ( corrJetMETL2.size() + corrJetMETL3.size() ) ) {
        std::cerr << "Linker Kompaktierer defekt. Kontaktieren Sie das Personal (68/111)." << std::endl;
        return 1;
    }


    // -- Global variables -------------------------
    std::vector<TF1*> fCorrKalibri(1+corrKalibri.size());
    std::vector<TF1*> fCorrJetMET(1+corrKalibri.size());
    TF1 * f = 0;
    char name[50];
    std::string outFileName = "comparisonL2L3_Summer08_IC5Calo.ps";
    bool printPar = false;


    // -- Create correction plots for JetMET -------
    std::cout << "Creating JetMET plots...\n";
    // L3 correction
    if( printPar ) std::cout << "  JetMETL3: " << std::flush;
    std::vector<double>& par = corrJetMETL3[0];
    f = new TF1("fJetMETL3",corrL3,par.at(2),par.at(3),4);
    for(int i = 0; i < 4; i++) {
        f->SetParameter(i,par.at(4+i));
        if( printPar ) std::cout << par.at(4+i) << " ";
    }
    if( printPar ) std::cout << "\n";
    fCorrJetMET.at(0) = f;

    // L2 correction
    for(unsigned int bin = 0; bin < corrJetMETL2.size(); bin++) {
        if( printPar ) std::cout << "  JetMETL2 " << bin << ": " << std::flush;
        par = corrJetMETL2[bin];
        sprintf(name,"fJetMETL2_%i",bin);
        f = new TF1(name,corrL2,par.at(2),par.at(3),3);
        sprintf(name,"Bin %i: %.3f <  #eta < %.3f",bin,par.at(0),par.at(1));
        f->SetTitle(name);
        for(int i = 0; i < 3; i++) {
            f->SetParameter(i,par.at(4+i));
            if( printPar ) std::cout << par.at(4+i) << " ";
        }
        if( printPar ) std::cout << "\n";
        fCorrJetMET.at(1+bin) = f;
    }

    // Nice plots
    for(unsigned i = 0; i < fCorrJetMET.size(); i++) {
        fCorrJetMET.at(i)->SetLineWidth(2);
        fCorrJetMET.at(i)->SetLineColor(4);
    }


    // -- Create correction plots for Kalibri ------
    std::cout << "Creating Kalibri plots...\n";
    int nTowerPar     = 0;
    int nJetPar       = 3;
    int nTrackPar     = 0;
    int nGlobalJetPar = 4;
    // Global correction
    if( printPar ) std::cout << "  Kalibri: " << std::flush;
    par = corrKalibri[0];
    f = new TF1("fKalibri_global",corrL3,par.at(2),par.at(3),4);
    for(int i = 0; i < nGlobalJetPar; i++) {
        f->SetParameter(i,par.at(4+nTowerPar+nJetPar+nTrackPar+i));
        if( printPar ) std::cout << par.at(4+nTowerPar+nJetPar+nTrackPar+i) << " ";
    }
//   f = new TF1("fKalibri_global",corrL3,4,2000,4);
//   f->SetParameter(0,0.998293);
//   f->SetParameter(1,5.43056);
//   f->SetParameter(2,3.3444);
//   f->SetParameter(3,2.39809);
    if( printPar ) std::cout << "\n";
    fCorrKalibri.at(0) = f;

    // Local correction
    double scale[3] = { 1., 0.1, 0.01 };
    for(unsigned int bin = 0; bin < corrKalibri.size(); bin++) {
        if( printPar ) std::cout << "  Kalibri " << bin << ": " << std::flush;
        par = corrKalibri[bin];
        sprintf(name,"fKalibri_%i",bin);
        f = new TF1(name,corrL2,par.at(2),par.at(3),3);
        sprintf(name,"Bin %i: %.3f <  #eta < %.3f",bin,par.at(0),par.at(1));
        f->SetTitle(name);
        for(int i = 0; i < nJetPar; i++) {
            f->SetParameter(i,scale[i]*par.at(4+nTowerPar+i));
            if( printPar ) std::cout << scale[i]*par.at(4+nTowerPar+i) << " ";
        }
        if( printPar ) std::cout << "\nDone";
        fCorrKalibri.at(1+bin) = f;
    }

    // Nice plots
    for(unsigned i = 0; i < fCorrKalibri.size(); i++) {
        fCorrKalibri.at(i)->SetLineWidth(2);
//.........这里部分代码省略.........

void FitDijetMass_Data() {

  
  TFile *inf  = new TFile("MassResults_ak7calo.root");
  TH1F *hCorMassDen     = (TH1F*) inf->Get("DiJetMass");
  hCorMassDen->SetXTitle("Corrected Dijet Mass (GeV)");
  hCorMassDen->SetYTitle("Events/GeV");
  hCorMassDen->GetYaxis()->SetTitleOffset(1.5);
  hCorMassDen->SetMarkerStyle(20);
  hCorMassDen->GetXaxis()->SetRangeUser(120.,900.);



  gROOT->ProcessLine(".L tdrstyle.C");
  setTDRStyle();
  tdrStyle->SetErrorX(0.5);
  tdrStyle->SetPadRightMargin(0.08);
  tdrStyle->SetLegendBorderSize(0);
  gStyle->SetOptFit(1111);
  tdrStyle->SetOptStat(0); 

  
  TCanvas* c2 = new TCanvas("c2","DijetMass", 500, 500);
  /////// perform 4 parameters fit
  TF1 *func = new TF1("func", "[0]*((1-x/7000.+[3]*(x/7000)^2)^[1])/(x^[2])", 
  100., 1000.);
  func->SetParameter(0, 1.0e+08);
  func->SetParameter(1, -1.23);
  func->SetParameter(2, 4.13);
  func->SetParameter(3, 1.0);

  func->SetLineColor(4);
  func->SetLineWidth(3);

  TVirtualFitter::SetMaxIterations( 10000 );
  TVirtualFitter *fitter;
  TMatrixDSym* cov_matrix;

  int fitStatus = hCorMassDen->Fit("func","LLI","",130.0, 800.0); // QCD fit
 
  TH1F *hFitUncertainty = hCorMassDen->Clone("hFitUncertainty");
  hFitUncertainty->SetLineColor(5);
  hFitUncertainty->SetFillColor(5);
  hFitUncertainty->SetMarkerColor(5);

  if (fitStatus == 0) {
    fitter = TVirtualFitter::GetFitter();
    double* m_elements = fitter->GetCovarianceMatrix();
    cov_matrix = new TMatrixDSym( func->GetNumberFreeParameters(),m_elements);
    cov_matrix->Print();
    double x, y, e;

    for(int i=0;i<hFitUncertainty->GetNbinsX();i++)
      {
	x = hFitUncertainty->GetBinCenter(i+1);
	y = func->Eval(x);
	e = QCDFitUncertainty( func, *cov_matrix, x);
	hFitUncertainty->SetBinContent(i+1,y);
	hFitUncertainty->SetBinError(i+1,e);
      }
  }

  hCorMassDen->Draw("ep");
  gPad->Update();
  TPaveStats *st = (TPaveStats*)hCorMassDen->FindObject("stats");
  st->SetName("stats1");
  st->SetX1NDC(0.3); //new x start position
  st->SetX2NDC(0.6); //new x end position
  st->SetTextColor(4);
  hCorMassDen->GetListOfFunctions()->Add(st);



  /////// perform 2 parameters fit
  TF1 *func2 = new TF1("func2", "[0]*(1-x/7000.)/(x^[1])", 100., 1000.);
  func2->SetParameter(0, 10000.);
  func2->SetParameter(1, 5.0);
  func2->SetLineWidth(3);

  fitStatus = hCorMassDen->Fit("func2","LLI","",130.0, 800.0); // QCD fit

  TH1F *hFitUncertainty2 = hCorMassDen->Clone("hFitUncertainty2");
  hFitUncertainty2->SetLineColor(kGray);
  hFitUncertainty2->SetFillColor(kGray);
  hFitUncertainty2->SetMarkerColor(kGray);

  if (fitStatus == 0) {
    fitter = TVirtualFitter::GetFitter();
    double* m_elements = fitter->GetCovarianceMatrix();
    cov_matrix = new TMatrixDSym( func2->GetNumberFreeParameters(),m_elements);
    cov_matrix->Print();
    double x, y, e;

    for(int i=0;i<hFitUncertainty2->GetNbinsX();i++)
      {
	x = hFitUncertainty2->GetBinCenter(i+1);
	y = func2->Eval(x);
	e = QCDFitUncertainty( func2, *cov_matrix, x);
	hFitUncertainty2->SetBinContent(i+1,y);
	hFitUncertainty2->SetBinError(i+1,e);
//.........这里部分代码省略.........

TF1 *fit(TTree *nt,TTree *ntMC,double ptmin,double ptmax, bool ispPb, int count){   
   //cout<<cut.Data()<<endl;
   //static int count=0;
   //count++;
   TCanvas *c= new TCanvas(Form("c%d",count),"",600,600);
   TH1D *h = new TH1D(Form("h%d",count),"",50,5,6);
   TH1D *hMC = new TH1D(Form("hMC%d",count),"",50,5,6);

   TString iNP="7.26667e+00*Gaus(x,5.10472e+00,2.63158e-02)/(sqrt(2*3.14159)*2.63158e-02)+4.99089e+01*Gaus(x,4.96473e+00,9.56645e-02)/(sqrt(2*3.14159)*9.56645e-02)+3.94417e-01*(3.74282e+01*Gaus(x,5.34796e+00,3.11510e-02)+1.14713e+01*Gaus(x,5.42190e+00,1.00544e-01))";
   TF1 *f = new TF1(Form("f%d",count),"[0]*([7]*Gaus(x,[1],[2])/(sqrt(2*3.14159)*[2])+(1-[7])*Gaus(x,[1],[8])/(sqrt(2*3.14159)*[8]))+[3]+[4]*x+[5]*("+iNP+")");
   nt->Project(Form("h%d",count),"mass",Form("%s&&pt>%f&&pt<%f",seldata_2y.Data(),ptmin,ptmax));   
   ntMC->Project(Form("hMC%d",count),"mass",Form("%s&&pt>%f&&pt<%f",seldata_2y.Data(),ptmin,ptmax));   
   clean0(h);

   TH1D *hraw = new TH1D(Form("hraw%d",count),"",50,5,6);
   clean0(hraw);
   hraw = (TH1D*)h->Clone(Form("hraw%d",count));


   h->Draw();
   f->SetParLimits(4,-1000,0);
   f->SetParLimits(2,0.01,0.05);
   f->SetParLimits(8,0.01,0.05);
   f->SetParLimits(7,0,1);
   f->SetParLimits(5,0,1000);

   f->SetParameter(0,setparam0);
   f->SetParameter(1,setparam1);
   f->SetParameter(2,setparam2);
   f->SetParameter(8,setparam3);
   f->FixParameter(1,fixparam1);
   h->GetEntries();

   hMC->Fit(Form("f%d",count),"q","",5,6);
   hMC->Fit(Form("f%d",count),"q","",5,6);
   f->ReleaseParameter(1);
   hMC->Fit(Form("f%d",count),"L q","",5,6);
   hMC->Fit(Form("f%d",count),"L q","",5,6);
   hMC->Fit(Form("f%d",count),"L q","",5,6);
   hMC->Fit(Form("f%d",count),"L m","",5,6);

   f->FixParameter(1,f->GetParameter(1));
   f->FixParameter(2,f->GetParameter(2));
   f->FixParameter(7,f->GetParameter(7));
   f->FixParameter(8,f->GetParameter(8));
   
   h->Fit(Form("f%d",count),"q","",5,6);
   h->Fit(Form("f%d",count),"q","",5,6);
   f->ReleaseParameter(1);
   h->Fit(Form("f%d",count),"L q","",5,6);
   h->Fit(Form("f%d",count),"L q","",5,6);
   h->Fit(Form("f%d",count),"L q","",5,6);
   h->Fit(Form("f%d",count),"L m","",5,6);
   h->SetMarkerSize(0.8);
   h->SetMarkerStyle(20);
   cout <<h->GetEntries()<<endl;

   // function for background shape plotting. take the fit result from f
   TF1 *background = new TF1(Form("background%d",count),"[0]+[1]*x");
   background->SetParameter(0,f->GetParameter(3));
   background->SetParameter(1,f->GetParameter(4));
   background->SetParameter(2,f->GetParameter(5));
   background->SetParameter(3,f->GetParameter(6));
   background->SetLineColor(4);
   background->SetRange(5,6);
   background->SetLineStyle(2);
   
   // function for signal shape plotting. take the fit result from f
   TF1 *Bkpi = new TF1(Form("fBkpi%d",count),"[0]*("+iNP+")");
   Bkpi->SetParameter(0,f->GetParameter(5));
   Bkpi->SetLineColor(kGreen+1);
   Bkpi->SetFillColor(kGreen+1);
//   Bkpi->SetRange(5.00,5.28);
   Bkpi->SetRange(5.00,6.00);
   Bkpi->SetLineStyle(1);
   Bkpi->SetFillStyle(3004);

   // function for signal shape plotting. take the fit result from f
   TF1 *mass = new TF1(Form("fmass%d",count),"[0]*([3]*Gaus(x,[1],[2])/(sqrt(2*3.14159)*[2])+(1-[3])*Gaus(x,[1],[4])/(sqrt(2*3.14159)*[4]))");
   mass->SetParameters(f->GetParameter(0),f->GetParameter(1),f->GetParameter(2),f->GetParameter(7),f->GetParameter(8));
   mass->SetParError(0,f->GetParError(0));
   mass->SetParError(1,f->GetParError(1));
   mass->SetParError(2,f->GetParError(2));
   mass->SetParError(7,f->GetParError(7));
   mass->SetParError(8,f->GetParError(8));
   mass->SetLineColor(2);
   mass->SetLineStyle(2);

//   cout <<mass->Integral(0,1.2)<<" "<<mass->IntegralError(0,1.2)<<endl;
   h->SetMarkerStyle(24);
   h->SetStats(0);
   h->Draw("e");
   h->SetXTitle("M_{B} (GeV/c^{2})");
   h->SetYTitle("Entries / (20 MeV/c^{2})");
   h->GetXaxis()->CenterTitle();
   h->GetYaxis()->CenterTitle();
   h->SetTitleOffset(1.5,"Y");
   h->SetAxisRange(0,h->GetMaximum()*1.2,"Y");
   Bkpi->Draw("same");
   background->Draw("same");   
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  sigma_vs_ring[2] = new TGraphErrors();
  sigma_vs_ring[2]->SetMarkerStyle(20);
  sigma_vs_ring[2]->SetMarkerSize(1);
  sigma_vs_ring[2]->SetMarkerColor(kBlue+2);
 
  /// Graph for scale vs ring EE+, EE- and folded

  TGraphErrors *scale_vs_ring[3];
  scale_vs_ring[0] = new TGraphErrors();
  scale_vs_ring[0]->SetMarkerStyle(20);
  scale_vs_ring[0]->SetMarkerSize(1);
  scale_vs_ring[0]->SetMarkerColor(kBlue+2);

  scale_vs_ring[1] = new TGraphErrors();
  scale_vs_ring[1]->SetMarkerStyle(20);
  scale_vs_ring[1]->SetMarkerSize(1);
  scale_vs_ring[1]->SetMarkerColor(kBlue+2);

  scale_vs_ring[2] = new TGraphErrors();
  scale_vs_ring[2]->SetMarkerStyle(20);
  scale_vs_ring[2]->SetMarkerSize(1);
  scale_vs_ring[2]->SetMarkerColor(kBlue+2);
    
  
  TF1 *fgaus = new TF1("fgaus","gaus",-10,10);
  int np[3] = {0};

  /// Gaussian fit for EE+ and EE-

  for (int k = 0; k < 2 ; k++){
    for (int iring = 0; iring < 40; iring++){
      if (hspread[k][iring]-> GetEntries() == 0) continue;
      float e     = 0.5*hcmap[k]-> GetYaxis()->GetBinWidth(1);
      fgaus->SetParameter(1,1);
      fgaus->SetParameter(2,hspread[k][iring]->GetRMS());
      fgaus->SetRange(1-5*hspread[k][iring]->GetRMS(),1+5*hspread[k][iring]->GetRMS());
      hspread[k][iring]->Fit("fgaus","QR");
      sigma_vs_ring[k]-> SetPoint(np[k],iring,fgaus->GetParameter(2));
      sigma_vs_ring[k]-> SetPointError(np[k], e ,fgaus->GetParError(2));
      scale_vs_ring[k]-> SetPoint(np[k],iring,fgaus->GetParameter(1));
      scale_vs_ring[k]-> SetPointError(np[k],e,fgaus->GetParError(1));
      np[k]++;    
    }
  }
  
    for (int iring = 0; iring < 40; iring++){
      if (hspreadAll[iring]-> GetEntries() == 0) continue;
      float e     = 0.5*hcmap[0]-> GetYaxis()->GetBinWidth(1);
      fgaus->SetParameter(1,1);
      fgaus->SetParameter(2,hspreadAll[iring]->GetRMS());
      fgaus->SetRange(1-5*hspreadAll[iring]->GetRMS(),1+5*hspreadAll[iring]->GetRMS());
      hspreadAll[iring]->Fit("fgaus","QR");
      sigma_vs_ring[2]-> SetPoint(np[2],iring,fgaus->GetParameter(2));
      sigma_vs_ring[2]-> SetPointError(np[2], e ,fgaus->GetParError(2));
      scale_vs_ring[2]-> SetPoint(np[2],iring,fgaus->GetParameter(1));
      scale_vs_ring[2]-> SetPointError(np[2],e,fgaus->GetParError(1));
      np[2]++;    
    }

  /// Intercalibration constant vs phi
 
 /* TGraphErrors *IC_vs_phi[2];
  IC_vs_phi[0] = new TGraphErrors();
  IC_vs_phi[0]->SetMarkerStyle(20);
  IC_vs_phi[0]->SetMarkerSize(1);
  IC_vs_phi[0]->SetMarkerColor(kBlue+2);

TF1* fitMass(TH1D* hData, TH1D* hMCSignal, TH1D* hMCSwapped)
{
  Double_t setparam0=100.;
  Double_t setparam1=1.865;
  Double_t setparam2=0.03;
  Double_t setparam10=0.005;
  Double_t setparam8=0.1;
  Double_t setparam9=0.1;
  Double_t fixparam1=1.865;
  Double_t minhisto=1.7;
  Double_t maxhisto=2.0;

  TF1* f = new TF1("fMass","[0]*([7]*([9]*Gaus(x,[1],[2]*(1+[11]))/(sqrt(2*3.1415927)*[2]*(1+[11]))+(1-[9])*Gaus(x,[1],[10]*(1+[11]))/(sqrt(2*3.1415927)*[10]*(1+[11])))+(1-[7])*Gaus(x,[1],[8]*(1+[11]))/(sqrt(2*3.1415927)*[8]*(1+[11])))+[3]+[4]*x+[5]*x*x+[6]*x*x*x", 1.7, 2.0);
  f->SetParLimits(4,-1000,1000);
  f->SetParLimits(10,0.005,0.05);
  f->SetParLimits(2,0.01,0.1);
  f->SetParLimits(8,0.02,0.2);
  f->SetParLimits(7,0,1);
  f->SetParLimits(9,0,1);

  f->SetParameter(0,setparam0);
  f->SetParameter(1,setparam1);
  f->SetParameter(2,setparam2);
  f->SetParameter(10,setparam10);
  f->SetParameter(9,setparam9);

  f->FixParameter(8,setparam8);
  f->FixParameter(7,1);
  f->FixParameter(1,fixparam1);
  f->FixParameter(3,0);
  f->FixParameter(4,0);
  f->FixParameter(5,0);
  f->FixParameter(6,0);
  f->FixParameter(11,0);

  hMCSignal->Fit("fMass","q","",minhisto,maxhisto);
  hMCSignal->Fit("fMass","q","",minhisto,maxhisto);
  f->ReleaseParameter(1);
  hMCSignal->Fit("fMass","L q","",minhisto,maxhisto);
  hMCSignal->Fit("fMass","L q","",minhisto,maxhisto);
  hMCSignal->Fit("fMass","L m","",minhisto,maxhisto);

  f->FixParameter(1,f->GetParameter(1));
  f->FixParameter(2,f->GetParameter(2));
  f->FixParameter(10,f->GetParameter(10));
  f->FixParameter(9,f->GetParameter(9));
  f->FixParameter(7,0);
  f->ReleaseParameter(8);
  f->SetParameter(8,setparam8);

  hMCSwapped->Fit("fMass","L q","",minhisto,maxhisto);
  hMCSwapped->Fit("fMass","L q","",minhisto,maxhisto);
  hMCSwapped->Fit("fMass","L q","",minhisto,maxhisto);
  hMCSwapped->Fit("fMass","L m","",minhisto,maxhisto);

  f->FixParameter(7,hMCSignal->Integral(0,1000)/(hMCSwapped->Integral(0,1000)+hMCSignal->Integral(0,1000)));
  f->FixParameter(8,f->GetParameter(8));
  f->ReleaseParameter(3);
  f->ReleaseParameter(4);
  f->ReleaseParameter(5);
  f->ReleaseParameter(6);

  f->SetLineColor(kRed);

  hData->Fit("fMass","q","",minhisto,maxhisto);
  hData->Fit("fMass","q","",minhisto,maxhisto);
  f->ReleaseParameter(1);
  f->SetParLimits(1,1.86,1.87);
  f->ReleaseParameter(11);
  f->SetParLimits(11,-0.2,0.2);
  hData->Fit("fMass","L q","",minhisto,maxhisto);
  hData->Fit("fMass","L q","",minhisto,maxhisto);
  hData->Fit("fMass","L q","",minhisto,maxhisto);
  hData->Fit("fMass","L m","",minhisto,maxhisto);

  TF1* background = new TF1("fBackground","[0]+[1]*x+[2]*x*x+[3]*x*x*x");
  background->SetParameter(0,f->GetParameter(3));
  background->SetParameter(1,f->GetParameter(4));
  background->SetParameter(2,f->GetParameter(5));
  background->SetParameter(3,f->GetParameter(6));
  background->SetLineColor(4);
  background->SetRange(minhisto,maxhisto);
  background->SetLineStyle(2);

  TF1* mass = new TF1("fSignal","[0]*([3]*([4]*Gaus(x,[1],[2]*(1+[6]))/(sqrt(2*3.1415927)*[2]*(1+[6]))+(1-[4])*Gaus(x,[1],[5]*(1+[6]))/(sqrt(2*3.1415927)*[5]*(1+[6]))))");
  mass->SetParameters(f->GetParameter(0),f->GetParameter(1),f->GetParameter(2),f->GetParameter(7),f->GetParameter(9),f->GetParameter(10),f->GetParameter(11));
  mass->SetParError(0,f->GetParError(0));
  mass->SetParError(1,f->GetParError(1));
  mass->SetParError(2,f->GetParError(2));
  mass->SetParError(3,f->GetParError(7));
  mass->SetParError(4,f->GetParError(9));
  mass->SetParError(5,f->GetParError(10));
  mass->SetFillColor(kOrange-3);
  mass->SetFillStyle(3002);
  mass->SetLineColor(kOrange-3);
  mass->SetLineWidth(3);
  mass->SetLineStyle(2);

  TF1* massSwap = new TF1("fBackground","[0]*(1-[2])*Gaus(x,[1],[3]*(1+[4]))/(sqrt(2*3.1415927)*[3]*(1+[4]))");
  massSwap->SetParameters(f->GetParameter(0),f->GetParameter(1),f->GetParameter(7),f->GetParameter(8),f->GetParameter(11));
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	in_tree->GetEntry(0);

	double* u = new double[sample_len];
	double* t = new double[sample_len];
	in_tree->SetBranchAddress("amplitude", u);
	in_tree->SetBranchAddress("time", t);
	const int nEvents = in_tree->GetEntries();

 	printf("[Pulse Shape] - Found tree with %i events, %i samples per event.\n", nEvents, sample_len);

	TFile* out_file = new TFile("out.root", "RECREATE");
	TTree* out_tree = new TTree("out_tree", "outtree");
	double riseTime, fallTime, pulseDuration;
	out_tree->Branch("risetime", &riseTime);
	out_tree->Branch("falltime", &fallTime);
	out_tree->Branch("pulseduration", &pulseDuration);
	
	TCanvas* c1 = new TCanvas();
	TGraph* pulse = new TGraph();
	pulse->SetTitle("Output pulse;t [s];U [V]");
	pulse->SetMarkerStyle(7);
	TGraph* rf = new TGraph();	// drawing rise and fall time points
	rf->SetMarkerStyle(8);
	rf->SetMarkerColor(46);
	TF1* bl = new TF1("baseline", "[0]", -100, 100);	// baseline
	bl->SetLineColor(38);
	
        // loop over data
	float uMax, lPos, hPos, lTime, hTime, rTime, buf;
	int maxEntry;
        for (int iEvent = 0; iEvent < nEvents; iEvent++){
                uMax = -1.;
                in_tree->GetEntry(iEvent);
                for (int i = 0; i < sample_len; i++){
                        u[i] *= -1;
                        pulse->SetPoint(i, t[i], u[i]);
                        // find Maximum by Hand because root apparently isnt able to do so
                        if (u[i] > uMax){
                                uMax = u[i];
                                maxEntry = i;
                        }
                }
		// get 10% and 90% amplitude
                lPos = 0.1*(uMax - baseline) + baseline;
                hPos = 0.9*(uMax - baseline) + baseline;
                
		// get rise time -> start at maximum and go left
		lTime = -1;
                hTime = -1;
                for (int i = maxEntry; (buf = pulse->GetY()[i]) > 0.9*lPos; i--){
                        if ( buf > hPos )
                                hTime = pulse->GetX()[i];
                        if ( buf > lPos ){
                                lTime = pulse->GetX()[i];
                        }
                }
                riseTime = hTime - lTime;
		rf->SetPoint(0, lTime, lPos);               
		rf->SetPoint(1, hTime, hPos);               
 
		// get fall time -> start at maximum and go right
		rTime = -1;
                hTime = -1;
                for (int i = maxEntry; (buf = pulse->GetY()[i]) > 0.9*lPos; i++){
                        if ( buf > hPos )
                                hTime = pulse->GetX()[i];
                        if ( buf > lPos ){
                                rTime = pulse->GetX()[i];
                        }
                }
                fallTime = rTime - hTime;
		pulseDuration = rTime - lTime; 
		rf->SetPoint(2, rTime, lPos);               
		rf->SetPoint(3, hTime, hPos);               

		out_tree->Fill();

                // draw & save every 500th event
                if (iEvent%100 == 0) {
 			printf("[Pulse Shape] - Risetime = %e s\n", riseTime);
			printf("[Pulse Shape] - Falltime = %e s\n", pulseDuration);
 			printf("[Pulse Shape] - Pulse duration = %e s\n", fallTime);
			
			bl->SetParameter(0, baseline);
                        pulse->Draw("A*");
			bl->Draw("SAME");
			rf->Draw("SAME*");
			c1->Write();
                }
	}




	// cleanup
	out_tree->Write();
	out_file->Close();
	in_file->Close();
	return 0;
}

//.........这里部分代码省略.........
      texPt->SetNDC();
      texPt->SetTextFont(42);
      texPt->SetTextSize(0.06);
      texPt->SetLineWidth(2);

      TLatex* texY = new TLatex(0.18,0.74,Form("|y| < 1.0"));
      texY->SetNDC();
      texY->SetTextFont(42);
      texY->SetTextSize(0.06);
      texY->SetLineWidth(2);

      c2->cd(1);

      hPtMD0Dca->GetZaxis()->SetRange(1,100);
      hPtMD0Dca->GetXaxis()->SetRangeUser(ptLow+0.001,ptHigh-0.001);
      hPtMD0DcaMCPSignal->GetXaxis()->SetRangeUser(ptLow+0.001,ptHigh-0.001);
      hPtMD0DcaMCPSwapped->GetXaxis()->SetRangeUser(ptLow+0.001,ptHigh-0.001);
      TH1D* hMData = (TH1D*)hPtMD0Dca->Project3D("y")->Clone(Form("hM_%1.1f_%1.1f", ptLow, ptHigh));
      TH1D* hMMCSignal = (TH1D*)hPtMD0DcaMCPSignal->Project3D("y");
      TH1D* hMMCSwapped = (TH1D*)hPtMD0DcaMCPSwapped->Project3D("y");

      setColorTitleLabel(hMData);
      setColorTitleLabel(hMMCSignal);
      setColorTitleLabel(hMMCSwapped);

      TF1* fMass = fitMass(hMData, hMMCSignal, hMMCSwapped);

      texCms->Draw();
      texCol->Draw();
      texPt->Draw();
      texY->Draw();

      TF1* fSignalAndSwapped = new TF1("fSignalAndSwapped","[0]*([3]*([5]*Gaus(x,[1],[2]*(1+[7]))/(sqrt(2*3.1415927)*[2]*(1+[7]))+(1-[5])*Gaus(x,[1],[6]*(1+[7]))/(sqrt(2*3.1415927)*[6]*(1+[7])))+(1-[3])*Gaus(x,[1],[4]*(1+[7]))/(sqrt(2*3.1415927)*[4]*(1+[7])))", 1.7, 2.0);      
      fSignalAndSwapped->SetParameter(0,fMass->GetParameter(0));
      fSignalAndSwapped->SetParameter(1,fMass->GetParameter(1));
      fSignalAndSwapped->SetParameter(2,fMass->GetParameter(2));
      fSignalAndSwapped->SetParameter(3,fMass->GetParameter(7));
      fSignalAndSwapped->SetParameter(4,fMass->GetParameter(8));
      fSignalAndSwapped->SetParameter(5,fMass->GetParameter(9));
      fSignalAndSwapped->SetParameter(6,fMass->GetParameter(10));
      fSignalAndSwapped->SetParameter(7,fMass->GetParameter(11));
  
      TF1* background = new TF1("fBackground","[0]+[1]*x+[2]*x*x+[3]*x*x*x");
      background->SetParameter(0,fMass->GetParameter(3));
      background->SetParameter(1,fMass->GetParameter(4));
      background->SetParameter(2,fMass->GetParameter(5));
      background->SetParameter(3,fMass->GetParameter(6));

      cout<<"MC signal width: "<<fMass->GetParameter(2)<<"   "<<fMass->GetParameter(10)<<endl;
      cout<<"MC swapped width: "<<fMass->GetParameter(8)<<endl;

      float massD = 1.8649;
      float massSignal1 = massD-0.025;
      float massSignal2 = massD+0.025;
      float massSideBand1 = massD-0.1;
      float massSideBand2 = massD-0.075;
      float massSideBand3 = massD+0.075;
      float massSideBand4 = massD+0.1;

      float scaleSideBandBackground = background->Integral(massSignal1, massSignal2)/(background->Integral(massSideBand1, massSideBand2)+background->Integral(massSideBand3, massSideBand4));
      cout<<"scaleSideBandBackground: "<<scaleSideBandBackground<<endl;
      totalYieldInvMassFit[i-1] = fMass->GetParameter(0)*fMass->GetParameter(7)/hMData->GetBinWidth(1);
      totalYieldInvMassFitError[i-1] = fMass->GetParError(0)*fMass->GetParameter(7)/hMData->GetBinWidth(1);
      cout<<"totalYieldInvMassFit: "<<totalYieldInvMassFit[i-1]<<" +- "<<totalYieldInvMassFitError[i-1]<<endl;
      float scaleSideBandMethodSignal = fSignalAndSwapped->GetParameter(0)*fSignalAndSwapped->GetParameter(3) / (fSignalAndSwapped->Integral(massSignal1, massSignal2)-fSignalAndSwapped->Integral(massSideBand1, massSideBand2)-fSignalAndSwapped->Integral(massSideBand3, massSideBand4));
      cout<<"scaleSideBandMethodSignal: "<<scaleSideBandMethodSignal<<endl;

TF1* fit(Double_t ptmin, Double_t ptmax)
{
  TCanvas* c = new TCanvas(Form("c_%.0f_%.0f",ptmin,ptmax),"",600,600);
  TFile* infile = new TFile(Form("%s_%s_%.0f_%.0f.root",infname.Data(),collisionsystem.Data(),ptmin,ptmax));
  TH1D* h = (TH1D*)infile->Get("h");                    h->SetName(Form("h_%.0f_%.0f",ptmin,ptmax));
  TH1D* hMCSignal = (TH1D*)infile->Get("hMCSignal");    hMCSignal->SetName(Form("hMCSignal_%.0f_%.0f",ptmin,ptmax));
  TH1D* hMCSwapped = (TH1D*)infile->Get("hMCSwapped");  hMCSwapped->SetName(Form("hMCSwapped_%.0f_%.0f",ptmin,ptmax));
  TF1* f = new TF1(Form("f_%.0f_%.0f",ptmin,ptmax),"[0]*([7]*([9]*Gaus(x,[1],[2])/(sqrt(2*3.14159)*[2])+(1-[9])*Gaus(x,[1],[10])/(sqrt(2*3.14159)*[10]))+(1-[7])*Gaus(x,[1],[8])/(sqrt(2*3.14159)*[8]))+[3]+[4]*x+[5]*x*x", 1.7, 2.0);

  f->SetParLimits(4,-1000,1000);
  f->SetParLimits(10,0.001,0.05);
  f->SetParLimits(2,0.01,0.1);
  f->SetParLimits(8,0.02,0.2);
  f->SetParLimits(7,0,1);
  f->SetParLimits(9,0,1);
  
  f->SetParameter(0,setparam0);
  f->SetParameter(1,setparam1);
  f->SetParameter(2,setparam2);
  f->SetParameter(10,setparam10);
  f->SetParameter(9,setparam9);

  f->FixParameter(8,setparam8);
  f->FixParameter(7,1);
  f->FixParameter(1,fixparam1);
  f->FixParameter(3,0);
  f->FixParameter(4,0);
  f->FixParameter(5,0);
  h->GetEntries();
  
  hMCSignal->Fit(Form("f_%.0f_%.0f",ptmin,ptmax),"q","",minhisto,maxhisto);
  hMCSignal->Fit(Form("f_%.0f_%.0f",ptmin,ptmax),"q","",minhisto,maxhisto);
  f->ReleaseParameter(1);
  hMCSignal->Fit(Form("f_%.0f_%.0f",ptmin,ptmax),"L q","",minhisto,maxhisto);
  hMCSignal->Fit(Form("f_%.0f_%.0f",ptmin,ptmax),"L q","",minhisto,maxhisto);
  hMCSignal->Fit(Form("f_%.0f_%.0f",ptmin,ptmax),"L m","",minhisto,maxhisto);
  
  f->FixParameter(1,f->GetParameter(1));
  f->FixParameter(2,f->GetParameter(2));
  f->FixParameter(10,f->GetParameter(10));
  f->FixParameter(9,f->GetParameter(9));
  f->FixParameter(7,0);
  f->ReleaseParameter(8);
  f->SetParameter(8,setparam8);
  
  hMCSwapped->Fit(Form("f_%.0f_%.0f",ptmin,ptmax),"L q","",minhisto,maxhisto);
  hMCSwapped->Fit(Form("f_%.0f_%.0f",ptmin,ptmax),"L q","",minhisto,maxhisto);
  hMCSwapped->Fit(Form("f_%.0f_%.0f",ptmin,ptmax),"L q","",minhisto,maxhisto);
  hMCSwapped->Fit(Form("f_%.0f_%.0f",ptmin,ptmax),"L m","",minhisto,maxhisto);
  
  f->FixParameter(7,hMCSignal->Integral(0,1000)/(hMCSwapped->Integral(0,1000)+hMCSignal->Integral(0,1000)));
  f->FixParameter(8,f->GetParameter(8));
  f->ReleaseParameter(3);
  f->ReleaseParameter(4);
  f->ReleaseParameter(5);

  f->SetLineColor(kRed);
  
  h->Fit(Form("f_%.0f_%.0f",ptmin,ptmax),"q","",minhisto,maxhisto);
  h->Fit(Form("f_%.0f_%.0f",ptmin,ptmax),"q","",minhisto,maxhisto);
  f->ReleaseParameter(1);
  //f->ReleaseParameter(2);                                     // you need to release these two parameters if you want to perform studies on the sigma shape
  //f->ReleaseParameter(10);                                   // you need to release these two parameters if you want to perform studies on the sigma shape
  h->Fit(Form("f_%.0f_%.0f",ptmin,ptmax),"L q","",minhisto,maxhisto);
  h->Fit(Form("f_%.0f_%.0f",ptmin,ptmax),"L q","",minhisto,maxhisto);
  h->Fit(Form("f_%.0f_%.0f",ptmin,ptmax),"L q","",minhisto,maxhisto);
  h->Fit(Form("f_%.0f_%.0f",ptmin,ptmax),"L m","",minhisto,maxhisto);
  
  TF1* background = new TF1(Form("background_%.0f_%.0f",ptmin,ptmax),"[0]+[1]*x+[2]*x*x");
  background->SetParameter(0,f->GetParameter(3));
  background->SetParameter(1,f->GetParameter(4));
  background->SetParameter(2,f->GetParameter(5));
  background->SetLineColor(4);
  background->SetRange(minhisto,maxhisto);
  background->SetLineStyle(2);
  
  TF1* mass = new TF1(Form("fmass_%.0f_%.0f",ptmin,ptmax),"[0]*([3]*([4]*Gaus(x,[1],[2])/(sqrt(2*3.14159)*[2])+(1-[4])*Gaus(x,[1],[5])/(sqrt(2*3.14159)*[5])))");
  mass->SetParameters(f->GetParameter(0),f->GetParameter(1),f->GetParameter(2),f->GetParameter(7),f->GetParameter(9),f->GetParameter(10));
  mass->SetParError(0,f->GetParError(0));
  mass->SetParError(1,f->GetParError(1));
  mass->SetParError(2,f->GetParError(2));
  mass->SetParError(3,f->GetParError(7));
  mass->SetParError(4,f->GetParError(9));
  mass->SetParError(5,f->GetParError(10));
  mass->SetFillColor(kOrange-3);
  mass->SetFillStyle(3002);
  mass->SetLineColor(kOrange-3);
  mass->SetLineWidth(3);
  mass->SetLineStyle(2);
  
  TF1* massSwap = new TF1(Form("fmassSwap_%.0f_%.0f",ptmin,ptmax),"[0]*(1-[2])*Gaus(x,[1],[3])/(sqrt(2*3.14159)*[3])");
  massSwap->SetParameters(f->GetParameter(0),f->GetParameter(1),f->GetParameter(7),f->GetParameter(8));
  massSwap->SetParError(0,f->GetParError(0));
  massSwap->SetParError(1,f->GetParError(1));
  massSwap->SetParError(2,f->GetParError(7));
  massSwap->SetParError(3,f->GetParError(8));
  massSwap->SetFillColor(kGreen+4);
  massSwap->SetFillStyle(3005);
  massSwap->SetLineColor(kGreen+4);
  massSwap->SetLineWidth(3);
//.........这里部分代码省略.........

TH1D* getYield(TTree* nt, TTree* ntMC, TString triggerpass, TString triggername, TString prescale, TString variable, TString varname, TString varlatex, Int_t BIN_NUM, Double_t BIN_MIN, Double_t BIN_MAX, TString addcut="")
{
  TH1D* hDistrib = new TH1D(Form("h%s_Distrib_%s",triggername.Data(),varname.Data()),"",BIN_NUM,BIN_MIN,BIN_MAX);
  for(float ivar=0;ivar<BIN_NUM;ivar++)
  {
    TCanvas* c = new TCanvas(Form("c%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"",500,500);
    TH1D* h = new TH1D(Form("h%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),";D^{0} mass (GeV/c^{2});Candidates",60,1.7,2.0);
    TH1D* hMC = new TH1D(Form("hMC%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"",60,1.75,1.95);
    TH1D* hSW = new TH1D(Form("hSW%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"",60,1.75,1.95);
    Float_t varmin = BIN_MIN+ivar*((BIN_MAX-BIN_MIN)/BIN_NUM);
    Float_t varmax = BIN_MIN+(ivar+1)*((BIN_MAX-BIN_MIN)/BIN_NUM);
    nt->Project(Form("h%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),Form("Dmass%s",prescale.Data()),Form("%s%s&&(%s>%f&&%s<%f)%s",prefilter.Data(),addcut.Data(),variable.Data(),varmin,variable.Data(),varmax,triggerpass.Data()));
    ntMC->Project(Form("hMC%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"Dmass",Form("%s%s&&(%s>%f&&%s<%f)%s",prefilterMC.Data(),addcut.Data(),variable.Data(),varmin,variable.Data(),varmax,"&&1"));
    ntMC->Project(Form("hSW%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"Dmass",Form("%s%s&&(%s>%f&&%s<%f)%s",prefilterSW.Data(),addcut.Data(),variable.Data(),varmin,variable.Data(),varmax,"&&1"));
    h->SetMaximum(h->GetMaximum()*1.20);
    h->Draw();

    TF1* f = new TF1(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"[0]*([7]*([9]*Gaus(x,[1],[2]*(1+[11]))/(sqrt(2*3.14159)*[2]*(1+[11]))+(1-[9])*Gaus(x,[1],[10]*(1+[11]))/(sqrt(2*3.14159)*[10]*(1+[11])))+(1-[7])*Gaus(x,[1],[8]*(1+[11]))/(sqrt(2*3.14159)*[8]*(1+[11])))+[3]+[4]*x+[5]*x*x+[6]*x*x*x", 1.7, 2.0);
    f->SetParLimits(4,-1000,1000);
    f->SetParLimits(10,0.001,0.05);
    f->SetParLimits(2,0.01,0.1);
    f->SetParLimits(8,0.02,0.2);
    f->SetParLimits(7,0,1);
    f->SetParLimits(9,0,1);

    f->SetParameter(0,setparam0);
    f->SetParameter(1,setparam1);
    f->SetParameter(2,setparam2);
    f->SetParameter(10,setparam10);
    f->SetParameter(9,setparam9);

    f->FixParameter(8,setparam8);
    f->FixParameter(7,1);
    f->FixParameter(1,fixparam1);
    f->FixParameter(3,0);
    f->FixParameter(4,0);
    f->FixParameter(5,0);
    f->FixParameter(6,0);
    f->FixParameter(11,0);
    h->GetEntries();

    hMC->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"q","",1.7,2.0);
    hMC->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"q","",1.7,2.0);
    f->ReleaseParameter(1);
    hMC->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"L q","",1.7,2.0);
    hMC->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"L q","",1.7,2.0);
    hMC->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"L m","",1.7,2.0);

    f->FixParameter(1,f->GetParameter(1));
    f->FixParameter(2,f->GetParameter(2));
    f->FixParameter(10,f->GetParameter(10));
    f->FixParameter(9,f->GetParameter(9));
    f->FixParameter(7,0);
    f->ReleaseParameter(8);
    f->SetParameter(8,setparam8);

    hSW->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"L q","",1.7,2.0);
    hSW->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"L q","",1.7,2.0);
    hSW->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"L q","",1.7,2.0);
    hSW->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"L m","",1.7,2.0);

    f->FixParameter(7,hMC->Integral(0,1000)/(hSW->Integral(0,1000)+hMC->Integral(0,1000)));
    f->FixParameter(8,f->GetParameter(8));
    f->ReleaseParameter(3);
    f->ReleaseParameter(4);
    f->ReleaseParameter(5);
    f->ReleaseParameter(6);

    f->SetLineColor(kRed);

    h->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"q","",1.7,2.0);
    h->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"q","",1.7,2.0);
    f->ReleaseParameter(1);
    f->SetParLimits(1,1.86,1.87);
    f->ReleaseParameter(11);
    f->SetParLimits(11,-1.,1.);
    h->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"L q","",1.7,2.0);
    h->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"L q","",1.7,2.0);
    h->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"L q","",1.7,2.0);
    h->Fit(Form("f%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"L m","",1.7,2.0);
    h->SetMarkerSize(0.8);
    h->SetMarkerStyle(20);

    TF1* background = new TF1(Form("background%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"[0]+[1]*x+[2]*x*x+[3]*x*x*x");
    background->SetParameter(0,f->GetParameter(3));
    background->SetParameter(1,f->GetParameter(4));
    background->SetParameter(2,f->GetParameter(5));
    background->SetParameter(3,f->GetParameter(6));
    background->SetLineColor(4);
    background->SetRange(1.7,2.0);
    background->SetLineStyle(2);

    TF1* mass = new TF1(Form("fmass%s_Fit_%s_%.0f",triggername.Data(),varname.Data(),ivar),"[0]*([3]*([4]*Gaus(x,[1],[2]*(1+[6]))/(sqrt(2*3.14159)*[2]*(1+[6]))+(1-[4])*Gaus(x,[1],[5]*(1+[6]))/(sqrt(2*3.14159)*[5]*(1+[6]))))");
    mass->SetParameters(f->GetParameter(0),f->GetParameter(1),f->GetParameter(2),f->GetParameter(7),f->GetParameter(9),f->GetParameter(10),f->GetParameter(11));
    mass->SetParError(0,f->GetParError(0));
    mass->SetParError(1,f->GetParError(1));
    mass->SetParError(2,f->GetParError(2));
    mass->SetParError(3,f->GetParError(7));
    mass->SetParError(4,f->GetParError(9));
    mass->SetParError(5,f->GetParError(10));
//.........这里部分代码省略.........

void compareCorrRCSpike(){
	

  char tmp[1000];
  setTDRStyle();
  // settings for purity TGraphAsymmetryErrors
  Double_t fBinsPtMidPoint[nPtBin]={0};
  Double_t fBinsPtError[nPtBin]={0};
  
  for(int ipt=0; ipt < nPtBin; ipt++)
    {
      fBinsPtMidPoint[ipt] = 0.5*(fBinsPt[ipt+1]+fBinsPt[ipt]);
      fBinsPtError[ipt] = 0.5*(fBinsPt[ipt+1]-fBinsPt[ipt]);
    }


  TH1F* hTemplate_S[nEtaBin][nPtBin];
  TH1F* hdata_S[nEtaBin][nPtBin];

  
  std::string dataFile     = "SBDataTemplate_131511_139239.root";
  std::string templateFile = "spike_131511_139239.root";

  TFile* inf_data = new TFile(dataFile.data());
  TFile* inf_template = new TFile(templateFile.data());

  
  char* dec[2] = {"EB","EE"};
  for(int ieta=0; ieta<1; ieta++){
    for(int ipt=0; ipt < nPtBin; ipt++){

      // getting histograms from data root file

      sprintf(tmp,"h_%s_comb3Iso_sig_pt%d",dec[ieta],(int)fBinsPt[ipt]);
      cout << "looking for histogram " << tmp << " in file " << 
	inf_data->GetName() << endl;
      hdata_S[ieta][ipt] = (TH1F*)inf_data->FindObjectAny(tmp);
      hdata_S[ieta][ipt]->SetName(Form("hdata_S_%d_%d",ieta,ipt));
      hdata_S[ieta][ipt]->Rebin(REBINNINGS);

      // getting histogram for template root file
      sprintf(tmp,"h_EB_comb3Iso_EGdata_SIG");
      cout << "looking for histogram " << tmp << " in file " << 
	inf_template->GetName() << endl;
      hTemplate_S[ieta][ipt] = (TH1F*)inf_template->FindObjectAny(tmp);
      hTemplate_S[ieta][ipt]->SetName(Form("hTemplate_S_%d_%d",ieta,ipt));
      hTemplate_S[ieta][ipt]->Rebin(REBINNINGS);

    }
  }


  TH1F* hfit_sig;
  TH1F* hfit_spike;

  TH1F* hsig = new TH1F("hsig","",120,-1,11);
  hsig->SetXTitle("Iso (GeV)");
  hsig->SetYTitle("A.U.");
  hsig->GetYaxis()->SetDecimals();
  hsig->SetYTitle("Iso (GeV)");

  hsig->SetLineColor(2);
  TH1F* hspike = new TH1F("hspike","",120,-1,11);
  hspike->SetLineColor(4);

  TCanvas* c1 = new TCanvas("c1","",500,500);

  for(int ieta=0; ieta< 1; ieta++){
    for(int ipt=0; ipt < nPtBin; ipt++){


      hsig->Reset();
      hspike->Reset();

      // first obtain fit for the random cone corrected values
      hfit_sig  = (TH1F*)hdata_S[ieta][ipt]->Clone();
      hfit_sig  -> SetName("hfit_sig");
      hfit_sig  -> Rebin(4);

      double sigPar[20] = {hfit_sig->GetMaximum(), 1., 0.5, 0.3,
			   1.,-.3,-1., 0.01, 0.5, 0.01, 1., 1.};

      TF1 *fsig = new TF1("fsig", exp_conv, -1., 11., 12);
      fsig->SetParameters(sigPar);
      fsig->SetNpx(2500);
      fsig->SetLineColor(2);
      hfit_sig->Fit(fsig,"","",-1,5.0);

      c1->Print(Form("hfit_sig_pt%d.gif",(int)fBinsPt[ipt]));
      
      fsig->SetParameter(0,2);
      fsig->SetParameter(1,fsig->GetParameter(1)*8.39614e-01/6.83080e-01);//correction from RC
      fsig->SetParameter(2,4.83182e-01);
      fsig->SetParameter(3,fsig->GetParameter(3)*2.33769e-01/2.26323e-01);
    
      cout << "fsig Integral = " << fsig->Integral(-1,11) << endl;
      for(int i=0;i<4;i++)
	cout << "fsig par " << i << "= " << fsig->GetParameter(i) << endl;

      // second obtain fit for the spikes
//.........这里部分代码省略.........