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

//.........这里部分代码省略.........
	kSProbab_ = 0.; jBProbab_ = 0.;
	kSValue_ = 0.; jBValue_ = 0.; 
	aDValue_= 0.; aDProbab_ = 0.;
	nBin_ = 0.; xMin_ = 0.; xMax_ = 0.;
	
	// basic info
	detid_ = *detid;
	fedKey_ = fedKey;
	fecCrate_ = *fecCrate;
	fecSlot_ = *fecSlot;
	fecRing_ = *fecRing;
	ccuAdd_  = *ccuAdd;
	ccuChan_ = *ccuChan;
	lldChannel_ = *lldChannel;
	fedId_   = *fedId;
	fedCh_   = *fedCh;
	apvId_   = apvID;
	stripId_ = stripID;
	
	// basic info of nioise distribution
	noiseMean_ = histoNoiseStrip->GetMean();
	noiseRMS_  = histoNoiseStrip->GetRMS();
	noiseSkewness_ = histoNoiseStrip->GetSkewness();
	noiseKurtosis_ = histoNoiseStrip->GetKurtosis();
	float integral = histoNoiseStrip->Integral();	
	noiseIntegral3Sigma_ = (histoNoiseStrip->Integral(histoNoiseStrip->FindBin(noiseMean_+noiseRMS_*3),histoNoiseStrip->GetNbinsX()+1) + histoNoiseStrip->Integral(0,histoNoiseStrip->FindBin(noiseMean_-noiseRMS_*3)))/integral;
	noiseIntegral4Sigma_ = (histoNoiseStrip->Integral(histoNoiseStrip->FindBin(noiseMean_+noiseRMS_*4),histoNoiseStrip->GetNbinsX()+1) + histoNoiseStrip->Integral(0,histoNoiseStrip->FindBin(noiseMean_-noiseRMS_*4)))/integral;
	noiseIntegral5Sigma_ = (histoNoiseStrip->Integral(histoNoiseStrip->FindBin(noiseMean_+noiseRMS_*5),histoNoiseStrip->GetNbinsX()+1) + histoNoiseStrip->Integral(0,histoNoiseStrip->FindBin(noiseMean_-noiseRMS_*5)))/integral;
	
	// make a gaussian fit	  	
	if(fitFunc == NULL or fitFunc == 0){
	  fitFunc = new TF1 ("fitFunc","gaus(0)",histoNoise->GetXaxis()->GetXmin(),histoNoise->GetXaxis()->GetXmax());
	}
	fitFunc->SetRange(histoNoise->GetXaxis()->GetXmin(),histoNoise->GetXaxis()->GetXmax());
	fitFunc->SetParameters(histoNoiseStrip->Integral(),noiseMean_,noiseRMS_);
	result = histoNoiseStrip->Fit(fitFunc,"QSR");

	if(result.Get()){
	    fitStatus_     = result->Status();
	    fitGausNormalization_  = fitFunc->GetParameter(0);
	    fitGausMean_   = fitFunc->GetParameter(1);
	    fitGausSigma_  = fitFunc->GetParameter(2);
	    fitGausNormalizationError_  = fitFunc->GetParError(0);
	    fitGausMeanError_  = fitFunc->GetParError(1);
	    fitGausSigmaError_ = fitFunc->GetParError(2);
	    fitChi2_           = result->Chi2();
	    fitChi2Probab_     = result->Prob();

	    noiseIntegral3SigmaFromFit_ = (histoNoiseStrip->Integral(histoNoiseStrip->FindBin(noiseMean_+fitGausSigma_*3),histoNoiseStrip->GetNbinsX()+1) + histoNoiseStrip->Integral(0,histoNoiseStrip->FindBin(noiseMean_-fitGausSigma_*3)))/histoNoiseStrip->Integral();
	    noiseIntegral4SigmaFromFit_ = (histoNoiseStrip->Integral(histoNoiseStrip->FindBin(noiseMean_+fitGausSigma_*4),histoNoiseStrip->GetNbinsX()+1) + histoNoiseStrip->Integral(0,histoNoiseStrip->FindBin(noiseMean_-fitGausSigma_*4)))/histoNoiseStrip->Integral();
	    noiseIntegral5SigmaFromFit_ = (histoNoiseStrip->Integral(histoNoiseStrip->FindBin(noiseMean_+fitGausSigma_*5),histoNoiseStrip->GetNbinsX()+1) + histoNoiseStrip->Integral(0,histoNoiseStrip->FindBin(noiseMean_-fitGausSigma_*5)))/histoNoiseStrip->Integral();
	    
	    jBValue_   = (histoNoiseStrip->Integral()/6)*(noiseSkewness_*noiseSkewness_+(noiseKurtosis_*noiseKurtosis_)/4);	  
	    jBProbab_  = ROOT::Math::chisquared_cdf_c(jBValue_,2);

	    if(randomHisto == 0 or randomHisto == NULL)
	      randomHisto = (TH1F*) histoNoiseStrip->Clone("randomHisto");	    	    
	    randomHisto->Reset();
	    randomHisto->SetDirectory(0);     
	
      
	    if(integral != 0){	      
	      if(generateRandomDistribution){
		randomHisto->FillRandom("fitFunc",histoNoiseStrip->Integral());	    
		kSValue_  = histoNoiseStrip->KolmogorovTest(randomHisto,"MN");
		kSProbab_ = histoNoiseStrip->KolmogorovTest(randomHisto,"N");	    

//-------------------------------------------------------------------
void Alignment_1(int cluster_length_leftCut = 2,int cluster_length_rightCut = 7,float TotSignal_Cut = 5){
  gStyle->SetOptStat(11111111);
  gStyle->SetOptFit(111111111);
  //gAnaCombine->LinkBranches();
  ofstream outputlog("./AMS/AMSAlignment_1.txt",ios::out | ios::app);
  outputlog<<"\nfile: "<<gRootIOSvc->GetInputFileName()<<".\tCluster length: "<<cluster_length_leftCut<<"~"<<cluster_length_rightCut<<".\tTotal signale cut "<<TotSignal_Cut<<endl;

  TH1D *h_ladder_stripID[5][2];
  TH1D *h_Refer[5][2];
  char name[80];
  for(short lid =0;lid<5;++lid){
    for(short sid=0;sid<2;++sid){
      snprintf(name,80,"ladder_%d-side_%d",lid,sid);
      h_ladder_stripID[lid][sid] = new TH1D(name,name,500,0,1100);
      snprintf(name,80,"ladder_%d-side_%d VS ladder 0",lid,sid);
      h_Refer[lid][sid] = new TH1D(name,name,200,-200,200);
    }
    snprintf(name,80,"ladder_%d",lid);
    h_ladder_stripID[lid][0]->SetTitle(name);
    snprintf(name,80,"Seed Address ladder_%d",lid);
    h_Refer[lid][0]->SetTitle(name);
  }

  TChain *tree = gAnaCombine->GetTree();
  long nEvt = tree->GetEntries();
  std::vector<Cluster*> clusterIndex[5][2];
  for(long ievt=0;ievt<nEvt;++ievt){
    tree->GetEntry(ievt);
    short nCluster = gAnaCombine->fEvtAMS->GetEntriesFast();
    for(int iclu=0;iclu<nCluster;++iclu){
      Cluster *acluster = (Cluster*)gAnaCombine->fEvtAMS->At(iclu);
      if(acluster->GetTotSig()>TotSignal_Cut && cluster_length_leftCut<=acluster->length && acluster->length <= cluster_length_rightCut){
        int ladder = acluster->ladder;
        if(ladder==5){
          ladder = 2;
        }else if(ladder == 4){
          ladder = 3;
        }else if(ladder == 3){
          ladder = 4;
        }
        clusterIndex[ladder][acluster->side].push_back(acluster);
      }
    }
    if(clusterIndex[0][0].size() == 1 && clusterIndex[0][1].size() == 1){
      for(short il=0;il<5;++il){
        for(short is=0;is<2;++is){
          if(clusterIndex[il][is].size() == 1){
            h_ladder_stripID[il][is]->Fill(clusterIndex[il][is][0]->GetSeedAdd());
            h_Refer[il][is]->Fill(clusterIndex[il][is][0]->GetSeedAdd() - clusterIndex[0][is][0]->GetSeedAdd());
          }
          clusterIndex[il][is].clear();
        }
      }
    }
  }
  snprintf(name,80,"AMS Alignment_1_a: cluster length %d~%d, total signal cut %f",cluster_length_leftCut,cluster_length_rightCut,TotSignal_Cut);
  TCanvas *c1 = new TCanvas(name,name);
  c1->Divide(2,3);
  snprintf(name,80,"AMS Alignment_1_b: cluster length %d~%d, totalt signal cut %f",cluster_length_leftCut,cluster_length_rightCut,TotSignal_Cut);

  TF1 *gausFit = new TF1("GausFit","gaus",0,150);
  TCanvas *c2 = new TCanvas(name,name);
  c2->Divide(2,3);
  double mean=0.,sigma =0.;
  for(short lid=0;lid<5;++lid){
    c1->cd(lid+1);
    h_ladder_stripID[lid][0]->Draw();
    h_ladder_stripID[lid][1]->SetLineColor(6);
    h_ladder_stripID[lid][1]->Draw("same");

    c2->cd(lid+1);
    h_Refer[lid][0]->Draw();
    mean = h_Refer[lid][0]->GetMean(); sigma = h_Refer[lid][0]->GetRMS();
    gausFit->SetRange(mean-2.5*sigma,mean+2.5*sigma);
    h_Refer[lid][0]->Fit(gausFit,"R0Q");
    gausFit->DrawCopy("lsame");
    outputlog<<lid<<"  "<<gausFit->GetParameter(1)<<"  "<<gausFit->GetParameter(2)<<"    ";

    h_Refer[lid][1]->SetLineColor(6);
    h_Refer[lid][1]->Draw("same");
    mean = h_Refer[lid][1]->GetMean(); sigma = h_Refer[lid][1]->GetRMS();
    gausFit->SetRange(mean-2.5*sigma,mean+2.5*sigma);
    h_Refer[lid][1]->Fit(gausFit,"R0Q");
    gausFit->DrawCopy("lsame");
    outputlog<<gausFit->GetParameter(1)<<"  "<<gausFit->GetParameter(2)<<endl;
  }
}

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");   
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  f->ReleaseParameter(6);

  f->SetLineColor(kRed);
  

  /*
  TCanvas*mycanvas=new TCanvas("mycanvas","mycanvas",1000,500);
  mycanvas->Divide(2,1);
  mycanvas->cd(1);
  hMCSignal->Draw();
  mycanvas->cd(2);
  hMCSwapped->Draw();
  mycanvas->SaveAs("mycanvas.pdf");
*/
  
  h->Fit(Form("f%d",count),"q","",minhisto,maxhisto);
  h->Fit(Form("f%d",count),"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%d",count),"L q","",minhisto,maxhisto);
  h->Fit(Form("f%d",count),"L q","",minhisto,maxhisto);
  h->Fit(Form("f%d",count),"L q","",minhisto,maxhisto);
  h->Fit(Form("f%d",count),"L m","",minhisto,maxhisto);
  
  std::cout<<"parameter 10="<<f->GetParameter(10)<<std::endl;

  TF1* background = new TF1(Form("background%d",count),"[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(Form("fmass%d",count),"[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%d",count),"[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);
  massSwap->SetLineStyle(1);
  
  h->SetXTitle("m_{#piK} (GeV/c^{2})");
  h->SetYTitle("Entries / (5 MeV/c^{2})");
  h->GetXaxis()->CenterTitle();
  h->GetYaxis()->CenterTitle();

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));
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  Float_t binsY[nBinY+1] = {-0.0734,-0.0562,-0.0428,-0.0320,-0.0236,-0.0170,-0.0118,-0.0078,-0.0046,-0.002,0.0,0.002,0.0046,0.0078,0.0118,0.0170,0.0236,0.0320,0.0428,0.0562,0.0734};



  for(int i=1; i<=nPtBins; i++)
    {
      pts[i-1] = 0.5*(ptBins[i-1]+ptBins[i]);
      ptErrors[i-1] = 0.5*(ptBins[i]-ptBins[i-1]);
      float ptLow = ptBins[i-1];
      float ptHigh = ptBins[i];
      cout<<endl<<"======================================="<<endl;
      cout<<"pT range: "<<ptLow<<" "<<ptHigh<<endl;

      TLatex* texPtY = new TLatex(0.32,0.82,Form("%.1f < p_{T} < %.1f GeV/c      |y| < 1.0",ptLow,ptHigh));
      texPtY->SetNDC();
      texPtY->SetTextFont(42);
      texPtY->SetTextSize(0.06);
      texPtY->SetLineWidth(2);

      TLatex* texPt = new TLatex(0.18,0.82,Form("%.1f < p_{T} < %.1f GeV/c",ptLow,ptHigh));
      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));

TF1 *fit(TTree *nt, TTree *ntMC, Double_t ptmin, Double_t ptmax, int isMC,bool isPbPb,TF1* &total,Float_t centmin, Float_t centmax)
{
	//cout<<cut.Data()<<endl;
	static Int_t count=0;
	count++;
	TCanvas* c= new TCanvas(Form("c%d",count),"",600,600);
	TH1D* h = new TH1D(Form("h%d",count),"",nbinsmasshisto,minhisto,maxhisto);

	//   TF1* f = new TF1(Form("f%d",count),"[0]*Gaus(x,[1],[2])/(sqrt(2*3.14159)*[2]) + [3]*Gaus(x,[4],[5])/(sqrt(2*3.14159)*[5]) + [6]*Gaus(x,[7],[8])/(sqrt(2*3.14159)*[8]) + [9]+[10]*x ");
	//   TF1* f = new TF1(Form("f%d",count),"[0]*Gaus(x,[1],[2])/(sqrt(2*3.14159)*[2]) + [3]*Gaus(x,[4],[5])/(sqrt(2*3.14159)*[5]) + [9]+[10]*x ");
	TF1* f = new TF1(Form("f%d",count),"[0]*TMath::Erf((x-[1])/[2]) + [0] + [9]+[10]*x ");

	if(isMC==1) nt->Project(Form("h%d",count),"Bmass",Form("%s*(%s&&Bpt>%f&&Bpt<%f)","1",seldata.Data(),ptmin,ptmax));   
	else nt->Project(Form("h%d",count),"Bmass",Form("(%s&&Bpt>%f&&Bpt<%f)",seldata.Data(),ptmin,ptmax));   

	clean0(h);

	h->Draw();
	f->SetParLimits(0, 1e-2, 1e4);
	f->SetParLimits(1, 5.02, 5.06);
	f->SetParLimits(2, 0.001, 0.1);

	f->SetParLimits(3, 1e-2, 1e4);
	f->SetParLimits(4, 5.06, 5.10);
	f->SetParLimits(5, 0.001, 0.1);

	f->SetParLimits(6, 0, 1e4);
	f->SetParLimits(7, 5.3, 5.4);
	f->SetParLimits(8, 0.001, 0.5);

	f->SetParLimits(9, 0, 1e5);
	f->SetParLimits(10, -500,  100);

	f->SetParameter(0,1e2);
	f->SetParameter(1,5.03);
	f->SetParameter(2,0.05);

	f->SetParameter(0,1e2);
	f->SetParameter(1,5.03);
	f->SetParameter(2,0.05);

	f->SetParameter(3,1e2);
	f->SetParameter(4,5.07);
	f->SetParameter(5,0.05);

	f->SetParameter(6,1e2);
	f->SetParameter(7,5.35);
	f->SetParameter(8,0.05);

	f->SetParameter(9,1e3);
	f->SetParameter(10,-1);

	//error fn
	f->SetParLimits(0, 1e1, 1e3);
	f->SetParLimits(1, 5., 5.3);
	f->SetParLimits(2, -10, 0);
	f->SetParameter(0,100);
	f->SetParameter(1,5.1);
	f->SetParameter(2,0);

	h->GetEntries();

	h->Fit(Form("f%d",count),"q","",minhisto,maxhisto);
	h->Fit(Form("f%d",count),"q","",minhisto,maxhisto);
	h->Fit(Form("f%d",count),"L q","",minhisto,maxhisto);
	h->Fit(Form("f%d",count),"L q","",minhisto,maxhisto);
	h->Fit(Form("f%d",count),"L q","",minhisto,maxhisto);
	h->Fit(Form("f%d",count),"L m","",minhisto,maxhisto);
	h->SetMarkerSize(0.8);
	h->SetMarkerStyle(20);

	TF1 *background = new TF1(Form("background%d",count),"[0]+[1]*x");
	background->SetParameter(0,f->GetParameter(9));
	background->SetParameter(1,f->GetParameter(10));
	background->SetLineColor(4);
	background->SetRange(minhisto,maxhisto);
	background->SetLineStyle(2);

	//   TF1 *mass = new TF1(Form("fmass%d",count),"[0]*Gaus(x,[1],[2])/(sqrt(2*3.14159)*[2]) + [3]*Gaus(x,[4],[5])/(sqrt(2*3.14159)*[5]) + [6]*Gaus(x,[7],[8])/(sqrt(2*3.14159)*[8])");
	//   mass->SetParameters(f->GetParameter(0),f->GetParameter(1),f->GetParameter(2),f->GetParameter(3),f->GetParameter(4),f->GetParameter(5),f->GetParameter(6),f->GetParameter(7),f->GetParameter(8));
	//TF1 *mass = new TF1(Form("fmass%d",count),"[0]*Gaus(x,[1],[2])/(sqrt(2*3.14159)*[2]) + [3]*Gaus(x,[4],[5])/(sqrt(2*3.14159)*[5])");
	//mass->SetParameters(f->GetParameter(0),f->GetParameter(1),f->GetParameter(2),f->GetParameter(3),f->GetParameter(4),f->GetParameter(5));
	TF1 *mass = new TF1(Form("fmass%d",count),"[0]*TMath::Erf((x-[1])/[2]) + [0]");
	mass->SetParameters(f->GetParameter(0),f->GetParameter(1),f->GetParameter(2));
	//mass->SetParError(0,f->GetParError(0));
	//mass->SetParError(1,f->GetParError(1));
	//mass->SetParError(2,f->GetParError(2));
	mass->SetLineColor(2);

	h->SetXTitle("m_{#mu#muK} (GeV/c^{2})");
	h->SetYTitle("Entries / (5 MeV/c^{2})");
	h->GetXaxis()->CenterTitle();
	h->GetYaxis()->CenterTitle();
	h->SetAxisRange(0,h->GetMaximum()*1.4*1.2,"Y");
	h->GetXaxis()->SetTitleOffset(1.3);
	h->GetYaxis()->SetTitleOffset(1.8);
	h->GetXaxis()->SetLabelOffset(0.007);
	h->GetYaxis()->SetLabelOffset(0.007);
	h->GetXaxis()->SetTitleSize(0.045);
	h->GetYaxis()->SetTitleSize(0.045);
//.........这里部分代码省略.........

void ClusterWidthAnalysisTreeMaker::FitProfiles(std::map<ULong64_t , std::vector<TProfile*> > &ProfVsAngle, string output_file){

  TFile * myFile = new TFile(output_file.c_str(), "recreate");

  ULong64_t detid;
  double voltage;
  double errvoltage;
  double Slope;
  double errSlope;
  double Origin;
  double errOrigin;
  double Chi2;
  int index;
  TTree *tree = new TTree("T", "summary information");

  tree->Branch("DetID",&detid, "DetID/l");
  tree->Branch("Voltage",&voltage,"Voltage/D");
  tree->Branch("Index",&index,"Index/I");
  tree->Branch("errVoltage",&errvoltage,"errVoltage/D");
  tree->Branch("Slope",&Slope,"Slope/D");
  tree->Branch("errSlope",&errSlope,"errSlope/D");
  tree->Branch("Origin",&Origin,"Origin/D");
  tree->Branch("errOrigin",&errOrigin,"errOrigin/D");
  tree->Branch("Chi2",&Chi2,"Chi2/D");


  //TCanvas* c1 = new TCanvas();
  TH1F* hChi2 = new TH1F("hChi2", "hChi2", 100, 0, 100);
  
  unsigned int nfitrm=0;

  for(std::map<ULong64_t , std::vector<TProfile*> >::iterator iter = ProfVsAngle.begin(); iter != ProfVsAngle.end(); ++iter){
    
	unsigned int i=0; // voltage index    
    std::set< int >::iterator itVolt;
	std::set< int > Voltage = VSmaker.getVoltageList();
    for( itVolt=Voltage.begin(); itVolt!=Voltage.end(); itVolt++){
      
      //std::cout<<"going through the measurement: " << i << std::endl;
            
      TString thestring;
      thestring.Form("DetID_%llu_prof_%u",iter->first,i);
 	  	  
	  if(i>=iter->second.size()) 
       { std::cout<<" Wrong number of voltage steps. "<<std::endl; i++; continue;}
 	  TProfile*  Prof = iter->second[i];
	  
	  if(!Prof) 
       { std::cout<<" Profile "<<thestring.Data()<<"_"<<i<<" not found."<<std::endl; i++; continue;}
 
      //if(Histo->GetEntries()) hNhits->Fill(Histo->Integral());
	  
	  /*if(Histo->Integral()<20) //0.1
	   { //std::cout<<" Not enought entries for histo "<<thestring.Data()<<std::endl;
	    i++; continue;}*/
 
	  
	  detid = iter->first;

	  bool rmfit=false;
	  
	  TF1* pol = new TF1("pol", "pol1", -3, 3);
	  pol->SetRange(-1,0);
	  Prof->Fit("pol", "qr");
	  double chi2 = pol->GetChisquare()/pol->GetNDF();
	  hChi2->Fill(chi2);
	  if(chi2>10) rmfit=true;
	  

      if( rmfit || 
      // TIB modules
          // TIB - 1.4.2.5
      detid==369121605 || detid==369121606 || detid==369121614 || 
      detid==369121613 || detid==369121610 || detid==369121609 ||
          // TIB - 1.2.2.1
      detid==369121390 || detid==369121382 || detid==369121386 || 
      detid==369121385 || detid==369121389 || detid==369121381 ||
          // others in TIB  
      detid==369121437 || detid==369142077 || detid==369121722 || 
      detid==369125534 || detid==369137018 || detid==369121689 ||
      detid==369121765 || detid==369137045 || detid==369169740 ||
      detid==369121689 ||
      // TOB modules 
	      // TOB + 4.3.3.8
      detid/10==436281512 || detid/10==436281528 || detid/10==436281508 ||
      detid/10==436281524 || detid/10==436281520 || detid/10==436281516 ||
          // others in TOB  
      detid/10==436228249 || detid/10==436232694 || detid/10==436228805 ||
      detid/10==436244722 || detid/10==436245110 || detid/10==436249546 ||
      detid/10==436310808 || detid/10==436312136 || detid/10==436315600 ||
	      // without 'sensors' option 
      detid==436281512 || detid==436281528 || detid==436281508 ||
      detid==436281524 || detid==436281520 || detid==436281516 ||
      detid==436228249 || detid==436232694 || detid==436228805 ||
      detid==436244722 || detid==436245110 || detid==436249546 ||
      detid==436310808 || detid==436312136 || detid==436315600 || 
      // TID modules
      detid==402664070 || detid==402664110 ||
	  // TEC modules in small scans
      detid==470148196 || detid==470148200 || detid==470148204 ||
//.........这里部分代码省略.........

void nudgeVpd( string vpdParams, string qaFile, double fitLow = -10.0, double fitHigh = 10.0 ) {

    Reporter rp( "VpdNudge.pdf" );
    rp.newPage();

    ifstream infile;
    infile.open( vpdParams.c_str() );

    ofstream outfile;
    outfile.open( "vpd_4DB.nudged.dat" );

    TFile * qa = new TFile( qaFile.c_str() );

    TH1D * dz = (TH1D*) qa->Get( "zvertexDelta" );
    dz->GetXaxis()->SetRangeUser( fitLow * 1.5, fitHigh * 1.5 );

    TF1 *g = new TF1( "gaus", "gaus" );
    g->SetRange( fitLow, fitHigh );

    dz->Draw();
    gPad->SetLogy(1);
    dz->Fit( g, "R" );
    rp.savePage();

    double mean = g->GetParameter( 1 );
    cout << "Adjust mean by " << -mean << " [cm]" << endl;
    const Float_t c_light = 29.9792458;
    double dt = 2 * (mean) / c_light;
    cout << "dt = " << dt << " [ns]" << endl;


    /**
     * Read in and write out new
     */

    for ( int i = 0; i < 38; i++ ) {

        int channel = -1, nBins = -1;
        infile >> channel;
        infile >> nBins;

        outfile << channel << endl;
        outfile << nBins << endl;

        //copy the bin edges
        for ( int iBin = 0; iBin <= nBins; iBin ++ ) {

            double bv = -1;
            infile >> bv;
            outfile << bv << " ";
        }
        outfile << endl;

        // copy the corrections
        for ( int iBin = 0; iBin <= nBins; iBin ++ ) {

            double bv = -1;
            infile >> bv;
            if ( channel < 20 ) // only apply correction to one side
                outfile << (bv + dt) << " ";
            else
                outfile << bv << " ";
        }
        outfile << endl;

    } // loop on channels




    infile.close();
    outfile.close();


    delete qa;
    delete g;



}

void PlotPhCal(const char *filename = "phCalibration_C0.dat", int col, int row)
{

  Init();

  int value, a, b, n;
  double x[10], y[10];
  char string[200];
  double aovergain, gain, par1;


  TF1 *pol2Fit = new TF1("pol2Fit", "pol2");
  pol2Fit->SetRange(50,500);
  TF1 *linFit = new TF1("linFit", "pol1");
  linFit->SetRange(150,600);
  TF1 *ursFit = new TF1("ursFit", Fitfcn, 40., 1500., 4);
  ursFit->SetParameter(0,0.00382);
  ursFit->SetParameter(1,0.886);
  ursFit->SetParameter(2,112.7);
  ursFit->SetParameter(3,113.0);


  TGraphErrors *graph = new TGraphErrors();
  TGraph *gSlope = new TGraph();

  FILE *inputFile = fopen(filename, "r");

  for (int i = 0; i < 4; i++) fgets(string, 200, inputFile);
  
  for (int icol = 0; icol < 52; ++icol)
    {
      for (int irow = 0; irow < 80; ++irow)
	{
	  n = 0;
	  for (int point = 0; point < 10; point++)
	    {
	      fscanf(inputFile, "%s", string);
	      
	      if (strcmp(string, "N/A") == 0)  value = 7777;
	      else 
		{
		  value  = atoi(string);
		}		      
	      if ((icol == col) && (irow == row))
		{

		  if (point == 0) x[n] = 50;
		  else if (point == 1) x[n] = 100;
		  else if (point == 2) x[n] = 150;
		  else if (point == 3) x[n] = 200;
		  else if (point == 4) x[n] = 250;
		  else if (point == 5) x[n] = 210;
		  else if (point == 6) x[n] = 350;
		  else if (point == 7) x[n] = 490;
		  else if (point == 8) x[n] = 560;
		  else if (point == 9) x[n] = 1400;
		  
		  if (value != 7777)
		    {
		      //  y[n] = value;  // in testboard adc range
		      y[n] = (value + 400) / 7;  // in 8 bit experiment adc range
		      graph->SetPoint(n, x[n], y[n]);	
		      //  graph->SetPointError(n, 2, 10);   // in testboard adc range
		      graph->SetPointError(n, 2, 10/6);   // in 8 bit experiment adc range

		      n++;
		    }
		}
	    }
	  fscanf(inputFile, "%s %2i %2i", string, &a, &b);  //comment
	}
    }
  graph->GetXaxis()->SetTitle("Vcal [low range DAC units]");
  graph->GetYaxis()->SetTitle("PH [ADC units]");
  graph->GetYaxis()->SetTitleOffset(1.2);


  graph->GetXaxis()->SetTitleSize(0.055);
  graph->GetYaxis()->SetTitleSize(0.055);
  graph->GetXaxis()->SetLabelSize(0.04);
  graph->GetYaxis()->SetLabelSize(0.05);
  graph->GetXaxis()->SetTitleOffset(1.15);
  graph->GetYaxis()->SetTitleOffset(1.2);
  graph->SetTitle("");
  graph->SetLineStyle(1);
  graph->SetLineWidth(2);

  //   graph->Fit("linFit","R");
//   gain = linFit->GetParameter(1);


  graph->Fit("ursFit","R");
  graph->Draw("A*");
  par1 = ursFit->GetParameter(1);

  
  NewLatex()->DrawLatex(0.55, 0.5, Form("p_{1} = %.1f ", par1));
  //  NewLatex()->DrawLatex(0.45, 0.4, "NON-LINEAR pixel");       // PlotPhCal("/home/l_tester/ptr/moduleDB/M0567-070115.08:23/T-10a/phCalibration_C12.dat",7,23)
  NewLatex()->DrawLatex(0.5, 0.4, "LINEAR pixel");                // PlotPhCal("/home/l_tester/ptr/moduleDB/M0493-061206.10:45/T-10a/phCalibration_C1.dat",7,23)
}

//.........这里部分代码省略.........
    }
    else if ((samples[is]=="zmmjet" || samples[is]=="zeejet") &&
	     methods[im]=="ptchs" && nozjptb) {
      tex->SetTextColor(kRed);
      tex->DrawLatex(0.20,0.63,"Z+jet p_{T}^{bal} excl. from fit");
      tex->SetTextColor(kBlack);
    }
    else if (samples[is]=="multijet") {
      g->SetMarkerColor(kGray+1);
      g->SetLineColor(kGray+1);
      if (methods[im]=="ptchs") gmjb = g;
      if (methods[im]=="mpfchs1") gmpf = g;
    }
    else
      mg->Add(g);
  } // for is
  } // for im
  
  if (nogjmpf) s += "_nogjmpf";
  if (nogjptb) s += "_nogptb";
  if (nozjptb) s += "_nozptb";
  if (mjvsjes) {
    s += "_mjvsjes";
    tex->SetTextColor(kBlack);
    tex->DrawLatex(0.20,0.58,"Multijet vs JES fit");
  }

  TF1 *fjes = new TF1("fjes",jesFit,30,2200,2);
  fjes->SetParameters(0.99,0.05);
  mg->Fit(fjes,"RN");
  fjes->SetLineColor(kBlack);
  fjes->SetLineStyle(kDashed);
  fjes->SetLineWidth(2);
  fjes->SetRange(10.,3500.);
  fjes->Draw("SAME");
  
  //TF1 *ft = new TF1("ft","1-[0]-[1]*pow(x,[2]) + ([3]+[4]*log(x))/x",30,2200);
  //ft->SetParameters(0,0.05,-0.5,1,0.1);
  //ft->FixParameter(3,0);

  // Logarithmic sigmoid
  //TF1 *ft = new TF1("ft","[0]+(1-[0])/(1. + exp(-(log(x)-log(abs([1])))"
  //	       "/(log(abs([2])+abs([1]))-log(abs([1])))))", 30,2200);
  //ft->SetParameters(0.98, 150, 50);
  TF1 *ft = new TF1("ft","[0]+(1-[0])/(1. + exp(-(log(x)-[1])/[2]))",30,2200);
  //ft->SetParameters(0.98,log(145),log(190)-log(145));
  //ft->SetParameters(0.982,4.967,0.271);
  //ft->SetParameters(0.976,5.040,0.370); // ENDCAP
  //ft->SetParameters(0.985,5.0,0.3);
  ft->SetParameters(0.985,5.025,0.3);
  //ft->FixParameter(1,5.03); // semi-weighted average of BCD and EF
  //ft->FixParameter(2,0.395); // combined fit to BCD+EF / G+H 

  // ( 12.9*5.055+6.8*5.000)/(12.9+6.8)
  ft->FixParameter(1,5.036); // semi-weighted average of BCD/GH and EF/GH
  // ( 12.9*0.344 + 6.8*0.455)/(12.9+6.8)
  ft->FixParameter(2,0.391); // combined fit to BCD+EF / GH 

  // Log-sigmoid + powerlaw
  //TF1 *ft = new TF1("ft","[0]+(1-[0])/(1. + exp(-(log(x)-[1])/[2]))"
  //	       "*(1-[3]*pow(x,[4]))",30,2200);
  //ft->SetParameters(0.982,4.967,0.271,0.1,-0.2);
  // Double powerlaw
  //TF1 *ft = new TF1("ft","[4]-[0]*pow(x,[1])-[2]*pow(x,[3])",30,2200);
  //ft->SetParameters(0.05,-0.15,0.01,-0.3,1);
  

TF1 *fit(TTree *nt, TTree *ntMC, double ptmin,double ptmax)
{   
   //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),"",30,5.03,5.93);
//   TH1D *hBck = new TH1D(Form("hBck%d",count),"",40,5,6);
   
   TH1D *hMC = new TH1D(Form("hMC%d",count),"",30,5.03,5.93);
   // Fit function
   //QM2014
   //TString iNP="6.71675e+00*Gaus(x,5.30142e+00,8.42680e-02)/(sqrt(2*3.14159)*8.42680e-02)+4.06744e+01*Gaus(x,5.00954e+00,8.11305e-02)/(sqrt(2*3.14159)*8.11305e-02)+5.99974e-01*(2.376716*Gaus(x,5.640619,0.095530)/(sqrt(2*3.14159)*0.095530)+3.702342*Gaus(x,5.501706,0.046222)/(sqrt(2*3.14159)*0.046222))+1.31767e-01*(61.195688*Gaus(x,5.127566,0.087439)/(sqrt(2*3.14159)*0.087439)+58.943919*Gaus(x,5.246471,0.041983)/(sqrt(2*3.14159)*0.041983))";
   TString iNP="2.28629e1*Gaus(x,5.02606,6.84e-2)/(sqrt(2*3.14159)*(6.84e-2))+3.85695*Gaus(x,5.27701,0.04305)/(sqrt(2*3.14159)*(0.04305))";


   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+[6]*("+iNP+")");

   nt->Project(Form("h%d",count),"mass",Form("%s&&pt>%f&&pt<%f",seldata_2y_kpi.Data(),ptmin,ptmax));   
   ntMC->Project(Form("hMC%d",count),"mass",Form("%s&&pt>%f&&pt<%f",seldata_2y_kpi.Data(),ptmin,ptmax));   
//   nt->Project(Form("hBck%d",count),"mass",Form("%s&&pt>%f&&pt<%f&&(gen==23333||gen==41000)",seldata.Data(),ptmin,ptmax));   
//   nt2->Project(Form("hBck%d",count),"mass",Form("%s&&pt>%f&&pt<%f&&(gen==23333||gen==41000)",seldata.Data(),ptmin,ptmax));   

//   cout <<"nsig = "<<hBck->GetEntries();
   clean0(h);
   h->Draw();
   f->SetParLimits(4,-1000,0);
   f->SetParLimits(2,0.01,0.08);
   f->SetParLimits(8,0.01,0.1);
   f->SetParLimits(7,0,1);
   f->SetParLimits(6,0,1000);

   f->SetParameter(0,setparam0);
   f->SetParameter(1,setparam1);
   f->SetParameter(2,setparam2);
   f->SetParameter(8,setparam3);
   f->FixParameter(1,fixparam1);
   f->FixParameter(6,0);
   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));
   f->ReleaseParameter(6);
   
   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;

   cout<<"======= chi2 ======="<<endl;
   cout<<f->GetChisquare()<<endl;
   cout<<"===== chi2 end ====="<<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->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",count),"[0]*("+iNP+")");
   Bkpi->SetParameter(0,f->GetParameter(6));
   Bkpi->SetLineColor(kGreen+1);
   Bkpi->SetFillColor(kGreen+1);
   Bkpi->SetRange(5.00,6.00);
   Bkpi->SetLineStyle(1);
   Bkpi->SetFillStyle(3005);

   // function for signal shape plotting. take the fit result from f
   TF1 *mass = new TF1(Form("fmass",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);
//.........这里部分代码省略.........

void combinedFit() {

  TH1D * hB = new TH1D("hB","histo B",100,0,100);
  TH1D * hSB = new TH1D("hSB","histo S+B",100, 0,100);

  TF1 * fB = new TF1("fB","expo",0,100);
  fB->SetParameters(1,-0.05);
  hB->FillRandom("fB");

  TF1 * fS = new TF1("fS","gaus",0,100);
  fS->SetParameters(1,30,5);

  hSB->FillRandom("fB",2000);
  hSB->FillRandom("fS",1000);

  // perform now global fit

  TF1 * fSB = new TF1("fSB","expo + gaus(2)",0,100);

  ROOT::Math::WrappedMultiTF1 wfB(*fB,1);
  ROOT::Math::WrappedMultiTF1 wfSB(*fSB,1);

  ROOT::Fit::DataOptions opt;
  ROOT::Fit::DataRange rangeB;
  // set the data range
  rangeB.SetRange(10,90);
  ROOT::Fit::BinData dataB(opt,rangeB);
  ROOT::Fit::FillData(dataB, hB);

  ROOT::Fit::DataRange rangeSB;
  rangeSB.SetRange(10,50);
  ROOT::Fit::BinData dataSB(opt,rangeSB);
  ROOT::Fit::FillData(dataSB, hSB);

  ROOT::Fit::Chi2Function chi2_B(dataB, wfB);
  ROOT::Fit::Chi2Function chi2_SB(dataSB, wfSB);

  GlobalChi2 globalChi2(chi2_B, chi2_SB);

  ROOT::Fit::Fitter fitter;

  const int Npar = 6;
  double par0[Npar] = { 5,5,-0.1,100, 30,10};

  // create before the parameter settings in order to fix or set range on them
  fitter.Config().SetParamsSettings(6,par0);
  // fix 5-th parameter
  fitter.Config().ParSettings(4).Fix();
  // set limits on the third and 4-th parameter
  fitter.Config().ParSettings(2).SetLimits(-10,-1.E-4);
  fitter.Config().ParSettings(3).SetLimits(0,10000);
  fitter.Config().ParSettings(3).SetStepSize(5);

  fitter.Config().MinimizerOptions().SetPrintLevel(0);
  fitter.Config().SetMinimizer("Minuit2","Migrad");

  // fit FCN function directly
  // (specify optionally data size and flag to indicate that is a chi2 fit)
  fitter.FitFCN(6,globalChi2,0,dataB.Size()+dataSB.Size(),true);
  ROOT::Fit::FitResult result = fitter.Result();
  result.Print(std::cout);

  TCanvas * c1 = new TCanvas("Simfit","Simultaneous fit of two histograms",
                             10,10,700,700);
  c1->Divide(1,2);
  c1->cd(1);
  gStyle->SetOptFit(1111);

  fB->SetFitResult( result, iparB);
  fB->SetRange(rangeB().first, rangeB().second);
  fB->SetLineColor(kBlue);
  hB->GetListOfFunctions()->Add(fB);
  hB->Draw();

  c1->cd(2);
  fSB->SetFitResult( result, iparSB);
  fSB->SetRange(rangeSB().first, rangeSB().second);
  fSB->SetLineColor(kRed);
  hSB->GetListOfFunctions()->Add(fSB);
  hSB->Draw();


}

TF1* fit(TH1D* h, TH1D* hMCSignal, TH1D* hMCSwapped, Float_t ptmin, Float_t ptmax, Int_t j)
{
  TCanvas* c = new TCanvas(Form("c_%.0f_%.0f_%d",ptmin,ptmax,j),"",600,600);
  TF1* f = new TF1(Form("f_%.0f_%.0f_%d",ptmin,ptmax,j),"[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.5);
  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();
  
  hMCSignal->Fit(Form("f_%.0f_%.0f_%d",ptmin,ptmax,j),"q","",minhisto,maxhisto);
  hMCSignal->Fit(Form("f_%.0f_%.0f_%d",ptmin,ptmax,j),"q","",minhisto,maxhisto);
  f->ReleaseParameter(1);
  hMCSignal->Fit(Form("f_%.0f_%.0f_%d",ptmin,ptmax,j),"L q","",minhisto,maxhisto);
  hMCSignal->Fit(Form("f_%.0f_%.0f_%d",ptmin,ptmax,j),"L q","",minhisto,maxhisto);
  hMCSignal->Fit(Form("f_%.0f_%.0f_%d",ptmin,ptmax,j),"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_%d",ptmin,ptmax,j),"L q","",minhisto,maxhisto);
  hMCSwapped->Fit(Form("f_%.0f_%.0f_%d",ptmin,ptmax,j),"L q","",minhisto,maxhisto);
  hMCSwapped->Fit(Form("f_%.0f_%.0f_%d",ptmin,ptmax,j),"L q","",minhisto,maxhisto);
  hMCSwapped->Fit(Form("f_%.0f_%.0f_%d",ptmin,ptmax,j),"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);
  
  h->Fit(Form("f_%.0f_%.0f_%d",ptmin,ptmax,j),"q","",minhisto,maxhisto);
  h->Fit(Form("f_%.0f_%.0f_%d",ptmin,ptmax,j),"q","",minhisto,maxhisto);
  f->ReleaseParameter(1);
  f->SetParLimits(1,1.86,1.87);
  /*
  f->ReleaseParameter(11);
  f->SetParLimits(11,-1.,1.);
  */
  h->Fit(Form("f_%.0f_%.0f_%d",ptmin,ptmax,j),"L q","",minhisto,maxhisto);
  h->Fit(Form("f_%.0f_%.0f_%d",ptmin,ptmax,j),"L q","",minhisto,maxhisto);
  h->Fit(Form("f_%.0f_%.0f_%d",ptmin,ptmax,j),"L q","",minhisto,maxhisto);
  h->Fit(Form("f_%.0f_%.0f_%d",ptmin,ptmax,j),"L m","",minhisto,maxhisto);
  
  TF1* background = new TF1(Form("background_%.0f_%.0f_%d",ptmin,ptmax,j),"[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(Form("mass_%.0f_%.0f_%d",ptmin,ptmax,j),"[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));
  mass->SetFillColor(kOrange-3);
  mass->SetFillStyle(3002);
  mass->SetLineColor(kOrange-3);
  mass->SetLineWidth(3);
  mass->SetLineStyle(2);
  
  TF1* massSwap = new TF1(Form("massSwap_%.0f_%.0f_%d",ptmin,ptmax,j),"[0]*(1-[2])*Gaus(x,[1],[3]*(1+[4]))/(sqrt(2*3.14159)*[3]*(1+[4]))");
  massSwap->SetParameters(f->GetParameter(0),f->GetParameter(1),f->GetParameter(7),f->GetParameter(8),f->GetParameter(11));
  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);
//.........这里部分代码省略.........

TF1* fit(TString variable, TString variableplot, TTree* nt, TTree* ntMC, Double_t ptmin, Double_t ptmax, int isMC)
{
  static int count=0;
  count++;

  TCanvas* c= new TCanvas(Form("c%d",count),"",600,600);
  TH1D* h = new TH1D(Form("h-%d",count),"",nbinsmasshisto,minhisto,maxhisto);
  TH1D* hMCSignal = new TH1D(Form("hMCSignal-%d",count),"",nbinsmasshisto,minhisto,maxhisto);
  TH1D* hMCSwapped = new TH1D(Form("hMCSwapped-%d",count),"",nbinsmasshisto,minhisto,maxhisto);
  
  TF1* f = new TF1(Form("f%d",count),"[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+[6]*x*x*x", 1.7, 2.0);
  TString mycut=Form("(%s&&(%s)>%f&&(%s)<%f)",seldata.Data(),variable.Data(),ptmin,variable.Data(),ptmax);
  if(isMC==1) nt->Project(Form("h-%d",count),"Dmass",TCut(weight)*TCut(mycut));   
  else nt->Project(Form("h-%d",count),"Dmass",Form("(%s&&(%s)>%f&&(%s)<%f)",seldata.Data(),variable.Data(),ptmin,variable.Data(),ptmax));   
  
  ntMC->Project(Form("hMCSignal-%d",count),"Dmass",Form("%s*(%s&&(%s)>%f&&(%s)<%f&&(Dgen==23333))",weight.Data(),selmc.Data(),variable.Data(),ptmin,variable.Data(),ptmax));   
  ntMC->Project(Form("hMCSwapped-%d",count),"Dmass",Form("%s*(%s&&(%s)>%f&&(%s)<%f&&(Dgen==23344))",weight.Data(),selmc.Data(),variable.Data(),ptmin,variable.Data(),ptmax));   
  
  
  TFile *fout=new TFile(Form("FitsFiles/Fits_%s_%d.root",collisionsystem.Data(),count),"recreate");
  fout->cd();
  hMCSignal->Write();
  hMCSwapped->Write();
  h->Write();  
  fout->Close();
  
  
  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);
  h->GetEntries();
  
  hMCSignal->Fit(Form("f%d",count),"q","",minhisto,maxhisto);
  hMCSignal->Fit(Form("f%d",count),"q","",minhisto,maxhisto);
  f->ReleaseParameter(1);
  hMCSignal->Fit(Form("f%d",count),"L q","",minhisto,maxhisto);
  hMCSignal->Fit(Form("f%d",count),"L q","",minhisto,maxhisto);
  hMCSignal->Fit(Form("f%d",count),"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%d",count),"L q","",minhisto,maxhisto);
  hMCSwapped->Fit(Form("f%d",count),"L q","",minhisto,maxhisto);
  hMCSwapped->Fit(Form("f%d",count),"L q","",minhisto,maxhisto);
  hMCSwapped->Fit(Form("f%d",count),"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);
  
  h->Fit(Form("f%d",count),"q","",minhisto,maxhisto);
  h->Fit(Form("f%d",count),"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%d",count),"L q","",minhisto,maxhisto);
  h->Fit(Form("f%d",count),"L q","",minhisto,maxhisto);
  h->Fit(Form("f%d",count),"L q","",minhisto,maxhisto);
  h->Fit(Form("f%d",count),"L m","",minhisto,maxhisto);
  
  TF1* background = new TF1(Form("background%d",count),"[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(Form("fmass%d",count),"[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));
//.........这里部分代码省略.........