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

void decon_ind(Int_t chan = 221){
  TFile *file = new TFile("Run_00009188.root");
  TH2* h1 = (TH2*)file->Get(Form("Ind_%d",chan));
  TH2* h3 = (TH2*)h1->Clone("h3");
  TH2* h4 = (TH2*)h1->Clone("h4");
  
  h3->Reset();
  h4->Reset();

  gStyle->SetOptStat(0);

  TCanvas *c1 = new TCanvas("c1","c1",800,800);
  c1->Divide(2,2);
  c1->cd(1);
  h1->Draw("COLZ");
  h1->GetZaxis()->SetRangeUser(-80,80);
  h1->GetXaxis()->SetRangeUser(0,3000);

  h1->SetTitle("Original Induction");
  h3->SetTitle("Noise Removed Induction");
  
  Int_t max_bin = h1->GetMaximumBin();
  Float_t max = h1->GetBinContent(max_bin);
  Int_t min_bin = h1->GetMinimumBin();
  Float_t min = h1->GetBinContent(min_bin);
  
  TH1F *h2 = new TH1F("h2","h2",Int_t(max - min+2), min-1, max+1);
  
  //cout << h1->GetXaxis()->GetNbins() << " " << h1->GetYaxis()->GetNbins() << endl;
  for (Int_t i=0;i!=h1->GetXaxis()->GetNbins();i++){
    h2->Reset();
    for (Int_t j=0;j!=h1->GetYaxis()->GetNbins();j++){
      h2->Fill(h1->GetBinContent(i,j));
    }
    Double_t xq = 0.5;
    Double_t par[10];
    h2->GetQuantiles(1,&par[0],&xq);
    
    for (Int_t j=0;j!=h1->GetYaxis()->GetNbins();j++){
      h3->SetBinContent(i,j,h1->GetBinContent(i,j)-par[0]);
    }
    //cout << par[0] << endl;
    
  }
  c1->cd(2);
  h3->Draw("COLZ");
  h3->GetXaxis()->SetRangeUser(0,3000);
  

  // do deconvolution for every channel ...
  c1->cd(3);
  TFile *file1 = new TFile("ave_res.root");
  TH1F *hva = (TH1F*)file1->Get("hu");
  TGraph *gva = new TGraph();
  for (Int_t i=0;i!=hva->GetNbinsX();i++){
    Double_t x = hva->GetBinCenter(i+1);
    Double_t y = hva->GetBinContent(i+1) * (-1);
    gva->SetPoint(i,x,y);
  }
  TH1F *h2t = new TH1F("h2t","h2t",h1->GetXaxis()->GetNbins(),0,h1->GetXaxis()->GetNbins());
  TH1F *h2r = (TH1F*)h2t->Clone("h2r");
  h2r->Reset();
  for (Int_t i=0;i!=h2r->GetNbinsX();i++){
     Double_t x = h2r->GetBinCenter(i+1)/2.-50;
     h2r->SetBinContent(i+1,gva->Eval(x));
  }
  TF1 *filter_v = new TF1("filter_v","(x>0.0)*gaus*exp(-0.5*pow(x/[3],[4]))");
  double par1[5]={1.74/0.941034, 1.46, 1.33, 0.23, 4.89};
  filter_v->SetParameters(par1);

  TVirtualFFT::SetTransform(0);
  TH1 *hmr = 0;
  TH1 *hpr = 0;
  Int_t n = 8192;
  TVirtualFFT *ifft = TVirtualFFT::FFT(1,&n,"C2R M K");
  Double_t value_re[8192];
  Double_t value_im[8192];
  TH1 *fb = 0;
  TH1 *hmv = 0;
  TH1 *hpv = 0;

  for (Int_t k=0;k!=h1->GetYaxis()->GetNbins();k++){
  // for (Int_t k=33;k!=33+1;k++){
    h2t->Reset();
    for (Int_t i=0;i!=8192;i++){
      float content = h3->GetBinContent(i+1,k+1);
      h2t->SetBinContent(i+1,content);
    }
    h2t->Draw();

    cout << h2t->Integral(500,1500) << endl;
    //cout << max << " " << min << endl;
    hmr = 0;
    hpr = 0;
    hmr = h2r->FFT(hmr,"MAG");
    hpr = h2r->FFT(hpr,"PH");
    
    hmv = 0;
    hpv = 0;
    hmv = h2t->FFT(hmv,"MAG");
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	  fabs(caloJet1Eta_) < 2.5 &&
	  fabs(caloJet2Eta_) < 2.5 &&
	  caloDeltaEta_ < 1.3)
	{
	  caloMjjSpectrum->Fill(caloMjj);
	  mjj200_eff->Fill((caloMjj>200 && l1Accept->at(L1scenario)==1) || hltAccept->at(HT250Calo)==1, caloMjj);
	  calo250_eff->Fill((hltAccept->at(HT250Calo)==1 && l1Accept->at(L1scenario)==1), caloMjj);

	  //references
	  HTT240_eff->Fill(l1Accept->at(HTT240)==1, caloMjj);
	  //l1 and hlt rates
	  for(unsigned int ii=0; ii<l1Names->size(); ++ii)
	    if (l1Accept->at(ii)==1)
	      l2->Fill(ii);
	}

      //PF analysis
      if (PFJet1Pt_ > 60. &&
	  PFJet2Pt_ > 30. &&
	  fabs(PFJet1Eta_) < 2.5 &&
	  fabs(PFJet2Eta_) < 2.5 &&
	  PFDeltaEta_ < 1.3)
	{
	  PFMjjSpectrum->Fill(PFMjj);
	  mjj450_eff->Fill((caloMjj>450 && l1Accept->at(L1scenario)==1) || hltAccept->at(HT410PF)==1, PFMjj);
	  pf410_eff->Fill((hltAccept->at(HT410PF)==1 && l1Accept->at(L1scenario)==1), PFMjj);
	}
    }

  mjj450_eff->Fit(f1,"r");
  mjj200_eff->Fit(f2,"r");


  caloMjjSpectrum->Scale(1./caloMjjSpectrum->GetBinContent(caloMjjSpectrum->GetMaximumBin()));
  PFMjjSpectrum->Scale(1./PFMjjSpectrum->GetBinContent(PFMjjSpectrum->GetMaximumBin()));

			      
  
  TLegend* leg0 = new TLegend(0.62, 0.78, 0.83, 0.89);
  leg0->AddEntry(mjj200_eff,"MJJ200Calo || HT250Calo","L");
  leg0->AddEntry(calo250_eff,"HT250_Calo","P");
  leg0->AddEntry(HTT240_eff,"HTT240","P");

  TLegend* leg1 = new TLegend(0.62, 0.78, 0.83, 0.89);
  leg1->AddEntry(mjj450_eff,"MJJ450PF || HT410PF","L");
  leg1->AddEntry(pf410_eff,"HT410_PF","P");

  TCanvas* c1 = new TCanvas();
  mjj200_eff->Draw();
  calo250_eff->Draw("sames");
  HTT240_eff->Draw("sames");
  caloMjjSpectrum->Draw("L,sames");
  leg0->Draw("sames");

  TCanvas* c2 = new TCanvas();
  mjj450_eff->Draw();
  pf410_eff->Draw("sames");
  PFMjjSpectrum->Draw("L,sames");
  leg1->Draw("sames");

  TCanvas* c3 = new TCanvas();
  //l1->Scale(PDRate/nentries);

  for(unsigned int ii=0; ii<l1Names->size(); ++ii)
    l1->GetXaxis()->SetBinLabel(ii+1,l1Names->at(ii).c_str());
  //l1->GetYaxis()->SetTitle("L1 Rate @4E33 [Hz]");

void analysis() {
    Int_t nbins = 800;
    Int_t j;
    char name[20];
    char title[100];
    TH1F *HistoEvent[2214];
    for (Int_t z=0;z<2214;z++) {
        sprintf(name,"HistoEvent%d",z-1);
        sprintf(title,"Event%d Histo", z-1);
        HistoEvent[z] = new TH1F(name,title,nbins, -0.1, 159.9);
    }
    TH1F *NewHistoEvent[2214];
    for (Int_t z=0;z<2214;z++) {
        sprintf(name,"NewHistoEvent%d",z-1);
        sprintf(title,"Event%d Histo", z-1);
        NewHistoEvent[z] = new TH1F(name,title,nbins, -0.1, 159.9);
    }
    TH1F *NewHistoEventFFT[2214];
    for (Int_t z=0;z<2214;z++) {
        sprintf(name,"NewHistoEventFFT%d",z-1);
        sprintf(title,"Event%d Histo", z-1);
        NewHistoEventFFT[z] = new TH1F(name,title,nbins, 0, 5);
    }
    Double_t mean;
    Double_t rms;
    Double_t meansum = 0;
    Double_t count = 0;
    Double_t meanrms = 0;
	TFile f("/home/marko/H4Analysis/ntuples/analysis_4443.root"); //ntuple generated by H4Analysis tool
	TFile f1("/home/marko/H4Analysis/ntuples/analysis_3905.root");
    TFile f2("/home/marko/Desktop/TB Timing Res/NormalizedSignalNoise.root", "read");
    TH1F* BestSignal = (TH1F*) f2.Get("BetterSignal");
    TFile outputfile("myoutput.root", "recreate");
    TCanvas* TimeandFreq = new TCanvas("TimeandFreq","Time and Frequency",1500,900);
    TCanvas* Freq = new TCanvas("Freq","Frequency",800,1200);
    TCanvas* TimeSignal = new TCanvas("TimeSignal","Pure Signal",800,1200);
	TimeandFreq->Divide(2,2);
	TTree* h4 = (TTree*) f.Get("h4");
    TTree* h4_2 = (TTree*) f1.Get("h4");
    TString plot;
    TString cut;
	TH2F* WavePulse = new TH2F ("WavePulse", "Wave Pulse", nbins, -0.1, 159.9, 850, -50, 800);
    TH2F* NoisePulse = new TH2F ("NoisePulse", "Noise", nbins, -0.1, 159.9, 100, -50, 50);
    TH1F* PulseTime = new TH1F ("PulseTime", "Original Wave Pulse", nbins, -0.1, 159.9); //nanoseconds
    TH2F* TempHisto = new TH2F ("TempHisto", "Temp Histo", nbins, -0.1, 159.9, 1000, -15, 15); //nanoseconds
	h4->Draw("WF_val:WF_time>>WavePulse", "WF_ch==2 && event==1 && spill==1");
    h4_2->Draw("WF_val:WF_time>>NoisePulse","WF_ch==APD1 && amp_max[APD3]<25 && b_rms[APD3]<5. && charge_tot[APD3]<20000 &&  amp_max[APD5]<25 && b_rms[APD5]<5. &&  amp_max[APD6]<25 && b_rms[APD6]<5. &&  amp_max[APD4]<25 && b_rms[APD4]<5. && amp_max[SiPM1]<20 && amp_max[SiPM2]<20 && amp_max[APD1]<40 && amp_max[APD2]<40 && b_rms[APD1]<5. && b_rms[APD2]<5. && WF_time<160");
    for (Int_t i=0; i<nbins; i++) {
		for (Int_t k=0; k<4096; k++) {
			if (WavePulse->GetBinContent(i+1, k) != 0) {
				PulseTime->SetBinContent(i+1,k-50);
			}
		}
	}
    TH1F *NoiseTime = new TH1F ("NoiseTime", "Noise", nbins, -0.1, 159.9);
    for (Int_t i=0; i<nbins; i++) {
        for (Int_t k=0; k<4096; k++) {
            if (NoisePulse->GetBinContent(i+1, k) != 0) {
                NoiseTime->SetBinContent(i+1,k-62.9087);
            }
        }
    }
    //TH1F* NormNoiseFFT = new TH1F ("NormNoiseFFT", "Normalized Noise FFT", nbins, 0, 5);
    //TStopwatch t;
    //t.Start(); //1 hour runtime
    //for (j=10;j<20;j++) {
    //    plot = "WF_val:WF_time>>TempHisto";
    //    cut = "WF_ch==APD1 && amp_max[APD3]<25 && b_rms[APD3]<5. && charge_tot[APD3]<20000 &&  amp_max[APD5]<25 && b_rms[APD5]<5. &&  amp_max[APD6]<25 && b_rms[APD6]<5. &&  amp_max[APD4]<25 && b_rms[APD4]<5. && amp_max[SiPM1]<20 && amp_max[SiPM2]<20 && amp_max[APD1]<40 && amp_max[APD2]<40 && b_rms[APD1]<5. && b_rms[APD2]<5. && WF_time<160 && event==";
    //    cut += j;
    //    h4_2->Draw(plot, cut, "goff");
    //    if (TempHisto->GetMaximum() == 0) {
    //        delete HistoEvent[j+1];
    //        continue;
    //    }
    //    for (Int_t i=0; i<nbins; i++) {
    //        for (Int_t k=0; k<1000; k++) {
    //            if (TempHisto->GetBinContent(i+1, k) != 0) {
    //                HistoEvent[j+1]->SetBinContent(i+1,k*0.03-15);
    //            }
    //        }
    //    }
    //    mean = TempHisto->GetMean(2);
    //    rms = TempHisto->GetRMS(2);
    //    for (Int_t q=0;q<nbins;q++) {
    //        NewHistoEvent[j+1]->SetBinContent(q+1, HistoEvent[j+1]->GetBinContent(q+1)-mean);
    //    }
    //    NewHistoEvent[j+1]->Scale(1/rms);
    //    NewHistoEvent[j+1]->FFT(NewHistoEventFFT[j+1], "MAG");
    //    NormNoiseFFT->Add(NormNoiseFFT, NewHistoEventFFT[j+1]);
    //    TempHisto->Write();
    //    NewHistoEvent[j+1]->Write();
    //    NewHistoEventFFT[j+1]->Write();
    //    cout << "Event " << j << ", Mean = " << mean << ", RMS = " << rms << endl;
    //    count += 1;
    //}
    //NormNoiseFFT->Scale(1/count);
    //NormNoiseFFT->Write();
    //t.Stop();
    //t.Print();
	new TFile("/home/marko/H4Analysis/ntuples/analysis_4443.root"); // ignore this reloading of the same file, it is required or else the plots do not show up (when I tried)
//.........这里部分代码省略.........

void makeTable(TString myVar, TString myCut, TString myName, TString myAxisNameX, TString myAxisNameY, 
               vector<const Sample*>& listOfSamples, vector<const Sample*> listOfDatasets, TString inFileName,
               bool isBlind, bool isLog, 
               int nBins, float xLow, float xHigh,
               float* xlowVec)
{
  // prepare the input file
  TFile* infile = new TFile(inFileName, "READ"); 
  infile -> cd();

  // prepare the necessary for the scale factor estimation
  float defScaleFactor = 1.3;
  float newScaleFactor = 1.;
  float varScaleFactor = 1.;
  float nBkg = 0;
  float nBkgErrSq = 0;
  float nSig = 0;
  float nSigErrSq = 0;
    
  // prepare the stack
  THStack *hs = new THStack("hs","");
  // prepare the histos pointers
  TH1F*   hist[20];
  // prepare the tree pointers
  TTree*  tree[20];
  // prepare the legend
  TLegend* leg = new TLegend(.7485,.7225,.9597,.9604);
  leg->SetFillColor(0);
  // prepare the colors
  Int_t col[20] = {46,2,12,5,3,4,9,7,47,49,49,50,51,52,53,54,55,56,57,58};
  // prepare the Y axis lable
  if (xlowVec != 0) myAxisNameY = "Events/" + myAxisNameY;
  else {
    float binWidth = (xHigh-xLow)/nBins;
    TString tempString;
    tempString.Form("%.2f ",binWidth); 
    myAxisNameY = "Events/" + tempString + myAxisNameY;
  }
  // prepare the legend strings
  vector<TString> theLegends;
  
  // loop through the datasets and produce the plots
  TH1F* hdata;
  //variable bin histo
  if (xlowVec != 0) hdata = new TH1F("hdata","",nBins,xlowVec);
  //fixed bin histo
  else hdata = new TH1F("hdata","",nBins,xLow,xHigh);
  
  TTree*  treedata[20];
  for (UInt_t iDatas=0; iDatas < listOfDatasets.size(); iDatas++) {
    //get the tree
    treedata[iDatas] = (TTree*) infile -> Get(listOfDatasets.at(iDatas)->Name()->Data());

    //fill the histogram
    if ( iDatas == 0 ) treedata[iDatas] -> Draw(myVar + " >> hdata","hlt_weight*evt_weight*kf_weight*pu_weight" + myCut,"goff");
    else treedata[iDatas] -> Draw(myVar + " >>+ hdata","hlt_weight*evt_weight*kf_weight*pu_weight" + myCut,"goff");
    
    if ( iDatas == 0 ) leg -> AddEntry(hdata, "DATA (19.8 fb^{-1})", "pl");    
    
  }//end loop on datasets
  if (xlowVec != 0) {
    for (int iBin = 1; iBin <= nBins; iBin++) hdata->SetBinError  (iBin,hdata->GetBinError(iBin)/hdata->GetBinWidth(iBin));
    for (int iBin = 1; iBin <= nBins; iBin++) hdata->SetBinContent(iBin,hdata->GetBinContent(iBin)/hdata->GetBinWidth(iBin));
  }
       
  int theHistCounter = 0;
  // loop through the samples and produce the plots
  for (UInt_t iSample=0; iSample < listOfSamples.size(); iSample++) {

    //determine if the histo is first of the series
    bool isFirstOfSerie = (*listOfSamples.at(iSample)->Legend()).CompareTo(" ");
    bool isLastOfSerie = false;
    if (iSample == listOfSamples.size() - 1) isLastOfSerie = true;
    if (iSample < listOfSamples.size() - 1 && (*listOfSamples.at(iSample+1)->Legend()).CompareTo(" ") != 0) isLastOfSerie = true;
    
    //get the tree
    tree[iSample] = (TTree*) infile -> Get(listOfSamples.at(iSample)->Name()->Data());
    //if sample first of the list create a new histogram
    if (isFirstOfSerie) {
       TString thisHistName = "h_" + *(listOfSamples.at(iSample)->Name());
       //variable bin histo
       if (xlowVec != 0) hist[theHistCounter] = new TH1F(thisHistName,thisHistName,nBins,xlowVec);
       //fixed bin histo
       else hist[theHistCounter] = new TH1F(thisHistName,thisHistName,nBins,xLow,xHigh);
       hist[theHistCounter] -> Sumw2();
       hist[theHistCounter] -> SetFillColor(col[theHistCounter]);
       hist[theHistCounter] -> SetFillStyle(1001);
       theLegends.push_back(*listOfSamples.at(iSample)->Legend());
       cout << *listOfSamples.at(iSample)->Legend() << " ";
    }

    //fill the histogram
    
    TString thisScale = Form("%f *", *(listOfSamples.at(iSample)->Scale()));
    if (isFirstOfSerie) tree[iSample] -> Draw(myVar + " >> " + TString(hist[theHistCounter] -> GetName()),thisScale + "hlt_weight*evt_weight*kf_weight*pu_weight" + myCut,"goff");
    else tree[iSample] -> Draw(myVar + " >>+ " + TString(hist[theHistCounter] -> GetName()),thisScale + "hlt_weight*evt_weight*kf_weight*pu_weight" + myCut,"goff");
        
    //add the histogram to the stack if the last of the series:
    //either last sample or ~ sample followed by non ~ sample
    if (isLastOfSerie) {
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	  if ( theEvent->energies[0]>RefScintThresh){
	    TOFEnergyNeutrons1PSL_L->Fill(theEvent->TOFEnergy);
	    TOFEnergyCFBinsL_L->Fill(theEvent->TOFEnergy);
	    EvsTOFEnergyL_L->Fill(theEvent->TOFEnergy,theEvent->energies[0]);
	  }
	}
	  

      
	if (theEvent->ShiftTOF < -2 && theEvent->ShiftTOF>-35){
	  if ( theEvent->energies[0]>RefScintThresh){
	    TOFEnergyRandomBkg1PSL_L->Fill(theEvent->TOFEnergy);
	    TOFEnergyBkgCFBinsL_L->Fill(theEvent->TOFEnergy);
	    EvsTOFEnergyBkgL_L->Fill(theEvent->TOFEnergy,theEvent->energies[0]);
	  }
	}
      }
    }



    Multiplicity->Fill(theEvent->N);


    if ( jentry % 100000 == 0){
      cout<<"ON "<<jentry<<endl;
    }

  }//End Main Loop over Everything

  ///Do analysis on whole historgrams 
  /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  for (int i=1;i<=TOFEnergyNeutrons->GetXaxis()->GetNbins();i++){
    Float_t N = TOFEnergyNeutrons->GetBinContent(i);
    Float_t NRandom = TOFEnergyRandomBkg->GetBinContent(i);
    TOFEnergyRBkgSubtracted->SetBinContent(i,N-NRandom);
    TOFEnergyRBkgSubtracted->SetBinError(i,TMath::Sqrt(N+NRandom));


    Float_t NL_L =TOFEnergyNeutrons1PSL_L->GetBinContent(i);
    Float_t NRandomL_L=TOFEnergyRandomBkg1PSL_L->GetBinContent(i);
    TOFEnergyRBkgSubtracted1PSL_L->SetBinContent(i,NL_L-NRandomL_L);
    TOFEnergyRBkgSubtracted1PSL_L->SetBinError(i,TMath::Sqrt(NL_L+NRandomL_L));
  }
  
  //Spectra that have the CF BINNING
  for (int i=1;i<=TOFEnergySubCFBins->GetXaxis()->GetNbins();i++){
    Float_t N =TOFEnergyCFBins->GetBinContent(i);
    Float_t NRandom=TOFEnergyBkgCFBins->GetBinContent(i);
    TOFEnergySubCFBins->SetBinContent(i,N-NRandom);
    TOFEnergySubCFBins->SetBinError(i,TMath::Sqrt(N+NRandom));
    
    Float_t NL_L = TOFEnergyCFBinsL_L->GetBinContent(i);
    Float_t NRandomL_L=TOFEnergyBkgCFBinsL_L->GetBinContent(i);
    TOFEnergySubCFBinsL_L->SetBinContent(i,NL_L-NRandomL_L);
    TOFEnergySubCFBinsL_L->SetBinError(i,TMath::Sqrt(NL_L+NRandomL_L));
    
    for (int a=0;a<LendaThreshs.size();a++){
      Float_t thisN=LendaThreshs[a]->GetBinContent(i);
      Float_t thisNRandom=LendaThreshsBkg[a]->GetBinContent(i);
      LendaThreshsResult[a]->SetBinContent(i,thisN-thisNRandom);
      LendaThreshsResult[a]->SetBinError(i,sqrt(thisN+thisNRandom));
    }
  }

  ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

///
/// Make a plot out of a 1D histogram holding a 1-CL curve.
/// The strategy is to always convert the 1-CL histogram (hCL) into
/// a TGraph. This way we can add known points (solutions, points
/// at end of scan range) and also have a filled area without line
/// smoothing. This is not possible with histograms due to a Root bug.
///
/// The function draws the TGraphs, and returns a pointer to the
/// TGraph object that can be used in the TLegend.
///
/// Markers are plotted if the method name of the scanner is "Plugin" or "BergerBoos" or "DatasetsPlugin".
/// One can plot a line instead of points even for the Plugin method by
/// using setPluginMarkers().
///
/// For the angle variables, a new axis is painted that is in Deg.
///
/// \param s The scanner to plot.
/// \param first
/// \param last
/// \param filled
///
TGraph* OneMinusClPlot::scan1dPlot(MethodAbsScan* s, bool first, bool last, bool filled, int CLsType)
{
	if ( arg->debug ){
		cout << "OneMinusClPlot::scan1dPlot() : plotting ";
		cout << s->getName() << " (" << s->getMethodName() << ")" << endl;
	}
	if ( m_mainCanvas==0 ){
		m_mainCanvas = newNoWarnTCanvas(name+getUniqueRootName(), title, 800, 600);
	}
	m_mainCanvas->cd();
	bool plotPoints = ( s->getMethodName()=="Plugin" || s->getMethodName()=="BergerBoos" || s->getMethodName()=="DatasetsPlugin" ) && plotPluginMarkers;
	TH1F *hCL = (TH1F*)s->getHCL()->Clone(getUniqueRootName());
	if (CLsType==1) hCL = (TH1F*)s->getHCLs()->Clone(getUniqueRootName());
  else if (CLsType==2) hCL = (TH1F*)s->getHCLsFreq()->Clone(getUniqueRootName());
	// fix inf and nan entries
	for ( int i=1; i<=s->getHCL()->GetNbinsX(); i++ ){
		if ( s->getHCL()->GetBinContent(i)!=s->getHCL()->GetBinContent(i)
				|| std::isinf(s->getHCL()->GetBinContent(i)) ) s->getHCL()->SetBinContent(i, 0.0);
	}

	// remove errors the hard way, else root ALWAYS plots them
	if ( !plotPoints ) hCL = histHardCopy(hCL, true, true);

	// disable any statistics box
	hCL->SetStats(0);

	// Convert the histogram into a TGraph so we can add the solution.
	// Also, the lf2 drawing option is broken in latest root versions.
	TGraph *g;
	if ( plotPoints ) g = new TGraphErrors(hCL->GetNbinsX());
	else              g = new TGraph(hCL->GetNbinsX());
	g->SetName(getUniqueRootName());
	for ( int i=0; i<hCL->GetNbinsX(); i++ ){
		g->SetPoint(i, hCL->GetBinCenter(i+1), hCL->GetBinContent(i+1));
		if ( plotPoints ) ((TGraphErrors*)g)->SetPointError(i, 0.0, hCL->GetBinError(i+1));
	}

	// add solution
	if ( ! s->getSolutions().empty() ){
		TGraphTools t;
		TGraph *gNew = t.addPointToGraphAtFirstMatchingX(g, s->getScanVar1Solution(0), 1.0);
		delete g;
		g = gNew;
	}

	// // set last point to the same p-value as first point by hand
	// // some angle plots sometimes don't manage to do it by themselves...
	// if ( arg->isQuickhack(XX) )
	// {
	//   Double_t pointx0, pointy0, err0;
	//   Double_t pointx1, pointy1, err1;
	//   g->GetPoint(0, pointx0, pointy0);
	//   g->GetPoint(g->GetN()-1, pointx1, pointy1);
	//   g->SetPoint(g->GetN()-1, pointx1, pointy0);
	//   if ( plotPoints ) err0 = ((TGraphErrors*)g)->GetErrorY(0);
	//   if ( plotPoints ) ((TGraphErrors*)g)->SetPointError(g->GetN()-1, 0.0, err0);
	// }

	// add end points of scan range
	if ( !plotPoints )
	{
		Double_t pointx0, pointy0;
		TGraph *gNew = new TGraph(g->GetN()+4);
		gNew->SetName(getUniqueRootName());
		for ( int i=0; i<g->GetN(); i++)
		{
			g->GetPoint(i, pointx0, pointy0);
			gNew->SetPoint(i+2, pointx0, pointy0);
		}

		// add origin
		gNew->SetPoint(0, hCL->GetXaxis()->GetXmin(), 0);

		// add a point at first y height but at x=origin.
		g->GetPoint(0, pointx0, pointy0);
		gNew->SetPoint(1, hCL->GetXaxis()->GetXmin(), pointy0);

		// add a point at last y height but at x=xmax.
		g->GetPoint(g->GetN()-1, pointx0, pointy0);
//.........这里部分代码省略.........

void GenerateInputFile_ge4b( double mgl=-1., double mlsp=-1., double target_susy_all0lep=-1. ) {

  TChain* dyTree = new TChain("treeZ") ;
  int nAdded = dyTree->Add("files15fb_8TeV_old1/DY.root") ;
  if ( nAdded <= 0 ) {
     printf("\n\n\n *** No treeZ in files15fb_8TeV_old1/DY.root\n\n\n") ;
     return ;
  }

  double t1bbbbWeight(0.) ;
  TChain chainT1bbbb("tree") ;
  char susycutstring[1000] ;
  sprintf( susycutstring, "&&mgluino==%.0f&&mlsp==%.0f", mgl, mlsp ) ;
  TString susycut( susycutstring ) ;
  if ( mgl>0. && mlsp>0. ) {
     nAdded = chainT1bbbb.Add("files5fb_MT/T1bbbb.root") ;
     if ( nAdded <= 0 ) {
        printf("\n\n\n *** No tree in files5fb_MT/T1bbbb.root\n\n\n") ;
        return ;
     }
     TFile f("referenceXSecs.root") ;
     TH1F* xsechist = (TH1F*) f.Get("gluino8TeV_NLONLL") ;
     if ( xsechist==0x0 ) { printf("\n\n *** can't find reference Xsec histogram in referenceXSecs.root.\n\n") ; return ; }
     int theBin = xsechist->FindBin( mgl ) ;
     if ( theBin <=0 || theBin > xsechist->GetNbinsX() ) {
        printf("\n\n *** can't find bin for mgl=%g.  Returned %d\n\n", mgl, theBin ) ;
        return ;
     }
     double xsec = xsechist->GetBinContent( theBin ) ;
     printf("\n\n T1bbbb xsec for mgl=%g is %g\n\n", mgl, xsec ) ;
     t1bbbbWeight = 1.5*xsec ;  //in pb. scan has 10k events, so nScan*1.5*xsec = events in 15fb-1
     //////  t1bbbbWeight = 0.1*xsec ;  //in pb. T1tttt scan has 50k events, so nScan*0.1*xsec = events in 5fb-1
     printf("\n\n Susy ttree cut: %s\n\n", susycutstring ) ;
  }


  TChain chainQCD("tree") ;
   //--- these have high weight
//chainQCD.Add("files15fb_8TeV/QCD-50to80.root");
//chainQCD.Add("files15fb_8TeV/QCD-80to120.root");
  chainQCD.Add("files15fb_8TeV/QCD-120to170.root");
  chainQCD.Add("files15fb_8TeV/QCD-170to300.root");
   //--- below here, these have weight less than one.
  chainQCD.Add("files15fb_8TeV/QCD-300to470.root");
  chainQCD.Add("files15fb_8TeV/QCD-470to600.root");
  chainQCD.Add("files15fb_8TeV/QCD-600to800.root");
  chainQCD.Add("files15fb_8TeV/QCD-800to1000.root");
  chainQCD.Add("files15fb_8TeV/QCD-1000to1400.root");
  chainQCD.Add("files15fb_8TeV/QCD-1400to1800.root");
  chainQCD.Add("files15fb_8TeV/QCD-1800.root");

  TChain chainTT("tree") ;
  chainTT.Add("files15fb_8TeV/TT.root") ;

  TChain chainZnn("tree") ;
  chainZnn.Add("files15fb_8TeV/Zinv-100to200.root") ;
  chainZnn.Add("files15fb_8TeV/Zinv-200to400.root") ;
  chainZnn.Add("files15fb_8TeV/Zinv-400.root") ;

  TChain chainWJets("tree") ;
  chainWJets.Add("files15fb_8TeV/WJets-250to300.root") ;
  chainWJets.Add("files15fb_8TeV/WJets-300to400.root") ;
  chainWJets.Add("files15fb_8TeV/WJets-400.root") ;
  chainWJets.Add("files15fb_8TeV/T-s.root") ;
  chainWJets.Add("files15fb_8TeV/T-t.root") ;
  chainWJets.Add("files15fb_8TeV/T-tW.root") ;
  chainWJets.Add("files15fb_8TeV/Tbar-s.root") ;
  chainWJets.Add("files15fb_8TeV/Tbar-t.root") ;
  chainWJets.Add("files15fb_8TeV/Tbar-tW.root") ;



      char qcdinputfile[9][1000] = {
        "files15fb_8TeV/QCD-120to170.root"
       ,"files15fb_8TeV/QCD-170to300.root"
       ,"files15fb_8TeV/QCD-300to470.root"
       ,"files15fb_8TeV/QCD-470to600.root"
       ,"files15fb_8TeV/QCD-600to800.root"
       ,"files15fb_8TeV/QCD-800to1000.root"
       ,"files15fb_8TeV/QCD-1000to1400.root"
       ,"files15fb_8TeV/QCD-1400to1800.root"
       ,"files15fb_8TeV/QCD-1800.root"
      } ;

      char qcdsamplename[9][100] = {
        "qcd_0120_to_0170"
       ,"qcd_0170_to_0300"
       ,"qcd_0300_to_0470"
       ,"qcd_0470_to_0600"
       ,"qcd_0600_to_0800"
       ,"qcd_0800_to_1000"
       ,"qcd_1000_to_1400"
       ,"qcd_1400_to_1800"
       ,"qcd_1800_to_9999"
      } ;



  gROOT->Reset();

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

//.........这里部分代码省略.........
  my_canvas2->cd();
  gPad->SetLogy();
  total->SetMaximum(1000);
  total->SetMarkerStyle(20);
  total->Draw("E");
  my_canvas2->Write();
  */
  TCanvas* ratio_c = new TCanvas("ratio","ratio_canvas",800,700);
  ratio_c->SetTopMargin(0.); 
  ratio_c->SetBottomMargin(0.); 
  ratio_c->Update();
  ratio_c->Divide(1, 2);
  ratio_c->cd(1);

  gPad->SetPad(.005, .30, .995, .995);
  TPad* pad1 = (TPad*)ratio_c->GetPad(1); 
  pad1->SetTopMargin(0.05);
  pad1->SetBottomMargin(0.05); 
  
  ratio_c->Update();
  gPad->SetLogy();
  
  hs->Draw("HIST");                                                                                                                                        
  hs->SetMaximum(100000);
  hs->SetMinimum(0.01);
  //  hs->GetXaxis()->SetLimits(40, 3200);
  hs->GetXaxis()->SetRangeUser(0, 1000);
  hs->GetXaxis()->SetTitle("");
  // total->Draw("SAME E2");
  mydata->Draw("SAME E1");
  tauPt_Stage1_Wprime_M4000->Draw("SAME HIST");
  CMS_text->Draw("same");
  CMS_text_2->Draw("same");
  lumiText->Draw("same");
  leg_example->Draw("same");
  ratio_c->cd(2);

  gPad->SetPad(.005, .08, .995, .28); 
  gStyle->SetOptStat(false);
  ratio_c->Update();
  TPad* pad2 = (TPad*)ratio_c->GetPad(2);
  pad2->SetTopMargin(0.04); 
  pad2->SetBottomMargin(0.30); 
  pad2->SetGridx();
  pad2->SetGridy();

  std::cout << "data bins = " << mydata->GetNbinsX() << " width=" << mydata->GetBinWidth(2);
  std::cout << " MC bins = " << total->GetNbinsX() << " width=" << total->GetBinWidth(2) << std::endl;

  int nbin=mydata->GetNbinsX() ;
  float width=mydata->GetBinWidth(2);
  TH1F *data_by_MC = new TH1F("h1", "ratio", nbin, 0, 4000);
  std::cout << "before D/MC bins=" << data_by_MC->GetNbinsX() << " width=" << data_by_MC->GetBinWidth(2) << std::endl;

  for (int i=0; i<nbin; i++) {
    float data = mydata->GetBinContent(i);
    float MC = total->GetBinContent(i);
    float ratio=0;
    if (MC>0) ratio=data/MC ;
    data_by_MC->SetBinContent(i,ratio);
    if (data>0) {
      float stat_err = sqrt(data)/data ;  
      data_by_MC->SetBinError(i,stat_err);
    }

    // data_by_MC->SetBinWidth(i,width);  
  }

  //  TH1F* data_by_MC = (TH1F*)mydata->Clone();
  // data_by_MC->Divide(total);
  data_by_MC->SetMarkerStyle(20);
  data_by_MC->Draw("E");
  data_by_MC->GetXaxis()->SetTitle("#tau pT [GeV]");
  data_by_MC->GetYaxis()->SetTitle("#frac{DATA}{MC}");
  data_by_MC->GetYaxis()->SetLabelSize(0.12);
  data_by_MC->GetXaxis()->SetLabelSize(0.12);
  data_by_MC->GetYaxis()->SetTitleSize(0.18);
  data_by_MC->GetXaxis()->SetTitleSize(0.18);
  data_by_MC->GetYaxis()->SetTitleOffset(0.20);
  data_by_MC->GetXaxis()->SetTitleOffset(0.75);
  data_by_MC->SetTitle("");
  data_by_MC->GetXaxis()->SetRangeUser(0, 1000);
  data_by_MC->SetMaximum(4);
  data_by_MC->SetMinimum(0);
  data_by_MC->GetYaxis()->SetNdivisions(4);



  TLine *l=new TLine(0,1,1000,1);
  l->SetLineColor(kRed);
  l->Draw("same");

  std::cout << "D/MC bins=" << data_by_MC->GetNbinsX() << " width=" << data_by_MC->GetBinWidth(2) << std::endl;
  ratio_c->Write();
  ratio_c->Print("tauPt_Stage1_ratio_DD.pdf");

 
  return 0;

}

void Skim(TString fname="ZnnH125", TString outputname = "")
{
    gROOT->LoadMacro("HelperFunctions.h" );  // make functions visible to TTreeFormula
    gROOT->SetBatch(1);

    TChain * chain  = new TChain("tree");
    TString dijet   = "";
    TString outdir  = "/afs/cern.ch/work/d/degrutto/public/MiniAOD/ZnnHbb_Phys14_PU20bx25/skimV11_v2/";
    TString prefix  = "skim_";
    TString suffix  = "tree*.root";
    


    TCut selection = baseline.c_str();
    // Different baseline for dimuons
    // MET filters
    selection += metfilter.c_str();

    // JSON & trigger
    if (fname.Contains("Data")) {
        TCut trigger = mettrigger.c_str();
        selection += trigger;
        //selection.Print();
      }
      /*
    } else if (process == "Data_METBTag_R" || process == "Data_METBTag_P") {
        TCut trigger = metbtagtrigger.c_str();
        selection += trigger;
        //selection.Print();
    } else if (process == "Data_SingleMu_R" || process == "Data_SingleMu_P") {
        TCut trigger = singlemutrigger.c_str();
        selection += trigger;
        //selection.Print();
    } else if (process == "Data_SingleEl_R" || process == "Data_SingleEl_P") {
        TCut trigger = singleeltrigger.c_str();
        selection += trigger;
        //selection.Print();
    }
      */


    chain->Add(fname);

    // Sum Count, CountWithPU, CountWithPUP, CountWithPUM
    TObjArray * files = chain->GetListOfFiles();
    TIter next(files);
    TChainElement * chainElem = 0;
    TFile * f2 = 0;
    TH1D * h1 = new TH1D("Count", ";Counts", 16, 0, 16);
    TH1F * htemp = 0;

    while ((chainElem = (TChainElement*) next())) {
      //#ifndef XROOTD
      //      f2 = TFile::Open("dcache:" + TString(chainElem->GetTitle()));
      //#else
      std::cout << "chainElem->GetTitle() " << chainElem->GetTitle() << std::endl; 
        f2 = TFile::Open( TString(chainElem->GetTitle()));
	//#endif
        htemp = (TH1F*) f2->Get("Count");
        h1->SetBinContent(1, h1->GetBinContent(1)+htemp->GetBinContent(1));
	/*
        htemp = (TH1F*) f2->Get("CountWithPU");
        h1->SetBinContent(2, h1->GetBinContent(2)+htemp->GetBinContent(1));
        htemp = (TH1F*) f2->Get("CountWithPUP");
        h1->SetBinContent(3, h1->GetBinContent(3)+htemp->GetBinContent(1));
        htemp = (TH1F*) f2->Get("CountWithPUM");
        h1->SetBinContent(4, h1->GetBinContent(4)+htemp->GetBinContent(1));
        htemp = (TH1F*) f2->Get("CountWithMCProd");
        h1->SetBinContent(5, h1->GetBinContent(5)+htemp->GetBinContent(1));
        htemp = (TH1F*) f2->Get("CountWithPUMCProd");
        h1->SetBinContent(6, h1->GetBinContent(6)+htemp->GetBinContent(1));
        htemp = (TH1F*) f2->Get("countWithSignalQCDcorrections");
        h1->SetBinContent(7, h1->GetBinContent(7)+htemp->GetBinContent(1));
        */
        std::clog << fname << ": skimmed from " << chainElem->GetTitle() << std::endl;
    }
    
    // LHE Count
    
    TH1D * h2 = new TH1D("LHECount", ";LHE Counts", 16, 0, 16);
    TString process_lhe = fname;
    if (process_lhe.BeginsWith("WJets"))
        process_lhe = "WJets";
    else if (process_lhe.BeginsWith("ZJets"))
        process_lhe = "ZJets";
    else 
        process_lhe = "";
    const std::vector<std::string>& lhecuts = GetLHECuts(process_lhe.Data());
    for (unsigned int i=0; i < lhecuts.size(); i++) {
        TCut cut2 = lhecuts.at(i).c_str();
        h2->SetBinContent(i+1, chain->GetEntries(cut2));
    }
    
    

    // Make output directory if it doesn't exist
    if (gSystem->AccessPathName(outdir))
        gSystem->mkdir(outdir);
    TString outname = outdir + prefix + Form("%s.root", outputname.Data());
    std::cout << "outname is " << outname << std::endl;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    qqH  ->Add(VH );
    ggH  ->Add(qqH);
#endif
#endif
  }

  /*
    Mass plot before and after fit
  */
  TCanvas *canv = MakeCanvas("canv", "histograms", 600, 600);

  canv->cd();
  if(log){ canv->SetLogy(1); }
#if defined MSSM
  if(!log){ data->GetXaxis()->SetRange(0, data->FindBin(345)); } else{ data->GetXaxis()->SetRange(0, data->FindBin(UPPER_EDGE)); };
#else
  data->GetXaxis()->SetRange(0, data->FindBin(345));
#endif
  data->SetNdivisions(505);
  data->SetMinimum(min);
#ifndef DROP_SIGNAL
  data->SetMaximum(max>0 ? max : std::max(std::max(maximum(data, log), maximum(Ztt, log)), maximum(ggH, log)));
#else
  data->SetMaximum(max>0 ? max : std::max(maximum(data, log), maximum(Ztt, log)));
#endif
  data->Draw("e");

  TH1F* errorBand = (TH1F*)Ztt ->Clone("errorBand");
  errorBand  ->SetMarkerSize(0);
  errorBand  ->SetFillColor(13);
  errorBand  ->SetFillStyle(3013);
  errorBand  ->SetLineWidth(1);
  for(int idx=0; idx<errorBand->GetNbinsX(); ++idx){
    if(errorBand->GetBinContent(idx)>0){
      std::cout << "Uncertainties on summed background samples: " << errorBand->GetBinError(idx)/errorBand->GetBinContent(idx) << std::endl;
      break;
    }
  }
  if(log){
    Ztt  ->Draw("histsame");
    ttbar->Draw("histsame");
    EWK  ->Draw("histsame");
    EWK1 ->Draw("histsame");
    Fakes->Draw("histsame");
#ifdef MSSM
    VH_SM125->Draw("histsame");
#endif   
    $DRAW_ERROR
#ifndef DROP_SIGNAL
    ggH  ->Draw("histsame");
#endif
  }
  else{
#ifndef DROP_SIGNAL
    ggH  ->Draw("histsame");
#endif
    Ztt  ->Draw("histsame");
    ttbar->Draw("histsame");
    EWK  ->Draw("histsame");
    EWK1 ->Draw("histsame");
    Fakes->Draw("histsame");
#ifdef MSSM
    VH_SM125->Draw("histsame");
#endif
    $DRAW_ERROR
  }

//.........这里部分代码省略.........
    {
        binLo = 1;
        binHi = nbins;
        xmin = hsherpa->GetBinLowEdge(1);
        xmax = hsherpa->GetBinLowEdge(nbins+1);
    }


    //   float scale_mc = (float)hsherpa->Integral(binLo,binHi)/(float)hmadgraph->Integral(binLo,binHi);
    //   cout << "binLo = " << binLo << ", binHi = " << binHi << endl;
    //   cout << "xmin = " << xmin << "xmax = " << xmax << endl;

    //   hmadgraph->Sumw2();
    //   hmadgraph->Scale(scale_mc);

    float scale_sherpa   = 1000.0*4.890*3048.0/4.71644910071102437e+06;

    float scale_madgraph = 1000.0*4.890*3048.0/3.59644320000000000e+07;

    hsherpa->Sumw2();
    hsherpa->Scale(scale_sherpa);

    hmadgraph->Sumw2();
    hmadgraph->Scale(scale_madgraph);

    cout << "hmadgraph integral = " << hmadgraph->Integral() << endl;
    cout << "hsherpa integral = "   << hsherpa->Integral() << endl;;

    // get the ratio
    double chi2 = 0;
    int realbin = 0;
    for(int i=1; i<= nbins; i++) {

        double nmc=hmadgraph->GetBinContent(i);
        double ndata=hsherpa->GetBinContent(i);
        double nmcerr=hmadgraph->GetBinError(i);
        double ndataerr=hsherpa->GetBinError(i);

        if(nmc<0 || ndata<0)continue;

        if(nmcerr==0 && ndataerr==0)continue;

        if(nmc==0 && ndata==0)continue;

        double chi2ndef = (nmc-ndata)*(nmc-ndata)/
                          ( nmcerr*nmcerr+ ndataerr*ndataerr);
        chi2 += chi2ndef;
        realbin++;

        cout << "Bin " << i << " : " << ndata << ", " << nmc;
        cout << " " << chi2ndef << endl;


        // now calculate the ratio
        if(nmc==0 || nmcerr==0 || ndata==0 || ndataerr==0)continue;
        cout << "Bin " << i << " ratio = " << ndata/nmc << endl;
        hscale->SetBinContent(i,ndata/nmc);
        double err = 0;
        err=
            (ndata/nmc)*sqrt(pow(nmcerr/nmc,2)+pow(ndataerr/ndata,2));
        hscale->SetBinError(i,err);

    }


    hsherpa->GetXaxis()->SetRangeUser(xmin,xmax);

void plotRatioDoubleratio(TString varname_="", TString vartex_="")
{
  gStyle->SetTextSize(0.05);
  gStyle->SetTextFont(42);
  gStyle->SetPadRightMargin(0.05);
  gStyle->SetPadLeftMargin(0.15);
  gStyle->SetPadTopMargin(0.1);
  gStyle->SetPadBottomMargin(0.145);
  gStyle->SetTitleX(.0f);
  gStyle->SetOptStat(0);
  gStyle->SetMarkerStyle(20);
  gStyle->SetMarkerSize(0.8);

  varname=varname_;
  vartex=vartex_;
  
  TFile* infPP = new TFile(Form("outfDoubleratio/fPP_%s_DoubleRatio.root",varname.Data()));
  TH1D* hPP = (TH1D*)infPP->Get("hDoubleRatio");
  hPP->SetName(Form("hPP_%.0f",varname.Data()));
  TFile* infPbPb = new TFile(Form("outfDoubleratio/fPbPb_%s_DoubleRatio.root",varname.Data()));
  TH1D* hPbPb = (TH1D*)infPbPb->Get("hDoubleRatio");
  hPbPb->SetName(Form("hPbPb_%.0f",varname.Data()));
  TFile* infPPMB = new TFile(Form("outfDoubleratio/fPPMB_%s_DoubleRatio.root",varname.Data()));
  TH1D* hPPMB = (TH1D*)infPPMB->Get("hDoubleRatio");
  hPPMB->SetName(Form("hPPMB_%.0f",varname.Data()));
  TFile* infPbPbMB = new TFile(Form("outfDoubleratio/fPbPbMB_%s_DoubleRatio.root",varname.Data()));
  TH1D* hPbPbMB = (TH1D*)infPbPbMB->Get("hDoubleRatio");
  hPbPbMB->SetName(Form("hPbPbMB_%.0f",varname.Data()));

  TH1F* hRatio = (TH1F*)hPbPb->Clone("hRatio");
  hRatio->Divide(hPP);
  TH1F* hRatioMB = (TH1F*)hPbPbMB->Clone("hRatioMB");
  hRatioMB->Divide(hPPMB);

  Int_t varbins = hRatio->GetNbinsX();
  Float_t varstep = hRatio->GetBinWidth(1);
  Float_t varmin = hRatio->GetBinCenter(1)-0.5*varstep;
  Float_t varmax = hRatio->GetBinCenter(varbins)-0.5*varstep;
  Float_t aRatio[varbins],aRatioErr[varbins];
  Float_t aRatioMB[varbins],aRatioMBErr[varbins];
  Float_t aX[varbins],aZero[varbins];
  for(int i=0;i<varbins;i++)
    {
      aX[i] = hRatio->GetBinCenter(i+1)-0.5*varstep;
      aZero[i] = 0;
      aRatio[i] = hRatio->GetBinContent(i+1);
      aRatioErr[i] = hRatio->GetBinError(i+1);
      aRatioMB[i] = hRatioMB->GetBinContent(i+1);
      aRatioMBErr[i] = hRatioMB->GetBinError(i+1);
    }

  TGraphErrors* gRatio = new TGraphErrors(varbins,aX,aRatio,aZero,aRatioErr);
  TH2F* hemptyRatio = new TH2F("hemptyRatio","",20,varmin-0.5*varstep,varmax+0.5*varstep,10.,0.7,1.3);
  hemptyRatio->GetXaxis()->SetTitle(Form("%s",vartex.Data()));
  hemptyRatio->GetYaxis()->SetTitle("DoubleRatio^{PbPb} / DoubleRatio^{PP}");
  hemptyRatio->GetXaxis()->SetTitleOffset(1.);
  hemptyRatio->GetYaxis()->SetTitleOffset(.9);
  hemptyRatio->GetXaxis()->SetTitleSize(0.045);
  hemptyRatio->GetYaxis()->SetTitleSize(0.045);
  hemptyRatio->GetXaxis()->SetTitleFont(42);
  hemptyRatio->GetYaxis()->SetTitleFont(42);
  hemptyRatio->GetXaxis()->SetLabelFont(42);
  hemptyRatio->GetYaxis()->SetLabelFont(42);
  hemptyRatio->GetXaxis()->SetLabelSize(0.04);
  hemptyRatio->GetYaxis()->SetLabelSize(0.04);
  TCanvas* cRatio = new TCanvas("cRatio","",600,600);
  hemptyRatio->Draw();
  gRatio->Draw("psame");
  cRatio->SaveAs(Form("plotRatios/cRatio_%s.pdf",varname.Data()));

  TGraphErrors* gRatioMB = new TGraphErrors(varbins,aX,aRatioMB,aZero,aRatioMBErr);
  TH2F* hemptyRatioMB = new TH2F("hemptyRatioMB","",20,varmin-0.5*varstep,varmax+0.5*varstep,10.,0.7,1.3);
  hemptyRatioMB->GetXaxis()->SetTitle(Form("%s",vartex.Data()));
  hemptyRatioMB->GetYaxis()->SetTitle("DoubleRatio^{PbPb} / DoubleRatio^{PP}");
  hemptyRatioMB->GetXaxis()->SetTitleOffset(1.);
  hemptyRatioMB->GetYaxis()->SetTitleOffset(.9);
  hemptyRatioMB->GetXaxis()->SetTitleSize(0.045);
  hemptyRatioMB->GetYaxis()->SetTitleSize(0.045);
  hemptyRatioMB->GetXaxis()->SetTitleFont(42);
  hemptyRatioMB->GetYaxis()->SetTitleFont(42);
  hemptyRatioMB->GetXaxis()->SetLabelFont(42);
  hemptyRatioMB->GetYaxis()->SetLabelFont(42);
  hemptyRatioMB->GetXaxis()->SetLabelSize(0.04);
  hemptyRatioMB->GetYaxis()->SetLabelSize(0.04);
  TCanvas* cRatioMB = new TCanvas("cRatioMB","",600,600);
  hemptyRatioMB->Draw();
  gRatioMB->Draw("psame");
  cRatioMB->SaveAs(Form("plotRatios/cRatioMB_%s.pdf",varname.Data()));

}

//Compare the btag histo files histograms bin-by-bin
void CompareBtagEffFiles( TString ind = "" , TString inf = "" ) {

  const int ngraph = 9;       // # plots to compare
  const int nbins = 19;       // # bins to compare in each plot
  const int taggedgraphs = 6; // # plots which are simple counts only

  TFile myFiled(ind);
  TFile myFilef(inf);
  if ( myFiled.IsZombie() || myFilef.IsZombie() ) { cout << "Files are zombies!" << endl; return; }
 
  string names_d[ngraph] = {"h_bjet","h_cjet","h_ljet","h_btag","h_ctag","h_ltag","h_btageff","h_ctageff","h_ltageff"};
  string names_f[ngraph] = {"h_bjet","h_cjet","h_ljet","h_btag","h_ctag","h_ltag","h_btageff","h_ctageff","h_ltageff"};
  double ptbins[18] = {20, 30, 40, 50, 60, 70, 80, 100, 120, 160, 210, 260, 320, 400, 500, 600, 800, 99999}; 
 
  TH1D * histd;
  TH1F * histf;

  cout << endl << "Comparing histograms in the btag eff histo files bin-by-bin (inc under/over-flow). " << endl;
  cout << "-------------------------------------------------------------" << endl;
  for( int i = 0; i < ngraph; i++ ) {
    //This 'hack' absolutely needed. Don't ask...
    char  hist_named[200], hist_namef[200]; 
    sprintf(hist_named,"%s",names_d[i].c_str());
    sprintf(hist_namef,"%s",names_f[i].c_str());

    cout << endl << "Now processing " << names_d[i] << " and " << names_f[i] << " ... "  << endl << endl;
    myFiled.cd("");
    histd = (TH1D*) gDirectory->Get(hist_named); 
    myFilef.cd("");
    histf = (TH1F*) gDirectory->Get(hist_namef); 
    if (histd->IsZombie() ) { cout << "There is something wrong with " << names_d[i] << "d. Please fix: now exiting." << endl; return; } 
    if (histf->IsZombie() ) { cout << "There is something wrong with " << names_f[i] << "f. Please fix: now exiting." << endl; return; }
 
    if ( i < taggedgraphs )
      cout << "  | TH1D       : TH1F       : D/F " << endl;
    else
      cout << "  | TH1D              : TH1F              : D/F " << endl;

    double totd = 0, totf = 0;
    for( int j = 0; j < nbins; j++ ) {
      if ( i < taggedgraphs ) { 
        printf("%2i| %10.0f : %-10.0f : %5.4f\n",j,histd->GetBinContent(j),histf->GetBinContent(j),histd->GetBinContent(j)/histf->GetBinContent(j));
        totd +=  histd->GetBinContent(j); 
        totf +=  histf->GetBinContent(j); 
      }
      else
        printf("%2i| %15.15f : %15.15f : %5.4f\n",j,histd->GetBinContent(j),histf->GetBinContent(j),histd->GetBinContent(j)/histf->GetBinContent(j));
    }//go through the bins..

    if ( i < taggedgraphs ) {
      cout << "  " << endl;
      printf("SUM %10.0f : %-10.0f : %5.4f\n",totd, totf, totd/totf);
    }
    cout << "-------------------------------------------------------------" << endl;

  }//end for  

  myFiled.Close();
  myFilef.Close();
 
  cout << endl << " .. Done " << endl;
} // end of function

//.........这里部分代码省略.........
  TLatex * CMSLabel = new TLatex (0.18, 0.96, "#bf{CMS}");
  CMSLabel->SetNDC ();
  CMSLabel->SetTextAlign (10);
  CMSLabel->SetTextFont (42);
  CMSLabel->SetTextSize (_tsize);
  CMSLabel->Draw ("same") ;


  TLatex * _lumiLabel = new TLatex (0.95, 0.96, "19.4fb#lower[0.3]{^{-1}} (8 TeV)");
  _lumiLabel->SetNDC ();
  _lumiLabel->SetTextAlign (30);
  _lumiLabel->SetTextFont (42);
  _lumiLabel->SetTextSize (_tsize);
  _lumiLabel->Draw ("same") ;


  // Draw also ratio
  //----------------------------------------------------------------------------
  if (drawRatio) {

     pad2->cd();
	
    TH1F* ratio_pow       = xsValue_Powheg->Clone("ratio");
    TH1F* ratio_mad       = xsValue_Madgraph->Clone("ratio");
    TH1F* ratio_mcnlo     = xsValue_MCnlo->Clone("ratio");
    TH1F* hratio_pow       = xsValue_Powheg->Clone("ratio");
    TH1F* hratio_mad       = xsValue_Madgraph->Clone("ratio");
    TH1F* hratio_mcnlo     = xsValue_MCnlo->Clone("ratio");
    TH1F* ratioErr        = xsValue->Clone("ratio");
    

    for (UInt_t ibin=1; ibin<=ratio->GetNbinsX(); ibin++) {

      Double_t powValue = xsValue_Powheg->GetBinContent(ibin);
      Double_t powError = xsValue_Powheg->GetBinError  (ibin);
      
      Double_t madValue = xsValue_Madgraph->GetBinContent(ibin);
      Double_t madError = xsValue_Madgraph->GetBinError  (ibin);
      
      Double_t mcnloValue = xsValue_MCnlo->GetBinContent(ibin);
      Double_t mcnloError = xsValue_MCnlo->GetBinError  (ibin);
      
      Double_t dataValue = xsValue->GetBinContent(ibin);
      Double_t statError = xsValue->GetBinError  (ibin);
      Double_t systError = systHisto->GetBinError(ibin);
      
      Double_t dataError = systError;
      
      Double_t ratioValue_pow           = (powValue > 0) ? powValue/dataValue : 0.0;
      Double_t ratioError_pow           = (powValue > 0) ? powError / dataValue : 0.0;
      
      Double_t ratioValue_mad           = (madValue > 0) ? madValue/dataValue : 0.0;
      Double_t ratioError_mad           = (madValue > 0) ? madError/dataValue : 0.0;
      
      Double_t ratioValue_mcnlo         = (mcnloValue > 0) ? mcnloValue/dataValue : 0.0;
      Double_t ratioError_mcnlo         = (mcnloValue > 0) ? mcnloError/dataValue : 0.0;
      
      Double_t uncertaintyError         = (dataValue > 0) ? dataError/dataValue : 0.0;
      

      //dataError/dataValue 
      ratio_pow->SetBinContent(ibin, ratioValue_pow);
      hratio_pow->SetBinContent(ibin, ratioValue_pow);
      hratio_pow->SetBinError  (ibin, ratioError_pow);
      
      ratio_mad->SetBinContent(ibin, ratioValue_mad);

void cetaflatHFM12(int nIterN=1, double Ethr1=10, double Ethr2=150) { // for HFM, L and S separately
  // nIterN - number of iterations
  // Ethr1 and Ethr2 - E thresholds within which E is estimated
  
  gStyle->SetOptLogz(0);
  gStyle->SetMarkerSize(0.7);
  gStyle->SetMarkerStyle(20);
  gStyle->SetPadGridX(0);
  gStyle->SetPadGridY(0);
  gStyle->SetTitleOffset(1.7,"Y");
  gStyle->SetTitleOffset(0.9,"X");
  //gStyle->SetPadRightMargin(0.12);
  gStyle->SetPadRightMargin(0.03);
  gStyle->SetPadLeftMargin(0.18);
  //gStyle->SetNdivisions(516);
  gStyle->SetStatH(0.025);
  gStyle->SetStatW(0.3);
  gStyle->SetTitleW(0.4);
  gStyle->SetTitleX(0.28);
  gStyle->SetOptStat(0);
  gROOT->ForceStyle();

  char ctit[245],ftit[245];
  float etaBounds[14] = {2.853,2.964,3.139,3.314,3.489,3.664,3.839,4.013,4.191,4.363,4.538,4.716,4.889,5.205};

  // ------Histos input: spectra of all channels-----------------------------------

  //sprintf(ftit,"%s","phi43val2012A");
  //sprintf(ftit,"%s","phi2012A_May");
  //sprintf(ftit,"%s","phiSym524_2012AB");
  //sprintf(ftit,"%s","phiSym524newGain_2012AB");
  //sprintf(ftit,"%s","phiSym524newGain_2012ABC");
  //sprintf(ftit,"%s","phisymNewCond2012Cval");
  //sprintf(ftit,"%s","phisymOldCond2012Cval");
  //sprintf(ftit,"%s","phiSym533Gain507_2012D");
  //  sprintf(ftit,"%s","phiSym533Corr45Gain507_2012D");  //  histo root file processed

  //sprintf(ctit,"/home/vodib/beam12/intercal/%s.root",ftit);
  sprintf(ctit,"hcal.root",ftit);    
TFile *fila = new TFile (ctit);
  cout<<"File= "<<ctit<<endl;

  TH1F *hcounter =   new TH1F(*((TH1F*)fila->Get("phaseHF/hcounter")));
  cout<<"Stat= "<<hcounter->GetBinContent(2)<<endl;
  cout<<"E within: "<<Ethr1<<" - "<<Ethr2<<endl;

  TH2F* hLmapP = new TH2F("hLmapP","E L HFM;i#eta;i#phi",13,-41.5,-28.5,36,0,72);
  TH2F* hSmapP = new TH2F("hSmapP","E S HFM;i#eta;i#phi",13,-41.5,-28.5,36,0,72);
  TH2F* hLmapM0 = new TH2F("hLmapM0","E0 L HFM;i#eta;i#phi",13,-41.5,-28.5,36,0,72);
  TH2F* hSmapM0 = new TH2F("hSmapM0","E0 S HFM;i#eta;i#phi",13,-41.5,-28.5,36,0,72);
  TH2F* hLmapPc = new TH2F("hLmapPc","corr L HFM;i#eta;i#phi",13,-41.5,-28.5,36,0,72);
  TH2F* hSmapPc = new TH2F("hSmapPc","corr S HFM;i#eta;i#phi",13,-41.5,-28.5,36,0,72);
  hLmapPc->Sumw2(); hSmapPc->Sumw2();
  //TH1F *hLcorr1D = new TH1F("hLcorr1D","Corr L",300,0.5,2);
  //TH1F *hScorr1D = new TH1F("hScorr1D","Corr S",300,0.5,2);
  TH1F *hLcorr1D = new TH1F("hLcorr1D","Corr L",180,0.7,1.5);
  TH1F *hScorr1D = new TH1F("hScorr1D","Corr S",180,0.7,1.5);
  TH1F *hLdatP[13][36], *hSdatP[13][36], *hLdatPx[13][36], *hSdatPx[13][36];
  for (int ii=0;ii<13;ii++) for (int jj=0;jj<36;jj++) {
    sprintf(ctit,"hL%d_%d",-29-ii,2*jj+1);
    hLdatP[ii][jj] = new TH1F(ctit,ctit,8000,0,250);
    sprintf(ctit,"hS%d_%d",-29-ii,2*jj+1);
    hSdatP[ii][jj] = new TH1F(ctit,ctit,8000,0,250);
  }
  TH1F *htL = new TH1F("htL","htL",20000,0,7e8/3.);
  TH1F *htS = new TH1F("htS","htS",20000,0,5e8/3.);
  //TH1F *htL = new TH1F("htL","htL",20000,0,4e8/40);
  //TH1F *htS = new TH1F("htS","htS",20000,0,2e8/40);
  TH1F *hLdatPx[13][36], *hSdatPx[13][36];

  TCanvas *cLx[200],*cSx[200];
  TSpline5 *ttL,*ttS;

  Double_t x,y,rPL,rPS,drPL,drPS,mLE,mSE,ermean,rms;
  Double_t xxL[1000],yyL[1000];
  Double_t xxS[1000],yyS[1000];
  Int_t nELP, nESP, nIter=0;
  Double_t mcorrL,scorrL,mcorrS,scorrS,erLP,erSP,rLP,drLP,rSP,drSP,corrL,corrS,dcorrL,dcorrS;
  double mLEphi[13],mSEphi[13],dmLEphi[13],dmSEphi[13];

  TCanvas *ccxx = new TCanvas("ccxx","ccxx",100,300,900,500);
  ccxx->Divide(2,1);

  // loop over channels, 
  for (int ii=0;ii<13;ii++) {  // loop over HFM ieta (rings)
  //for (int ii=1;ii<2;ii++) {
    int ieta=-ii-29;

    mLE=mSE=0;   // ------------------for initial condition
    int nmLE=0, nmSE=0;
    htL->Reset(); htS->Reset();
    for (int ll=1;ll<=72;ll+=2) { // loop over channels within HFM ieta (ring) to get data
      int iphi=ll;
      if (abs(ieta)>39 && (iphi-1)%4==0) continue;
      //if (ieta==29 && iphi==67) continue;
      hSmapPc->SetBinContent(13-ii,ll/2+1,1);  //  initial correction set to 1
      hLmapPc->SetBinContent(13-ii,ll/2+1,1);
      hSmapPc->SetBinError(13-ii,ll/2+1,1.e-6);  //  non-zero err required
      hLmapPc->SetBinError(13-ii,ll/2+1,1.e-6);
      sprintf(ctit,"phaseHF/espec/E_%d_%d_1",ieta,iphi);
//.........这里部分代码省略.........