C++ TH1F::GetMaximum方法代码示例

void testLD(double scale=1.0){

  RooRealVar* zzmass = new RooRealVar("zzmass","m_{ZZ}",110,180);
  RooRealVar* z1mass = new RooRealVar("z1mass","m_{Z1}",0,180);
  RooRealVar* z2mass = new RooRealVar("z2mass","m_{Z2}",0,120);
  RooRealVar* costheta1 = new RooRealVar("costheta1","cos#theta_{1}",-1,1);  
  RooRealVar* costheta2 = new RooRealVar("costheta2","cos#theta_{2}",-1,1);
  RooRealVar* phi= new RooRealVar("phi","#Phi",-3.1415,3.1415);

  RooArgSet *argSet = new RooArgSet(*z1mass,*z2mass,*costheta1,*costheta2,*phi);

  vector<RooRealVar*> meas;
  meas.push_back(z1mass);
  meas.push_back(z2mass);
  meas.push_back(costheta1);
  meas.push_back(costheta2);
  meas.push_back(phi);

  cout << "Initializing PDFs" << endl;

  AngularPdfFactory SMHiggs(z1mass,z2mass,costheta1,costheta2,phi,zzmass);
  SMHiggs.makeSMHiggs();
  SMHiggs.makeParamsConst(true);

  AngularPdfFactory PSHiggs(z1mass,z2mass,costheta1,costheta2,phi,zzmass);
  PSHiggs.makePSHiggs();
  PSHiggs.makeParamsConst(true);

  LikelihoodToolbox<RooXZsZs_5D,RooXZsZs_5D> LD(SMHiggs.PDF,PSHiggs.PDF,NULL,meas);

  TFile* file1 = new TFile("../datafiles/SMHiggs_125_JHU.root");
  TTree* tree1 = (TTree*) file1->Get("angles");
  RooDataSet *data1 = new RooDataSet("data1","data1",tree1,RooArgSet(*zzmass,*z1mass,*z2mass,*costheta1,*costheta2,*phi));

  LD.LD_lineColor=4;
  LD.LD_bins=100;
  TH1F *SMLD =new TH1F("SMLD","SMLD",100,0,1);
  LD.makeLDdistribution(*SMLD,data1);

  TFile* file2 = new TFile("../datafiles/PSHiggs_125_JHU.root");
  TTree* tree2 = (TTree*) file2->Get("angles");
  RooDataSet *data2 = new RooDataSet("data2","data2",tree2,RooArgSet(*zzmass,*z1mass,*z2mass,*costheta1,*costheta2,*phi));

  LD.LD_lineColor=2;
  TH1F *PSLD = new TH1F("PSLD","PSLD",100,0,1);
  LD.makeLDdistribution(*PSLD,data2);

  if(SMLD->GetMaximum()>PSLD->GetMaximum()){
    SMLD->Draw();
    PSLD->Draw("SAME");
  }else{
    PSLD->Draw();
    SMLD->Draw("SAME");
  }
  
}

void drawHistos(TCanvas * C, TString filename, TString category, TTree* Tmine, TTree* Tother,TString var, int nbins, float xmin, float xmax, TString selection, TString myGroup, TString myRootFile, TString group, TString groupRootFile,TString mySel="1",TString groupSel="1"){
  TH1F* Hmine = new TH1F(TString("Hmine")+var,"",nbins,xmin,xmax); 
  Hmine->GetYaxis()->SetTitle(category);
  Hmine->GetXaxis()->SetTitle(var);
  Hmine->SetLineColor(1);
  Hmine->SetStats(0);
  TH1F* Hother = new TH1F(TString("Hother")+var,"",nbins,xmin,xmax); 
  Hother->GetYaxis()->SetTitle(category);
  Hother->GetXaxis()->SetTitle(var);
  Hother->SetLineColor(2);
  Hother->SetStats(0);

  TText TXmine;
  TXmine.SetTextColor(1);
  TXmine.SetTextSize(.04);
  TText TXother;
  TXother.SetTextColor(2);
  TXother.SetTextSize(.04);

  Tmine->Draw(var+">>"+Hmine->GetName(),selection+"*("+mySel+")");
  Tother->Draw(var+">>"+Hother->GetName(),selection+"*("+groupSel+")");

  ////Draw one histogram on top of the other
  C->Clear();
  //Hmine->Scale(1./Hmine->Integral());
  //Hother->Scale(1./Hother->Integral()); 
  //Hother->Scale(968134./688134.); //GGH e-tau
  if(Hmine->GetMaximum()>Hother->GetMaximum())
    Hmine->GetYaxis()->SetRangeUser(0,Hmine->GetMaximum()*1.1);
  else Hmine->GetYaxis()->SetRangeUser(0,Hother->GetMaximum()*1.1);
  Hmine->Draw("hist");
  Hother->Draw("histsame");

//   ///Draw the difference of the historgrams
//   TH1F*HDiff=(TH1F*)Hmine->Clone("HDiff");
//   HDiff->Add(Hother,-1);
//   int max= abs(HDiff->GetMaximum())>abs( HDiff->GetMinimum()) ?   abs(HDiff->GetMaximum()): abs( HDiff->GetMinimum());
//   HDiff->GetYaxis()->SetRangeUser(-2*(max>0?max:1),2*(max>0?max:1));
//   HDiff->Draw("hist");
//   TLine line;
//   line.DrawLine(HDiff->GetXaxis()->GetXmin(),0,HDiff->GetXaxis()->GetXmax(),0);

  //Print the integrals of the histograms a the top
  //TXmine.DrawTextNDC(.2,.965,myGroup+"_"+myRootFile+": "+(long)(Hmine->Integral(0,Hmine->GetNbinsX()+1)));
  //TXother.DrawTextNDC(.2,.93,group+"_"+groupRootFile+": "+(long)(Hother->Integral(0,Hother->GetNbinsX()+1)));
  TXmine.DrawTextNDC(.2,.965,myGroup+" : "+(long)(Hmine->Integral(0,Hmine->GetNbinsX()+1)));
  TXother.DrawTextNDC(.2,.93,group+": "+(long)(Hother->Integral(0,Hother->GetNbinsX()+1)));
  C->Print(filename);

  delete Hmine;
  delete Hother;
}

// get FWHHM
vector<double> getFWHM(RooRealVar *mass, RooAbsPdf *pdf, RooDataSet *data, double wmin=110., double wmax=130., double step=0.0004) {
 
  cout << "Computing FWHM...." << endl;
  double nbins = (wmax-wmin)/step;
  TH1F *h = new TH1F("h","h",int(floor(nbins+0.5)),wmin,wmax);
  if (data){
    pdf->fillHistogram(h,RooArgList(*mass),data->sumEntries());
  }
  else {
    pdf->fillHistogram(h,RooArgList(*mass));
  }
  
  double hm = h->GetMaximum()*0.5;
  double low = h->GetBinCenter(h->FindFirstBinAbove(hm));
  double high = h->GetBinCenter(h->FindLastBinAbove(hm));

  cout << "FWHM: [" << low << "-" << high << "] Max = " << hm << endl;
  vector<double> result;
  result.push_back(low);
  result.push_back(high);
  result.push_back(hm);
  result.push_back(h->GetBinWidth(1));
  
  delete h;
  return result;
}

void CompareGammaPoisson(float N=100)
{ 

    gStyle->SetOptStat(111111110);
    gStyle->SetStatW(0.3);                
    gStyle->SetStatH(0.25);                

    TRandom1 *frand = new TRandom1();
    TRandom3 rand(1234); // gamma

    int Nbins = (int)(N+TMath::Sqrt(N)*10);

    TH1F *hpoisson    = new TH1F("hpoisson",    "Poisson",   Nbins, -0.5,  Nbins-0.5);
    TH1F *hgamma      = new TH1F("hgamma",      "Gamma",     Nbins, -0.5,  Nbins-0.5);
    hpoisson->Sumw2();
    hgamma->Sumw2();

    for(int i=0; i<Npseudo; i++)  
    { 
        int checkpoint = (int)Npseudo/10;
        //if((i%checkpoint)==0) cout << "Generated " << i << "/" << Npseudo << " experiments " << endl;
        
        float Npoisson      = frand->Poisson(N);
        float Ngamma        = gsl_ran_gamma(2*N+1,0.5,rand);;

        hpoisson->Fill(Npoisson);
        hgamma->Fill(Ngamma);

    }
    hpoisson->Scale(1./hpoisson->Integral());
    hgamma->Scale(1./hgamma->Integral());

    h1cosmetic(hpoisson, Form("Poisson(black) and Gamma(red), N=%i",(int)N), kBlack, 3, 0, "");
    h1cosmetic(hgamma,   "Gamma", kRed, 3, 0, "");

    TCanvas *c = new TCanvas("c","c",400,400);
    c->cd(1);
    hpoisson->SetMaximum(hpoisson->GetMaximum()>hgamma->GetMaximum()?hpoisson->GetMaximum()*1.2:hgamma->GetMaximum()*1.2);
    hpoisson->Draw("HIST");
    hgamma->Draw("HIST SAME");
    c->Print(Form("fig/ComparePoissonGamma_N%i.pdf", (int)N));
}

void plotTrk() {
  gStyle->SetOptStat(0);
  gStyle->SetOptTitle(1);
  char c[50];
  c1->Clear();  
  c1->Divide(2,3); 
  for(int hid=0; hid<6; hid++){
    TVirtualPad* pad = c1->cd(hid+1); 
    pad->SetLogy(lTrkHist[hid]);
    double xmin, xmax, ymin=0.0, ymax=0.0;
    if(lTrkHist[hid]==1) ymin=0.1;
    for(int quad=0; quad<kFgtNumQuads; quad++){
      sprintf(c,"Quad%1s-%s",cquad[quad],cTrkHist[hid]);
      //printf("Getting %s\n",c);
      TH1F *h = histTrk[quad][hid] = (TH1F*)file->Get(c);
      xmin=h->GetXaxis()->GetXmin();
      xmax=h->GetXaxis()->GetXmax();
      double m=h->GetMaximum();
      if(ymax<m) ymax=m;
      //printf("quad=%d max=%6.1f ymax=%6.1f xmin=%6.1f xmax=%6.1f\n",quad,m,ymax,xmin,xmax);
    }
    ymax*=1.2; if(lTrkHist[hid]==1){ymax*=20.0;}
    sprintf(c,"%s",cTrkHist[hid]);
    TH2F *frame = new TH2F(c,c,1,xmin,xmax,1,ymin,ymax); frame->SetStats(0); frame->Draw();
    for(int quad=0; quad<kFgtNumQuads; quad++){
      TH1F *h=histTrk[quad][hid];
      h->SetLineColor(colorQuad[quad]); h->SetLineWidth(3); 
      if(hid<5){
	h->Draw("SAME");  
	float mean=h->GetMean();
	if(hid==1 || hid==4) {sprintf(c,"%s mean=%6.4f",cquad[quad],mean);}
	else                 {sprintf(c,"%s mean=%6.2f",cquad[quad],mean);}
      }else{
	h->SetMarkerColor(colorQuad[quad]); h->SetMarkerStyle(20); h->SetMarkerSize(1);
	h->Draw("SAME PL");
	sprintf(c,"Quad%s",cquad[quad]);
      }
      TText *t1;
      float x1= 0.2, x2= 0.55;
      float y1=0.8 - 0.07*quad;
      float y2=0.8 - 0.07*(quad-2);
      if(quad<2) { t1 = new TText(x1,y1,c); }
      else       { t1 = new TText(x2,y2,c); }
      t1->SetNDC();
      t1->SetTextSize(0.06); 
      t1->SetTextColor(colorQuad[quad]); 
      t1->Draw();
    }
  }
  c1->Update();
  save("trk");
}

void fithist() {
    //fit function ftotal to signal + background

    histgen();

    TFile *f = new TFile("background.root");
    background = (TH1F*)f->Get("background"); //pointer used in ftotal
    TH1F *result = (TH1F*)f->Get("result");

    TF1 *ftot = new TF1("ftot",ftotal,0,10,4);
    Double_t norm = result->GetMaximum();
    ftot->SetParameters(0.5*norm,5,.2,norm);
    ftot->SetParLimits(0,.3*norm,norm);

    result->Fit("ftot","b");
}

void plotMttResolution() {

  setTDRStyle();

  TCanvas *c = new TCanvas("c", "c", 800, 800);

  TH1F *mtt = (TH1F*) _file0->Get("mtt_resolution_mva");
  mtt->SetTitle("");
  mtt->SetName("m_{t#bar{t}}");
  mtt->GetXaxis()->SetTitle("m_{t#bar{t}} - m_{t#bar{t}}^{gen} (GeV)");
  mtt->SetLineColor(TColor::GetColor("#542437"));
  mtt->SetLineColor(TColor::GetColor("#8A9B0F"));
  mtt->SetLineColor(TColor::GetColor("#C02942"));
  mtt->SetFillColor(mtt->GetLineColor());
  mtt->SetFillStyle(3004);
  mtt->GetXaxis()->SetTitleOffset(1.2);

  TF1 *voigt = new TF1("voigt", "gaus", -100, 100);
	voigt->SetParName(0, "amp");
	voigt->SetParName(2, "mean");
	voigt->SetParName(1, "sigma");

  voigt->SetParameter(0, mtt->GetMaximum());
	voigt->SetParameter(2, 0);
	voigt->SetParameter(1, 100);

  mtt->Fit(voigt, "RVN");

  voigt->SetLineColor(TColor::GetColor("#8A9B0F"));
  voigt->SetLineWidth(2);

  mtt->Draw();
  voigt->Draw("same");

  TLatex* latex = new TLatex();
  latex->SetTextFont(42);
  latex->SetTextSize(0.033);
  latex->DrawLatexNDC(0.25, 0.84, TString::Format("mean = %.2f", voigt->GetParameter(1)));
  latex->DrawLatexNDC(0.25, 0.8, TString::Format("#sigma = %.2f", voigt->GetParameter(2)));

  gPad->Modified();
  gPad->Update();

  c->Print("mtt_resolution.pdf");
}

void test_TF1_FitErr2(){

  gStyle->SetOptFit(1111);
  gStyle->SetOptStat(0);

   histgen();
   TFile *f = new TFile("background.root");
   background = (TH1F*)f->Get("background"); //pointer used in ftotal
   TH1F *result = (TH1F*)f->Get("result");
   TF1 *ftot = new TF1("ftot",ftotal,0,10,5);
   Double_t norm = result->GetMaximum();
//   ftot->SetParameters(0.5*norm,5,.2,norm);
   ftot->SetParameters(0.5*norm,6,.5,5,-0.5);
   ftot->SetParLimits(0,.3*norm,norm);
//	 ftot->FixParameter(1,6);
//	 ftot->FixParameter(2,0.5);
	 ftot->SetParLimits(1,5.5,6.5);
	 ftot->SetParLimits(2,0.25,1);
//	 ftot->SetParLimits(4,-0.01,-50);
   result->Fit("ftot","b");
   result->Fit("ftot","b");
   result->Fit("ftot","b");
   result->Fit("ftot");
   result->Fit("ftot");
   result->Fit("ftot");

	 TFitResult fitResult=result->Fit("ftot");

	 double 


	 TF1 *bgf = TF1("bgf","pol1",0,10);
	 



}

void toyMC(int nsteps = 1e6) {
    
    Float_t LA = 9.2;
    Float_t LB = 10.3;
    
    Float_t SF = 492./LB;
    Float_t eSF = TMath::Sqrt(23.*23.+19.7*19.7)/LB;
//     Float_t OF = 358./LA;
//     Float_t eOF = 27./LA;
    Float_t OF = 358./LB;
    Float_t eOF = 31./LB;
    
    Float_t SigB = 188.+238-414;
    
    
    TH1F* hSig = new TH1F("hSig ; SF-R_{SF/OF}#timesOF ; ","Signal component",600,-100.,500.);
    
    hSig->SetLineColor(kRed+2);
    
    TRandom3* ran = new TRandom3();
    
    for ( int i=0; i<nsteps; ++i ) {
        
        Float_t nBSF = ran->Gaus(SF*LB,eSF*LB);
        Float_t nBOF = ran->Gaus(OF*LB,eOF*LB);
        Float_t rsfof = ran->Gaus(1.0,0.05);
        
        hSig->Fill(nBSF-nBOF*rsfof);
        
    }
    
    
    TCanvas* mycan = new TCanvas("mycan","",100,10,900,600);
    mycan->SetLogy();
    TH1F* hSigNorm = hSig->DrawNormalized("");
    hSigNorm->SetMinimum(1e-5);
    hSigNorm->Draw();
    
    // Find 95% CL range
    float integral = 0;
    int binStart = -1;
    while ( integral <= 0.025 ) {
        ++binStart;
        integral += hSigNorm->GetBinContent(binStart);
    }
    std::cout << integral << " up to " << hSigNorm->GetBinCenter(binStart) << std::endl;
    integral = 0;
    int binEnd = hSigNorm->GetNbinsX()+2;
    while ( integral <= 0.025 ) {
        --binEnd;
        integral += hSigNorm->GetBinContent(binEnd);
    }
    std::cout << integral << " up to " << hSigNorm->GetBinCenter(binEnd) << std::endl;
    
    // Draw 95% CL
    TBox* range95 = new TBox(hSigNorm->GetBinCenter(binStart),hSigNorm->GetMinimum(),hSigNorm->GetBinCenter(binEnd),1.2*hSigNorm->GetMaximum());
    range95->SetFillColor(kBlue-9);
    range95->SetFillStyle(1000);
    range95->SetLineColor(range95->GetFillColor());
    range95->Draw();
    
    hSigNorm->SetTitle("hSigNorm; \"S\" #equiv SF - R_{SF/OF}#timesOF ; pdf");
    hSigNorm->Draw("same");
    
    std::cout << "Integrating from 0 to " << SigB << ": " << std::endl;
    std::cout << hSigNorm->Integral(0,hSigNorm->FindBin(SigB)) <<std::endl;
    
    TLegend* legend = new TLegend(0.6,0.7,0.95,0.9,"","brNDC");
    legend->SetBorderSize(0);
    legend->AddEntry(hSigNorm,"Expected \"S\" for block B","l");
    legend->AddEntry(range95,"95% region","f");
    legend->Draw();

    mycan->RedrawAxis();
    mycan->SaveAs("toyMCexp.pdf");

    TArrow* a = new TArrow(SigB,hSigNorm->GetMaximum(),SigB,hSigNorm->GetMinimum()*1.1,0.02);
    a->SetLineColor(kBlue+2);
    a->Draw();


    TLegend* legend2 = new TLegend(0.6,0.6,0.95,0.7,"","brNDC");
    legend2->SetBorderSize(0);
    legend2->AddEntry(a,"Observed (p-value 0.6%)","l");
    legend2->Draw();
    

    mycan->SaveAs("toyMC.pdf");
    
}

void th2polyEurope(Int_t npoints=500000)
{
   Int_t i,j;
   Double_t lon1 = -25;
   Double_t lon2 =  35;
   Double_t lat1 =  34;
   Double_t lat2 =  72;
   Double_t R = (lat2-lat1)/(lon2-lon1);
   Int_t W    = 800;
   Int_t H    = (Int_t)(R*800);
   gStyle->SetStatX(0.28);
   gStyle->SetStatY(0.45);
   gStyle->SetStatW(0.15);

   // Canvas used to draw TH2Poly (the map)
   TCanvas *ce = new TCanvas("ce", "ce",0,0,W,H);
   ce->ToggleEventStatus();
   ce->SetGridx();
   ce->SetGridy();

   // Real surfaces taken from Wikipedia.
   const Int_t nx = 36;
   // see http://en.wikipedia.org/wiki/Area_and_population_of_European_countries
   const char *countries[nx] = {
      "france",     "spain",  "sweden",  "germany",       "finland",
      "norway",     "poland", "italy",   "yugoslavia",    "united_kingdom",
      "romania",    "belarus","greece",  "czechoslovakia","bulgaria",
      "iceland",    "hungary","portugal","austria",       "ireland",
      "lithuania",  "latvia", "estonia", "denmark",       "netherlands",
      "switzerland","moldova","belgium", "albania",       "cyprus",
      "luxembourg", "andorra","malta",   "liechtenstein", "san_marino", "monaco" };
   Float_t surfaces[nx] = {
      547030,        505580,   449964,      357021,        338145,
      324220,        312685,   301230,      255438,        244820,
      237500,        207600,   131940,      127711,        110910,
      103000,         93030,    89242,       83870,         70280,
      65200,         64589,    45226,       43094,         41526,
      41290,         33843,    30528,       28748,          9250,
      2586,           468,      316,         160,            61, 2};

   TH1F *h = new TH1F("h","Countries surfaces (in km^{2})",3,0,3);
   for (i=0; i<nx; i++) h->Fill(countries[i], surfaces[i]);
   h->LabelsDeflate();

   TFile::SetCacheFileDir(".");
   TFile *f;
   f = TFile::Open("http://root.cern.ch/files/europe.root","cacheread");

   if (!f) {
      printf("Cannot access europe.root. Is internet working ?\n");
      return;
   }

   TH2Poly *p = new TH2Poly(
             "Europe",
             "Europe (bin contents are normalized to the surfaces in km^{2})",
             lon1,lon2,lat1,lat2);
   p->GetXaxis()->SetNdivisions(520);
   p->GetXaxis()->SetTitle("longitude");
   p->GetYaxis()->SetTitle("latitude");

   p->SetContour(100);

   TMultiGraph *mg;
   TKey *key;
   TIter nextkey(gDirectory->GetListOfKeys());
   while ((key = (TKey*)nextkey())) {
      TObject *obj = key->ReadObj();
      if (obj->InheritsFrom("TMultiGraph")) {
         mg = (TMultiGraph*)obj;
         p->AddBin(mg);
      }
   }

   TRandom r;
   Double_t longitude, latitude;
   Double_t x, y, pi4 = TMath::Pi()/4, alpha = TMath::Pi()/360;

   gBenchmark->Start("Partitioning");
   p->ChangePartition(100, 100);
   gBenchmark->Show("Partitioning");

   // Fill TH2Poly according to a Mercator projection.
   gBenchmark->Start("Filling");
   for (i=0; i<npoints; i++) {
      longitude = r.Uniform(lon1,lon2);
      latitude  = r.Uniform(lat1,lat2);
      x         = longitude;
      y         = 38*TMath::Log(TMath::Tan(pi4+alpha*latitude));
      p->Fill(x,y);
   }
   gBenchmark->Show("Filling");

   Int_t nbins = p->GetNumberOfBins();
   Double_t maximum = p->GetMaximum();


   // h2 contains the surfaces computed from TH2Poly.
   TH1F *h2 = (TH1F *)h->Clone("h2");
   h2->Reset();
//.........这里部分代码省略.........

void pionContaminationVeta(){
  
  
	gStyle->SetOptStat(0);
	gStyle->SetOptFit(1);
	gStyle->SetPalette(1);
	
	char hist  [100];
	char hfile [100];
  
	
	sprintf(hfile,"~/Desktop/Research/2012IFF/rootFiles/nsigma20bintest35mil.root");  
  
	TFile* file = TFile::Open(hfile,"read");
	assert(file);
	
	
  double totalYieldWhole = 0;
  double pionYieldWhole = 0;
  double kpYieldWhole = 0;
  double electYieldWhile = 0;
	
	
	double pionFrac  [20];
	double kpFrac    [20];
	double electFrac [20];
	
	double pionCent  [20];
	double kpCent    [20];
	double electCent [20];
	
	double pionMax  [20];
	double kpMax    [20];
	double electMax [20];
	
	double pionSig  [20];
	double kpSig    [20];
	double electSig [20];
	
	
	double etastart[20];
	double etaend[20];
  
	double eta[20];
	
	
	double piYield[20];
	double highestYield = 0;
	
  double totalYeild[20];
	
	
	
  TH1F* hTotal = new TH1F("hTotal","hTotal",300,-10,10);

  
	for (int ihist=0; ihist<20; ihist++)
	{
		/*
     if (ihist != 0 || ihist != 1 || ihist != 2 || ihist != 3 ||)
     {
     cou
     }
     */
		sprintf(hist,"nSigmaPionBin_%i",ihist);
		cout << hist << endl;
		
		TH1F *hNsigma = (TH1F*) file->Get(hist);
	  assert(hNsigma);
    
    hTotal->Add(hNsigma, 1);
    cout << hTotal->GetEntries() << endl;
		
		TF1 *fitFunc = new TF1("fitFunc","gaus(0)+gaus(3)+gaus(6)",-10,10);
    TF1 *piFunc = new TF1("piFunc","gaus(0)",-10,10);
    TF1 *kpFunc = new TF1("kpFunc","gaus(0)",-10,10);
    TF1  *eFunc = new TF1( "eFunc","gaus(0)",-10,10);
    fitFunc->SetParName(0, "A_{#pi}");
    fitFunc->SetParName(1, "#LTn_{#pi}#GT");
    fitFunc->SetParName(2, "#sigma_{#pi}");
    fitFunc->SetParName(3, "A_{K/p}");
    fitFunc->SetParName(4, "#LTn_{K/p}#GT");
    fitFunc->SetParName(5, "#sigma_{K/p}");
    fitFunc->SetParName(6, "A_{e}");
    fitFunc->SetParName(7, "#LTn_{e}#GT");
    fitFunc->SetParName(8, "#sigma_{e}");
		
    fitFunc->SetParameter(0, hNsigma->GetMaximum());
    fitFunc->SetParameter(1, 0.);
    fitFunc->SetParameter(2, 1.);
		
    fitFunc->SetParameter(3, 0.05*hNsigma->GetMaximum());
    fitFunc->SetParameter(4,-3.);
    fitFunc->SetParameter(5, 1.5);
		
    fitFunc->SetParameter(6, 0.05*hNsigma->GetMaximum());
    fitFunc->SetParameter(7, 4.); //was at 5. but it seems like 4. works better
    fitFunc->SetParameter(8, 3.);
		
		hNsigma->Fit("fitFunc");
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  HistoMC_SUM->GetXaxis()->SetTitle(variableX.c_str());
  HistoMC_SUM->GetYaxis()->SetTitle(variableY.c_str());
  
  HistoRatioMC_SUM->Add(HistoRatioMC[iMC]);
  HistoRatioMC_SUM->GetXaxis()->SetTitle(Form("(%s)/(%s)",variableY.c_str(),variableX.c_str()));
  
  hsMC->Add(HistoMC[iMC]);
  if (iMC == 0) hsMC->Add(HistoMC[iMC]);
  leg->AddEntry(HistoMC[iMC],HistoMC[iMC]->GetTitle(),"f");
  legDown->AddEntry(HistoMC[iMC],HistoMC[iMC]->GetTitle(),"f");
 }
  
  
 ///==== make fit ====
 
 TH1D* htemp_prof = smartGausProfileXSQRTN(HistoDATA,2);  
 TGraphErrors grH = buildGEfromH_Personalized(*htemp_prof);
 grH.SetMarkerStyle(22);
 grH.SetMarkerSize(0.7);
 grH.SetMarkerColor(kBlue);
 grH.SetFillColor(kBlue);
 grH.SetFillStyle(3335);
 grH.SetLineWidth(1);
 grH.SetLineColor(kBlue);
 grH.GetXaxis()->SetTitle(variableX.c_str());
 grH.GetYaxis()->SetTitle(variableY.c_str());
 

 ///----------------------
 ///---- Plot results ----
 ///----------------------
 std::cerr << " MC expected : DATA = " << MC_Expected << " : " << HistoDATA->GetEntries() << std::endl;
 
 std::cerr << " HistoDATA->GetMaximum() = " << HistoDATA->GetMaximum() << std::endl;
 std::cerr << " HistoDATA->GetEntries() = " << HistoDATA->GetEntries() << std::endl;
 outFile->cd();
 
 TCanvas cResultDistro("cResultDistro","cResultDistro",900,900);
 cResultDistro.Divide(2,2);
 cResultDistro.cd(1);
 HistoDATA->DrawClone("colz");
 gPad->SetGrid();
 TString tLumiName = Form("#int L = %.4f pb^{-1}",luminosity);
 TLatex tLumi(1.5 * HistoDATA->GetMean(1),0.2 * HistoDATA->GetMaximum(),tLumiName.Data());
 tLumi.DrawClone();
 
 
 cResultDistro.cd(2);
 HistoRatioMC_SUM->Draw();
 HistoRatioDATA->Draw("E1same");
 TString tEleName = Form("%d electrons",(int) HistoDATA->GetEntries());
 TLatex tEle(2.0 * HistoRatioDATA->GetMean(1),0.3 * HistoRatioDATA->GetMaximum(),tEleName.Data());
 tEle.DrawClone();
 gPad->SetGrid();

 
 cResultDistro.cd(3);
 grH.Draw("APL");
 double minFit = gConfigParser -> readDoubleOption("Options::minFit");
 double maxFit = gConfigParser -> readDoubleOption("Options::maxFit");
 std::cout << ">>>>> Options::minFit      " << minFit  << std::endl;
 std::cout << ">>>>> Options::maxFit      " << maxFit  << std::endl;
 
 TF1* funz = new TF1 ("funz","pol1",minFit,maxFit);
//  funz->FixParameter(0,0);
 funz->SetParameter(0,0);

//.........这里部分代码省略.........
				for(int i=0; i<nPar; i++) par2[i] = par[i];
				fval2 = fval;
				}
				*/
		}

		/*
			printf("forcing new fit..\n");
			for(int i=0; i<nPar; i++) {
			printf("old,new = %10.5f %10.5f\n",parMin[i],par2[i]);
			fits[ind]->SetParameter(i,par2[i]);
			}
			*/

		// restore original fit
		fval = 0.0;
		gMinuit->Eval(nPar,0,fval,parMin,0);
		for(int i=0; i<nPar; i++) fits[ind]->SetParameter(i,parMin[i]);



		//extract fit error matrix
		gMinuit->mnemat(matrix.GetMatrixArray(),nPar);
		matrix.Print();

		for(int i=0; i<nPar; i++) {
			for(int j=0; j<nPar; j++) {
				printf("%10.5f",matrix(i,j)/sqrt(matrix(i,i)*matrix(j,j)));
			}
			printf("\n");
		}
		//matrix.Draw("text");

		float hm = h->GetMaximum();
		if(mode==0) {
			//TAxis* ax = h->GetXaxis();
			//ax->SetRangeUser(24.1,199.9);
			h->SetMaximum(1.2*hm);
			//h->SetMinimum(0.0);
		} else if(mode==1) {
			TAxis* ax = h->GetXaxis();
			ax->SetRangeUser(20.0,500.0);
			h->SetMaximum(1.15*hm);
			h->SetMinimum(0.0);
		}


		h->Draw();
		fits[ind]->SetLineColor(1);
		fits[ind]->SetLineWidth(2.0);
		fits[ind]->Draw("SAME");
		// find chi2's and KS's
		//AnaChiKs(h,fits[ind]);



		TAxis* ax,*ay;
		ax = h->GetXaxis(); 
		ay = h->GetYaxis();
		ax->SetTitle("m(#gamma#gamma) (GeV/c^{2})"); 
		ay->SetTitle("Entries/2 GeV/c^{2}");
		ax->CenterTitle(); ay->CenterTitle();
		ax->SetTitleOffset(0.9);
		ay->SetTitleOffset(1.0);
		ax->SetTitleSize(0.08);
		ay->SetTitleSize(0.07);

//.........这里部分代码省略.........
    float xdata, xdata1, xdata2; // combined isolation, pt, eta

    int flag = 1;
    while (flag!=-1){
      flag =fscanf(infile,"%f %f %f",&xdata, &xdata1, &xdata2);
      if( xdata1 >= ptbin && xdata1 < ptmax && xdata<20.) {
	if((strcmp(EBEE,"EB")==0 && TMath::Abs(xdata2)<1.45) ||
	   (strcmp(EBEE,"EE")==0 && TMath::Abs(xdata2)<2.5 && TMath::Abs(xdata2)>1.7) ) {
 	  dataColl.push_back(xdata);
	  hdata->Fill(xdata);
 	}
      } 
    }// keep reading files as long as text exists
    ndata = dataColl.size();

    printf("test print data 2 %2.3f \n", dataColl[2]);    
//     cout << "ndata in dataColl = " << ndata << endl;
    if ( ndata == 0 ) {
      printf(" no data to fit \n");
      return fitted;
    }
  }

  if(ndata==0) {
    printf(" ---  no events in the fit \n");
    return fitted;
  }
    
  //test fit the template and get PDFs
  TCanvas *c10 = new TCanvas("c10","c10",1000,500);
  c10->Divide(2,1);
  c10->cd(1);

  double par[20] = {hsig->GetMaximum(), 1., 0.6, 0.3,
		    hbkg->GetMaximum(), -.45, -0.05, 0.03, 1., 1., 1., 1.};

  if(strcmp(EBEE,"EE")==0) { par[2]=-0.1, par[3]=0.2; par[6]=-0.15; par[7]=0.02; };
  int fit_status;

  TF1 *f1 = new TF1("f1", exp_conv, -1., 20., 11);
  TF1 *fmcsigfit = new TF1("fmcsigfit", exp_conv, -1., 20., 11);
  fmcsigfit->SetLineColor(4);
  fmcsigfit->SetLineWidth(2);

  f1->SetNpx(10000);
  f1->SetParameters(par);
  f1->SetLineWidth(2);
  c10->cd(1);
  fit_status = hsig->Fit(f1,"","",-1., 5.);

  hsig->Draw();
  f1->Draw("same");
  if ( fit_status > 0 ) {
     printf("fit signal template failed. QUIT \n");
     return fitted;
  }
  if(para_index>0 && para_index<4){
    double tmppar = f1->GetParameter(para_index);
    f1->SetParameter(para_index, tmppar+para_sigma*f1->GetParError(para_index));
  }

  TF1 *fmcsig = (TF1*)f1->Clone();
  TF1 *fmcsigcorr = (TF1*)f1->Clone();
  fmcsig->SetNpx(10000);
  fmcsigcorr->SetNpx(10000);
  fmcsigfit->SetNpx(10000);

//------------------------------------------------------------------------------
// Draw projections and residuals
//------------------------------------------------------------------------------
void DrawWithRatio(TCanvas *canvas, char *cTitle,
                  TH1F *hNum, TH1F *hDen)
{
 
 // sanity check
 if (hNum->GetNbinsX() != hDen->GetNbinsX()){
   std::cout<< " *** Error: binning not consistent between data"
            << " and MC -> Exit!\n";
   return;
 }

 
 hNum->Sumw2();
 hDen->Sumw2();

 
 TH1F *hPull = (TH1F*)hNum ->Clone("hPull");
 hPull->Sumw2();
 hPull->Divide(hDen);
 
 //----------------------------------------------------------------------------
 // Create the pads
 //----------------------------------------------------------------------------
 TPad* pad1;
 TPad* pad2;
 
 pad1 = new TPad("pad1","This is pad1",0.02,0.30,0.98,0.98,0);
 pad2 = new TPad("pad2","This is pad2",0.02,0.01,0.98,0.29,0);
 
 pad1->SetLogx();
 pad2->SetLogx();
 pad1->SetBottomMargin(0.01);
 pad2->SetBottomMargin(0.33);
 pad2->SetTopMargin   (0.10);
 
 pad1->Draw(); // Projections pad
 pad2->Draw(); // Residuals   pad
        
        _leg2 = new TLegend(.68,.77,.98,.93);
        _leg2->AddEntry(hDen,"Gen ","l");
        _leg2->AddEntry(hNum,"Unfolded Reco","p");
 pad1->cd();
 hDen->Draw("histo");
 hNum->Draw("pe same");
 _leg2->Draw(); 
 PrintItLog(pad1,cTitle);

//   TLegend* leg = SetLegend(0.73, 0.7, 0.92, 0.89);
//   leg -> AddEntry(hDen," no mass cut","f");
//   leg -> AddEntry(hNum," 60 < M^{#mu #mu} < 120","p");
//   leg ->Draw("same");
 //----------------------------------------------------------------------------
 // Residuals pad
 //----------------------------------------------------------------------------
 pad2->cd();
 
 TAxis *xPull = NULL;
 TAxis *yPull = NULL;
 char xAxisName[200];
 sprintf(xAxisName,"%s",hDen->GetXaxis()->GetTitle());
 axis1F(hPull,xPull,yPull,xAxisName,"ratio");
 
 if (hPull->GetMaximum() > 100) {
   hPull->SetMinimum(-100);
   hPull->SetMaximum( 100);
 }
 
 hPull->GetXaxis()->SetLabelOffset(0.005);
 hPull->GetXaxis()->SetLabelSize  (0.11);
 hPull->GetXaxis()->CenterTitle(1);
 hPull->GetXaxis()->SetTitleOffset(1.10);
 hPull->GetXaxis()->SetTitleSize  (0.12);
 hPull->GetXaxis()->SetNdivisions(7);

 hPull->GetYaxis()->SetLabelSize  (0.09);
 hPull->GetYaxis()->CenterTitle(1);
 hPull->GetYaxis()->SetTitleOffset(0.5);
 hPull->GetYaxis()->SetTitleSize  (0.12);

 hPull->SetMaximum(1.5);
 hPull->SetMinimum(0.5);
 hPull->Draw("pe");

 pad2->Update();
 pad2->GetFrame()->DrawClone();

}