C++ TH2D::GetNbinsX方法代码示例

void Yield1S(){
  
  TFile *f = new TFile("Upsilon_2D_10ptbins.root");
  TH2D *H;
  H = (TH2D*)gFile->Get("fS1Yield");
  double tot(0.); double totE(0.); 
  
  ofstream OUT("Yield_1S.tex");
  OUT << "% ----------------------------------------------------------------------" << endl;
  OUT << "% -- Yields" << endl;
  
  for ( int x = 1; x <= H->GetNbinsX(); ++x ){
    OUT << Form("\\vdef{%iS_Rap_bin%iLowEdge} {\\ensuremath{ {%.1f } } }",1, x, H->GetXaxis()->GetBinLowEdge(x) ) << endl;
    if ( x == H->GetNbinsX() ) OUT << Form("\\vdef{%iS_Rap_bin%iHighEdge} {\\ensuremath{ {%.1f } } }",1 , x, H->GetXaxis()->GetBinUpEdge(x) ) << endl;
    //OUT << Form("\\vdef{%iS_bin%iContent} {\\ensuremath{ {%.0f } } }",1, x, H->GetBinContent(x) ) << endl;
    //OUT << Form("\\vdef{%iS_bin%iError} {\\ensuremath{ {%.0f } } }",1, x, H->GetBinError(x) ) << endl;
  }	
  
  for ( int x = 1; x <= H->GetNbinsY(); ++x ){
    OUT << Form("\\vdef{%iS_Pt_bin%iLowEdge} {\\ensuremath{ {%.1f } } }",1, x, H->GetYaxis()->GetBinLowEdge(x) ) << endl;
    if ( x == H->GetNbinsY() ) OUT << Form("\\vdef{%iS_Pt_bin%iHighEdge} {\\ensuremath{ {%.1f } } }",1 , x, H->GetYaxis()->GetBinUpEdge(x) ) << endl;
  }  
  
  for ( int x = 1; x <= H->GetNbinsX(); ++x ){
    for ( int y = 1; y <= H->GetNbinsY(); ++y ){
      OUT << Form("\\vdef{Yield%iS_bin%i_%iContent} {\\ensuremath{ {%.0f } } }",1, x, y, H->GetCellContent(x,y) ) << endl;
      OUT << Form("\\vdef{Yield%iS_bin%i_%iError} {\\ensuremath{ {%.0f } } }",1, x, y, H->GetCellError(x,y) ) << endl;
    }  
  }
  
  OUT.close();
  
}

TMatrixD convertHistToMatrix(const TH2D& from)
{
    TMatrixD matrix(from.GetNbinsX(),from.GetNbinsY());
    for (int row=0; row<from.GetNbinsX();++row){
        for (int col=0; col<from.GetNbinsY(); ++col) {
            (matrix)[row][col]=(from).GetBinContent(row+1,col+1);
        }
    }
    return matrix;
}

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

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

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

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

void Matching::postProcessing()
{
  for (unsigned int nsens = 0; nsens < _dutDevice->getNumSensors(); nsens++)
  {
    Mechanics::Sensor* sensor = _dutDevice->getSensor(nsens);

    TH2D* tot = _inPixelTot.at(nsens);
    TH2D* hits = _inPixelTracks.at(nsens);

    assert(hits->GetNbinsX() == tot->GetNbinsX() &&
           hits->GetNbinsY() == tot->GetNbinsY() &&
           "Matching: histograms should have the same number of bins");

    for (int nx = 1; nx <= tot->GetNbinsX(); nx++)
    {
      for (int ny = 1; ny <= tot->GetNbinsY(); ny++)
      {
        const unsigned int num = hits->GetBinContent(nx, ny);
        if (num == 0) continue;
        const double average = tot->GetBinContent(nx, ny) / (double)num;
        tot->SetBinContent(nx, ny, average);
      }
    }

    for (unsigned int axis = 0; axis < 2; axis++)
    {
      TH1D* dist = axis ? _matchDistX.at(nsens) :
                          _matchDistY.at(nsens);

      double pixelWidth = axis ? sensor->getPitchX() :
                                 sensor->getPitchY();
      double beamSigma = axis ? _refDevice->getSensor(0)->getPitchX() :
                                _refDevice->getSensor(0)->getPitchY();

      TF1* fit = Processors::fitPixelBeam(dist, pixelWidth, beamSigma, false);

      std::stringstream ss;
      ss << _dutDevice->getName() << sensor->getName() << "PixelBeamFit"
         << ( axis ? "X" : "Y") << _nameSuffix;
      _plotDir->WriteObject(fit, ss.str().c_str());
    }
  }
}

void mytestreweight(){


  TFile*finput=new TFile("fworking.root");
  TH1D*hBPtFONLL=(TH1D*)finput->Get("hBPt");
  TH1D*hDPtFONLL=(TH1D*)finput->Get("hDPt");
  TH2D*hD=(TH2D*)finput->Get("hD");
  
   hD->SetXTitle("B p_{T} (GeV/c)");
   hD->SetYTitle("D^{0} p_{T} (GeV/c)");
   hBPtFONLL->SetXTitle("B p_{T} (GeV/c)");
   hDPtFONLL->SetXTitle("D p_{T} (GeV/c)");
  
  TH2D *hDreweight = (TH2D*)hD->Clone("hDreweight");
  TH1D *hBPtPythia = hD->ProjectionX("hBPtPythia");
  TH1D *hDPtPythia = hD->ProjectionY("hDPtPythia");

   for (int x=1;x<=hDreweight->GetNbinsX()+1;x++){  //loop over the B pt bins
      if (hBPtPythia->GetBinContent(x)==0) continue;
      double ratio = hBPtFONLL->GetBinContent(x)/hBPtPythia->GetBinContent(x);  // in each pt bin we calculate the ratio of pythiaB/FONLLB 
       
      for (int y=1;y<=hDreweight->GetNbinsY()+1;y++){   //loop over the D pt bins
         double ratio2 = ratio;                                              //copy the ratio for each B bin
         double val = hDreweight->GetBinContent(x,y)*ratio2;
         double valError = hDreweight->GetBinError(x,y)*ratio2;
	 hDreweight->SetBinContent(x,y,val);
	 hDreweight->SetBinError(x,y,valError);
      }   
   }

  TH1D *hBPtPythiaReweight = hDreweight->ProjectionX("hBPtPythiaReweight");
  TH1D *hDPtPythiaReweight = hDreweight->ProjectionY("hDPtPythiaReweight");
  
  hBPtPythiaReweight->Divide(hBPtFONLL);
   
  TCanvas*canvas=new TCanvas("canvas","canvas",1000,500);  
  canvas->SetLogy();
  canvas->Divide(3,1);
  canvas->cd(1);
  hBPtFONLL->Draw();
  canvas->cd(2);
  hBPtPythia->Draw();
  canvas->cd(3);
  hBPtPythiaReweight->Draw();

}

void HitInfo::postProcessing()
{
  for (unsigned int nsens = 0; nsens < _device->getNumSensors(); nsens++)
  {
    TH2D* map = _totMap.at(nsens);
    TH2D* count = _totMapCnt.at(nsens);

    for (Int_t x = 1; x <= map->GetNbinsX(); x++)
    {
      for (Int_t y = 1; y <= map->GetNbinsY(); y++)
      {
        const double average = map->GetBinContent(x, y) / count->GetBinContent(x, y);
        map->SetBinContent(x, y, average);
      }
    }
  }
}

TH2D* frameTH2D(TH2D *in){

	Double_t xw = in->GetXaxis()->GetBinWidth(0);
	Double_t yw = in->GetYaxis()->GetBinWidth(0);

	Int_t nx = in->GetNbinsX();
	Int_t ny = in->GetNbinsY();

	Double_t x0 = in->GetXaxis()->GetXmin();
	Double_t x1 = in->GetXaxis()->GetXmax();

	Double_t y0 = in->GetYaxis()->GetXmin();
	Double_t y1 = in->GetYaxis()->GetXmax();

	TH2D *framed = new TH2D(
			Form("%s framed",in->GetName()),
			Form("%s framed",in->GetTitle()),
			nx + 2, x0-xw, x1+xw,
			ny + 2, y0-yw, y1+yw
			);

	//Copy over the contents
	for(int ix = 1; ix <= nx ; ix++){
		for(int iy = 1; iy <= ny ; iy++){
			framed->SetBinContent(1+ix, 1+iy, in->GetBinContent(ix,iy));
		}
	}
	//Frame with huge values
	nx = framed->GetNbinsX();
	ny = framed->GetNbinsY();
	for(int ix = 1; ix <= nx ; ix++){
		framed->SetBinContent(ix,  1, 1000.);
		framed->SetBinContent(ix, ny, 1000.);
	}
	for(int iy = 2; iy <= ny-1 ; iy++){
		framed->SetBinContent( 1, iy, 1000.);
		framed->SetBinContent(nx, iy, 1000.);
	}

	return framed;
}

void XS1S(){
  
  TFile *f = new TFile("Upsilon_2D_10ptbins.root");
  TH2D *H;
  H = (TH2D*)gFile->Get("fXS");
  double tot(0.); double totE(0.); 
  
  ofstream OUT("XS_1S.tex");
  OUT << "% ----------------------------------------------------------------------" << endl;
  OUT << "% -- XSections" << endl;
  
  for ( int x = 1; x <= H->GetNbinsX(); ++x ){
    for ( int y = 1; y <= H->GetNbinsY(); ++y ){
      OUT << Form("\\vdef{XS%iS_bin%i_%iContent} {\\ensuremath{ {%.4f } } }",1, x, y, H->GetCellContent(x,y) ) << endl;
      OUT << Form("\\vdef{XS%iS_bin%i_%iError} {\\ensuremath{ {%.4f } } }",1, x, y, H->GetCellError(x,y) ) << endl;
    }  
  }
  
  OUT.close();
  
}

std::pair<Matrix, Matrix> roothist_to_matrix(const TH2D & hist, bool transpose){
    int ngen = hist.GetNbinsY();
    int nreco = hist.GetNbinsX();
    if(transpose){
        swap(ngen, nreco);
    }
    
    Matrix m = Matrix(nreco, ngen);
    Matrix m_e = Matrix(nreco, ngen);
    
    for(int i=0; i<nreco; ++i){
        for(int j=0; j<ngen; ++j){
            int nx = i+1;
            int ny = j+1;
            if(transpose) swap(nx, ny);
            m(i,j) = hist.GetBinContent(nx, ny);
            m_e(i,j) = hist.GetBinError(nx, ny);
        }
    }
    return pair<Matrix, Matrix>(move(m), move(m_e));
}

void drawJetFragBalance_DRDiff(
			       TString inFileName=    "plot/jfhCorrEtaPtBin4RBin20v2_HCPR_J50U_djcalo_Cent0to30_Aj0to100_SubEtaRefl.root",
			       TString inFileNameHyPy="plot/jfhCorrEtaPtBin4RBin20v2_Hydjet_djcalo_Cent0to30_Aj0to100_SubEtaRefl.root",
    TString title = "test",
    Int_t drawMode=1,
			       Int_t doLeg=1,
bool cumulative = 1
    ) {
  TFile *f = new TFile(inFileName);
  TString inFileNameStrip(inFileName); inFileNameStrip.ReplaceAll(".root","");

  TH2D * hPtPNDR = (TH2D*) f->Get("hPtPNDR");
  TH2D * hPtPADR = (TH2D*) f->Get("hPtPADR");
  TH2D * hPtPNDRBg = (TH2D*) f->Get("hPtPNDRBg");
  TH2D * hPtPADRBg = (TH2D*) f->Get("hPtPADRBg");
  TH2D * hPtPNDRBgSub = (TH2D*)hPtPNDR->Clone(inFileNameStrip+"hPtPNDRBgSub");
  TH2D * hPtPADRBgSub = (TH2D*)hPtPADR->Clone(inFileNameStrip+"hPtPADRBgSub");
  hPtPNDRBgSub->Add(hPtPNDR,hPtPNDRBg,1,-1);
  hPtPADRBgSub->Add(hPtPADR,hPtPADRBg,1,-1);

  TFile *fhypy = new TFile(inFileNameHyPy);
  TH2D * hPtPNDRHyPy = (TH2D*) fhypy->Get("hPtPNDR");
  TH2D * hPtPADRHyPy = (TH2D*) fhypy->Get("hPtPADR");
  TH2D * hPtPNDRBgHyPy = (TH2D*) fhypy->Get("hPtPNDRBg");
  TH2D * hPtPADRBgHyPy = (TH2D*) fhypy->Get("hPtPADRBg");
  TH2D * hPtPNDRBgSubHyPy = (TH2D*)hPtPNDR->Clone(inFileNameStrip+"hPtPNDRBgSubHyPy");
  TH2D * hPtPADRBgSubHyPy = (TH2D*)hPtPADR->Clone(inFileNameStrip+"hPtPADRBgSubHyPy");
  hPtPNDRBgSubHyPy->Add(hPtPNDRHyPy,hPtPNDRBgHyPy,1,-1);
  hPtPADRBgSubHyPy->Add(hPtPADRHyPy,hPtPADRBgHyPy,1,-1);

  // Get Pt info
  Int_t numBinsPt=hPtPNDR->GetNbinsX();
  Int_t numBinsDR=hPtPNDR->GetNbinsY();
  TH1D * hPt = (TH1D*)hPtPNDR->ProjectionX("hPt");
  cout << "Pt bins: " << numBinsPt << endl;

  Double_t totPtBgSubNr=hPtPNDRBgSub->Integral();
  Double_t totPtBgSubAw=hPtPADRBgSub->Integral();

  Double_t totPtBgSubNrHyPy=hPtPNDRBgSubHyPy->Integral();
  Double_t totPtBgSubAwHyPy=hPtPADRBgSubHyPy->Integral();

  int ptUp = 2;
  TH1D * hDRBgSubNr = (TH1D*)hPtPNDRBgSub->ProjectionY(inFileNameStrip+"hDRBgSubNr",1,ptUp);
  TH1D * hDRBgSubAw = (TH1D*)hPtPADRBgSub->ProjectionY(inFileNameStrip+"hDRBgSubAw",1,ptUp);
  TH1D * hDRBgSubNrHyPy = (TH1D*)hPtPNDRBgSubHyPy->ProjectionY(inFileNameStrip+"hDRBgSubNrHyPy",1,ptUp);
  TH1D * hDRBgSubAwHyPy = (TH1D*)hPtPADRBgSubHyPy->ProjectionY(inFileNameStrip+"hDRBgSubAwHyPy",1,ptUp);
  // Print
  cout << Form("%.1f < p_{T} < %.1f GeV/c: ",hPt->GetBinLowEdge(1),hPt->GetBinLowEdge(ptUp+1)) << " SigSubBkg Integral - Nr: " << endl;
  cout << " Data - Nr: " << hDRBgSubNr->Integral() << " Aw: " << hDRBgSubAw->Integral() << endl;
  cout << " Pythia+Hydjet - Nr: " << hDRBgSubNrHyPy->Integral() << " Aw: " << hDRBgSubAwHyPy->Integral() << endl;

  if (drawMode==1) TCanvas * c6 = new TCanvas("c6","c6",500,500);
  hDRBgSubNr->Scale(1./totPtBgSubNr);
  hDRBgSubAw->Scale(1./totPtBgSubAw);
  hDRBgSubNrHyPy->Scale(1./totPtBgSubNrHyPy);
  hDRBgSubAwHyPy->Scale(1./totPtBgSubAwHyPy);

  if(cumulative){
     hDRBgSubNr = IntegrateFromLeft(hDRBgSubNr);
     hDRBgSubAw = IntegrateFromLeft(hDRBgSubAw);
     hDRBgSubNrHyPy = IntegrateFromLeft(hDRBgSubNrHyPy);
     hDRBgSubAwHyPy = IntegrateFromLeft(hDRBgSubAwHyPy);
  }

  // Set Styles
  hDRBgSubNr->SetMarkerStyle(kOpenSquare);
  mcStyle1(hDRBgSubNrHyPy);
  mcStyle2(hDRBgSubAwHyPy);
  hDRBgSubNrHyPy->SetMarkerStyle(0);
  hDRBgSubAwHyPy->SetMarkerStyle(0);
  // Draw
  hDRBgSubNrHyPy->SetTitle(";#DeltaR_{max};F(#DeltaR<#DeltaR_{max})");
  hDRBgSubNrHyPy->SetAxisRange(0,0.784,"X");
  hDRBgSubNrHyPy->SetAxisRange(0,0.7,"Y");
  if(!cumulative){
     hDRBgSubNrHyPy->SetAxisRange(0,0.1,"Y");
  }
  fixedFontHist(hDRBgSubNrHyPy);
  hDRBgSubNrHyPy->DrawCopy("Ehist");
  hDRBgSubAwHyPy->DrawCopy("Ehistsame");
  hDRBgSubNr->DrawCopy("Esame");
  hDRBgSubAw->DrawCopy("Esame");

  if (doLeg==1) {
     TLegend *leg = new TLegend(0.302407,0.67,0.7536548,0.9324599);
    leg->SetFillStyle(0);
    leg->SetBorderSize(0);
    leg->SetTextFont(63);
    leg->SetTextSize(16);
    leg->AddEntry(hDRBgSubNr,Form("%.1f < p_{T} < %.1f GeV/c",hPt->GetBinLowEdge(1),hPt->GetBinLowEdge(ptUp+1)),"");
    leg->AddEntry(hDRBgSubNr,"Data Leading Jet","pl");
    leg->AddEntry(hDRBgSubAw,"Data SubLeading Jet","pl");
    leg->AddEntry(hDRBgSubNrHyPy,"MC Leading Jet","l");
    leg->AddEntry(hDRBgSubAwHyPy,"MC SubLeading Jet","l");
    leg->Draw();
  }
}

//.........这里部分代码省略.........
  ifs.open(conf.Data());
  assert(ifs.is_open());
  string line;
  while(getline(ifs,line)) {
    if(line[0]=='#') continue;
    
    string fname;
    Int_t color, linesty;
    stringstream ss(line);
    ss >> fname >> color >> linesty;
    string label = line.substr(line.find('@')+1);
    fnamev.push_back(fname);
    labelv.push_back(label);
    colorv.push_back(color);
    linev.push_back(linesty);
  }
  ifs.close();

  // Create output directory
  gSystem->mkdir(outputDir,kTRUE);
  
  //
  // Get efficiency
  //
  TFile *dataHLTEffFile = new TFile(dataHLTEffName);
  CEffUser2D dataHLTEff;
  TH2D *hHLTErr=0, *hHLTErrB=0, *hHLTErrE=0;
  if(dataHLTEffName) {
    dataHLTEff.loadEff((TH2D*)dataHLTEffFile->Get("hEffEtaPt"),
                       (TH2D*)dataHLTEffFile->Get("hErrlEtaPt"),
                       (TH2D*)dataHLTEffFile->Get("hErrhEtaPt"));
    
    TH2D* h =(TH2D*)dataHLTEffFile->Get("hEffEtaPt");
    hHLTErr  = new TH2D("hHLTErr", "",h->GetNbinsX(),h->GetXaxis()->GetXmin(),h->GetXaxis()->GetXmax(),
                                      h->GetNbinsY(),h->GetYaxis()->GetXmin(),h->GetYaxis()->GetXmax());
    hHLTErrB = new TH2D("hHLTErrB","",h->GetNbinsX(),h->GetXaxis()->GetXmin(),h->GetXaxis()->GetXmax(),
                                      h->GetNbinsY(),h->GetYaxis()->GetXmin(),h->GetYaxis()->GetXmax());
    hHLTErrE = new TH2D("hHLTErrE","",h->GetNbinsX(),h->GetXaxis()->GetXmin(),h->GetXaxis()->GetXmax(),
                                      h->GetNbinsY(),h->GetYaxis()->GetXmin(),h->GetYaxis()->GetXmax());
  }
  
  TFile *zeeHLTEffFile = new TFile(zeeHLTEffName);
  CEffUser2D zeeHLTEff;
  if(zeeHLTEffName) {
    zeeHLTEff.loadEff((TH2D*)zeeHLTEffFile->Get("hEffEtaPt"),
                      (TH2D*)zeeHLTEffFile->Get("hErrlEtaPt"),
                      (TH2D*)zeeHLTEffFile->Get("hErrhEtaPt"));
  }
  
  TFile *dataGsfSelEffFile = new TFile(dataGsfSelEffName);
  CEffUser2D dataGsfSelEff;
  TH2D *hGsfSelErr=0, *hGsfSelErrB=0, *hGsfSelErrE=0;
  if(dataGsfSelEffName) {
    dataGsfSelEff.loadEff((TH2D*)dataGsfSelEffFile->Get("hEffEtaPt"),
                       (TH2D*)dataGsfSelEffFile->Get("hErrlEtaPt"),
                       (TH2D*)dataGsfSelEffFile->Get("hErrhEtaPt"));
    
    TH2D* h =(TH2D*)dataGsfSelEffFile->Get("hEffEtaPt");
    hGsfSelErr  = new TH2D("hGsfSelErr", "",h->GetNbinsX(),h->GetXaxis()->GetXmin(),h->GetXaxis()->GetXmax(),
                                      h->GetNbinsY(),h->GetYaxis()->GetXmin(),h->GetYaxis()->GetXmax());
    hGsfSelErrB = new TH2D("hGsfSelErrB","",h->GetNbinsX(),h->GetXaxis()->GetXmin(),h->GetXaxis()->GetXmax(),
                                      h->GetNbinsY(),h->GetYaxis()->GetXmin(),h->GetYaxis()->GetXmax());
    hGsfSelErrE = new TH2D("hGsfSelErrE","",h->GetNbinsX(),h->GetXaxis()->GetXmin(),h->GetXaxis()->GetXmax(),
                                      h->GetNbinsY(),h->GetYaxis()->GetXmin(),h->GetYaxis()->GetXmax());
  }
  

//.........这里部分代码省略.........
        myPadSetUp(myPad4,0.15,0.04,0.04,0.15);
        myPad4->Draw();
        myPad4->SetLogy();

        myPad4->cd();

        TH1D* matpp = (TH1D*)matp->ProjectionX("zxc3",0,100);
        matpp->Draw("");

        cout << "material in cut- " << matpp->Integral(191,211) << endl;
        cout << "material all- " << matpp->Integral(1,400) << endl;

        postprocess(myCan4,Form("DCAxy%sProX",status),rWrite,rPerformance,system);

        cout  << "in cut: " << endl;
        cout << "prim - " << primpp->Integral(191,211) / (primpp->Integral(191,211)+weakpp->Integral(191,211)+matpp->Integral(191,211)) << endl;
        cout << "weak - " << weakpp->Integral(191,211) / (primpp->Integral(191,211)+weakpp->Integral(191,211)+matpp->Integral(191,211)) << endl;
        cout << "mat - " << matpp->Integral(191,211) / (primpp->Integral(191,211)+weakpp->Integral(191,211)+matpp->Integral(191,211)) << endl;

        cout << endl << "in cut / all " << endl;
        cout << "prim - " << primpp->Integral(191,211) / primpp->Integral(1,400) << endl;
        cout << "weak - " << weakpp->Integral(191,211) / weakpp->Integral(1,400) << endl;
        cout << "mat - " << matpp->Integral(191,211) / matpp->Integral(1,400) << endl;


        // _____sum____

        TH2D* psum = new TH2D("psum","",400, -2.0, 2.0, 100,0.0,2.0);
        psum->GetXaxis()->SetTitle("DCA_{xy} (cm)");
        psum->GetYaxis()->SetTitle("Number of Entries (normalized)");
        psum->GetXaxis()->SetLimits(-2,2);


        for (int i = 0; i < primp->GetNbinsX(); i++) {
            for (int j = 0; j < primp->GetNbinsY(); j++) {
                psum->SetBinContent(i,j,primp->GetBinContent(i,j)+weakp->GetBinContent(i,j)+matp->GetBinContent(i,j));
            }
        }

        gStyle->SetOptStat(0);

        TCanvas *cansum = new TCanvas("cansum",cStamp1,600,400);
        cansum->Draw();
        cansum->cd();

        TPad *padsum = new TPad("padsum", "The pad4",0,0,1,1);
        myPadSetUp(padsum,0.15,0.04,0.04,0.15);
        padsum->Draw();
        padsum->SetLogy();

        padsum->cd();

        // TCanvas* cansum = new TCanvas("cansum","cansum");
        // cansum->SetLogy();
        //psum->Draw("colz");
        TH1D* asd0 = (TH1D*)psum->ProjectionX("zxc",0,100);

        // asd0->Scale(1./asd0->Integral());
        // TFile* fout = new TFile("dca.root","update");
        // asd0->SetName(Form("dcaxyMC%d",isMC));
        // asd0->Write();

        asd0->GetXaxis()->SetTitle("DCA_{xy} (cm)");
        asd0->GetYaxis()->SetTitle("Number of Entries (normalized)");
        //asd0->SetMaximum(5000);
        //asd0->SetMinimum(0.00008);

TCanvas* pHitSpecPosGen( )
{
  TCanvas* c = new TCanvas("cHitSpecPosGen","cHitSpecPosGen",1200,600);
  c->Divide(2,1);
  TVirtualPad* p; TH2D *h; TH2D *hAtPhi0;
  TText t; t.SetTextColor(4);

  p =c->cd(1); p->SetLogy(); p->SetGrid(0,1);
  h = (TH2D*)gROOT->FindObject("hHitSpec_PtVsPhiGen_Minus"); 
  hAtPhi0 = (TH2D*)gROOT->FindObject("hHitSpec_PtVsPhiGenAtPhi0_Minus"); 
  h->GetYaxis()->SetRange(4,33);
  h->GetXaxis()->SetNdivisions(505); h->GetXaxis()->SetLabelSize(0.04);
  h->DrawCopy("box");
  hAtPhi0->SetLineColor(2);
  hAtPhi0->DrawCopy("box same");
  t.DrawTextNDC(0.17,0.15, "BARREL MU MINUS");
  std::cout <<h->GetTitle() << std::endl;
  for (int iy = 1; iy <=h->GetNbinsY(); iy++) {
  std::cout <<" pt: " << h->GetYaxis()->GetBinLowEdge(iy);
  double xmin=100.;
  double xAtMax = 0.; double valAtMax = 0.;
  double xmax=-100.;
  for (int ix = 1; ix <=h->GetNbinsX(); ix++) {
     double val = h->GetBinContent(ix,iy); 
     double xbmin = h->GetXaxis()->GetBinLowEdge(ix);
     double xcent = h->GetXaxis()->GetBinCenter(ix);
     double xbmax = h->GetXaxis()->GetBinUpEdge(ix);
     if (val > valAtMax) { valAtMax = val; xAtMax = xcent; }
     if (val > 1 && (xbmin < xmin)) xmin = xbmin;    
     if (val > 1 && (xbmax > xmax)) xmax = xbmax;    
  }
  std::cout <<"  set DPHI0 = "<<1.025-xAtMax 
            <<"; set DPHI_MARGIN = " <<std::max( (xAtMax-xmin), (xmax-xAtMax))
<<";" << std::endl;
  }

  p =c->cd(2); p->SetLogy(); p->SetGrid(1,1);
  h = (TH2D*)gROOT->FindObject("hHitSpec_PtVsPhiGen_Plus"); 
  hAtPhi0 = (TH2D*)gROOT->FindObject("hHitSpec_PtVsPhiGenAtPhi0_Plus"); 
  h->GetYaxis()->SetRange(4,32);
  h->GetXaxis()->SetNdivisions(505); h->GetXaxis()->SetLabelSize(0.04);
  h->DrawCopy("box");
  hAtPhi0->SetLineColor(2);
  hAtPhi0->DrawCopy("box same");
  t.DrawTextNDC(0.17,0.15, "BARREL MU PLUS");
  std::cout <<h->GetTitle() << std::endl;
  for (int iy = 1; iy <=h->GetNbinsY(); iy++) {
  std::cout <<" pt: " << h->GetYaxis()->GetBinLowEdge(iy);
  double xmin=100.;
  double xAtMax = 0.; double valAtMax = 0.;
  double xmax=-100.;
  for (int ix = 1; ix <=h->GetNbinsX(); ix++) {
     double val = h->GetBinContent(ix,iy); 
     double xbmin = h->GetXaxis()->GetBinLowEdge(ix);
     double xcent = h->GetXaxis()->GetBinCenter(ix);
     double xbmax = h->GetXaxis()->GetBinUpEdge(ix);
     if (val > valAtMax) { valAtMax = val; xAtMax = xcent; }
     if (val > 1 && (xbmin < xmin)) xmin = xbmin;    
     if (val > 1 && (xbmax > xmax)) xmax = xbmax;    
  }
  std::cout <<"  set DPHI0 = "<<1.025-xAtMax 
            <<"; set DPHI_MARGIN = " <<std::max( (xAtMax-xmin), (xmax-xAtMax))
<<";" << std::endl;
  }

  return c;
}

void computeCleoPdf()
{
  //TString fName = "CLEO_K3PiScan_2009";
  //TString hName = "deltachisq";

  //TString fName = "CLEO_K3PiScan_2014";
  //TString hName = "deltachisq";

  //TString fName = "CLEO_KSKpiScan_2012_FullDP";
  //TString hName = "loglik";

  //TString fName = "CLEO_KSKpiScan_2012_Kst";
  //TString hName = "loglik";

  //TString fName = "CLEO_KPiPi0Scan_2014";
  //TString hName = "deltachisq";

  TString fName = "CLEO_LHCb_K3PiScan_2015";
  TString hName = "deltachisq";

  //
  // START THE CONVERSION
  //
  TFile *fExpNll = new TFile(fName+".root", "ro");
  TH2D *hExpNll = (TH2D*)fExpNll->Get(hName);
  assert(hExpNll);

  int nX = hExpNll->GetNbinsX();
  int nY = hExpNll->GetNbinsY();

  // initialize a new histogram
  TH2D *hNew = new TH2D("hPdf", "hPdf",
	nX, hExpNll->GetXaxis()->GetXmin(), hExpNll->GetXaxis()->GetXmax(),
	nY, hExpNll->GetYaxis()->GetXmin()/180.*TMath::Pi(), hExpNll->GetYaxis()->GetXmax()/180.*TMath::Pi());
  for ( int ii=1; ii<nX; ii++ )
  for ( int jj=1; jj<nY; jj++ )
  {
    hNew->SetBinContent(ii, jj, 0);
  }

  TCanvas *c1 = new TCanvas("c1", "c1", 800, 400);
  c1->Divide(2,1);
  c1->cd(1);
  hExpNll->Draw("colz");

  // time the macro
  TStopwatch t;
  t.Start();

  float max = 0;
  float xMax = 0;
  float yMax = 0;
  float minEval = 1e6;

  for ( int i=1; i<=nX; i++ )
  for ( int j=1; j<=nY; j++ )
  {
    // exponentionate to turn chi2 into likelihood
    float eVal = hExpNll->GetBinContent(i,j);
    // hack for bad histogram
    if (eVal>9) eVal=9;
    if (eVal>0){
	    eVal = TMath::Exp(-eVal/2.);
	    hNew->SetBinContent(i, j, eVal);
    }
    // find coordinates of maximum in the physical region, R>0
    if ( max < eVal && hNew->GetXaxis()->GetBinCenter(i)>0 )
    {
	    xMax = hNew->GetXaxis()->GetBinCenter(i);
	    yMax = hNew->GetYaxis()->GetBinCenter(j);
	    max = eVal;
    }
    // find minimum likeilihood value in the tails - we'll set any zero bin to this minimum, usually around 1e-25
    if ( 0.0 < eVal && eVal < minEval ) minEval = eVal;
  }

  // set zero bins to minimum likelihood value found
  for ( int i=1; i<=nX; i++ )
  for ( int j=1; j<=nY; j++ )
  {
    if ( hNew->GetBinContent(i,j)<=0.0 ) hNew->SetBinContent(i, j, minEval);
  }

  hNew->Scale(1./hNew->Integral());

  t.Stop();
  t.Print();

  cout << "xMax = " << xMax << " yMax = " << yMax << endl;

  c1->cd(2);
  hNew->Draw("colz");

  cout << "saving histogram " << fName+"_pdf.root" << endl;
  TFile *fNew = new TFile(fName+"_pdf.root", "recreate");
  hNew->Write();
  fNew->Close();
}

void run_HMx(){
  gROOT->LoadMacro("Fitter.cpp");

  // mass range
  double xmin = 50., xmax = 130.;

  // data
  TFile *data=new TFile("result_data_50_130.root");
  gDirectory->cd("Eff_Hist");

  const Int_t NBins = 31;
  const Int_t NBins_1D = 60;

  // get histogram before cut
  TH2D *before = (TH2D *)gROOT->FindObject("ZMass_vs_DetEta_Before_HMx_CC");
  TH2D *after = (TH2D *)gROOT->FindObject("ZMass_vs_DetEta_After_HMx_CC");

  Int_t NBins_X = before->GetNbinsX();
  Double_t Axis_XMin = before->GetXaxis()->GetXmin();
  Double_t Axis_XMax = before->GetXaxis()->GetXmax();

  if(NBins_X != NBins) {
    cout<<"NBins="<<NBins<<" while NBins_X="<<NBins_X<<" make sure they are the same"<<endl;
    exit;
  }

  Int_t NBins_Y = before->GetNbinsY();
  Double_t Axis_YMin = before->GetYaxis()->GetXmin();
  Double_t Axis_YMax = before->GetYaxis()->GetXmax();

  cout<<"HMx bin "<<NBins_X<<" from "<<Axis_XMin<<" to "<<Axis_XMax<<endl;
  cout<<"Mass bin "<<NBins_Y<<" from "<<Axis_YMin<<" to "<<Axis_YMax<<endl;

  // 1D histogram
  TH1D *before_mass[NBins], *after_mass[NBins], *diff_mass[NBins];
  double variable[NBins], variable_err[NBins];

  char name[50];
  for (int xbin = 0; xbin < NBins_X; xbin ++){
    sprintf(name, "%s%d", "before_", xbin);
    before_mass[xbin] = new TH1D(name, name, NBins_Y, Axis_YMin, Axis_YMax);
    sprintf(name, "%s%d", "after_", xbin);
    after_mass[xbin] = new TH1D(name, name, NBins_Y, Axis_YMin, Axis_YMax);
    sprintf(name, "%s%d", "diff_", xbin);
    diff_mass[xbin] = new TH1D(name, name, NBins_Y, Axis_YMin, Axis_YMax);

    variable[xbin] = before->GetXaxis()->GetBinCenter(xbin);
    variable_err[xbin] = 0.;

    for(int ybin=0; ybin < NBins_Y; ybin ++) {
      before_mass[xbin] -> SetBinContent(ybin, before->GetBinContent(xbin, ybin));
      after_mass[xbin] -> SetBinContent(ybin, after->GetBinContent(xbin, ybin));
      diff_mass[xbin] -> SetBinContent(ybin, before->GetBinContent(xbin, ybin)-after->GetBinContent(xbin, ybin));
    }
  }

  // signal
  TFile *mc = new TFile("pmcs.root");
  gDirectory->cd("smeared");
  TH2D *Sig_All = (TH2D *)gROOT->FindObject("smeared_ZMass_vs_DetEta_Before_HMx_CC");
  TH1D *Sig_before[NBins], *Sig_after[NBins];

  Int_t Sig_NBins_X = Sig_All->GetNbinsX();
  Double_t Sig_Axis_XMin = Sig_All->GetXaxis()->GetXmin();
  Double_t Sig_Axis_XMax = Sig_All->GetXaxis()->GetXmax();
  Int_t Sig_NBins_Y = before->GetNbinsY();
  Double_t Sig_Axis_YMin = before->GetYaxis()->GetXmin();
  Double_t Sig_Axis_YMax = before->GetYaxis()->GetXmax();

  if((NBins_X != Sig_NBins_X) || (NBins_Y != Sig_NBins_Y) || 
     (Axis_XMin != Sig_Axis_XMin) || (Axis_YMin != Sig_Axis_YMin) || 
     (Axis_XMax != Sig_Axis_XMax) || (Axis_YMax != Sig_Axis_YMax)) {
    cout<<"Different bin size between data and MC signal histogram, will exit now"<<endl;
    exit;
  }
  
  for (int xbin = 0; xbin < Sig_NBins_X; xbin ++){
    sprintf(name, "%s%d", "Sig_before_", xbin);
    Sig_before[xbin] = new TH1D(name, name, Sig_NBins_Y, Sig_Axis_YMin, Sig_Axis_YMax);
    sprintf(name, "%s%d", "Sig_after_", xbin);
    Sig_after[xbin] = new TH1D(name, name, Sig_NBins_Y, Sig_Axis_YMin, Sig_Axis_YMax);
    
    for(int ybin=0; ybin < Sig_NBins_Y; ybin ++) {
      Sig_before[xbin] -> SetBinContent(ybin, Sig_All->GetBinContent(xbin, ybin));
      Sig_after[xbin] -> SetBinContent(ybin, Sig_All->GetBinContent(xbin, ybin));
    }
  }

  cout<<"Sig HMx bin "<<Sig_NBins_X<<" from "<<Sig_Axis_XMin<<" to "<<Sig_Axis_XMax<<endl;
  cout<<"Sig Mass bin "<<Sig_NBins_Y<<" from "<<Sig_Axis_YMin<<" to "<<Sig_Axis_YMax<<endl;

  //
  // create a background histogram from the function
  //
  TF1 *func = new TF1("func", "[0] + [1]*x + [2]*x*x + [3]*x*x*x + [4]*x*x*x*x + [5]*x*x*x*x*x", Axis_YMin, Axis_YMax);
  func -> SetParameter(0, -28900.202091);
  func -> SetParameter(1, 1544.304929);
  func -> SetParameter(2, -31.665620);
  func -> SetParameter(3, 0.317247);
  func -> SetParameter(4, -0.001565);
//.........这里部分代码省略.........