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

// ============================================================================
StatusCode Gaudi::Utils::Histos::fromXml
( TH1F& result , const std::string& input )
{
  //
  result.Reset() ;                                 // RESET old histogram
  //
  _Xml<TH1F> _xml ;
  std::auto_ptr<TH1F> histo =  _xml ( input ) ;
  if ( 0 == histo.get() ) { return StatusCode::FAILURE ; }        // RETURN
  //
  result.Reset() ;
  histo->Copy ( result ) ;
  //
  return StatusCode::SUCCESS ;
}

double NewChi2Func(const double *xx ){
 const Double_t scale = xx[0];
 fileInMC->cd();
 TTree* MyTreeMC = (TTree*) fileInMC->Get(treeNameMC.c_str());
 outFile->cd();
 TString nameDATA = Form("hDATA_%d_%d_%.5f",Data_or_MC,nIter,ScaleTrue);
 TH1F* hDATA = (TH1F*) outFile->Get(nameDATA);
 
 TString NameMC = Form("hMC_Chi2_%.5f",scale);
 TH1F* hMC;
 if (!gROOT->FindObject(NameMC.Data())){
  hMC = new TH1F(NameMC,NameMC,numBINS,minBINS,maxBINS);
  hMC->Reset();
  
  MyTreeMC->SetEntryList(0); 
  MyTreeMC->Draw(">> myListMCTot",(AdditionalCut + Form(" && (ET * (1+(%f)))>%f",scale,minET)).Data(),"entrylist");
  TEntryList *mylistMCTot = (TEntryList*)gDirectory->Get("myListMCTot");
  MyTreeMC->SetEntryList(mylistMCTot);
  
  TString DrawMC = Form("(%s * (1+(%f)))>>%s",variableName.c_str(),scale,NameMC.Data());
  MyTreeMC->Draw(DrawMC);
//   MyTreeMC->Draw(DrawMC,(AdditionalCut+Form("&& (ET * (1+(%f)))>%f",scale,minET)).Data());
  hMC->Sumw2();
  hMC->Scale(hDATA->GetEffectiveEntries()/hMC->GetEffectiveEntries());
  outFile->cd();
  hMC->Write();
 } 
 else {
//   std::cerr << " NewChi2 old " << NameMC.Data() << std::endl;
  hMC = (TH1F*) outFile->Get(NameMC.Data());
 }
 
 double result = myChi2(hMC,hDATA);
 return result;
}

TH1F* Prediction::CalcPrediction(TH2F* prediction_raw) {

   TH1F* prediction = new TH1F();
   prediction = (TH1F*) prediction_raw->ProjectionX();
   prediction->Reset();

   for (int i = 0; i <= prediction_raw->GetXaxis()->GetNbins() + 1; ++i) {
   
      TH1F h = *((TH1F*) prediction_raw->ProjectionY("py", i, i));
      
      double summ = 0;
      double sumv = 0;
      int N = 0;
      
      //// Calculate mean
      for (int j = 1; j <= (int) h.GetNbinsX(); ++j) {
         summ += h.GetBinContent(j);
         ++N;
      }
      double mean = summ / N;
      
      //// Calculated variance
      for (int j = 1; j <= (int) h.GetNbinsX(); ++j) {
         sumv += pow(mean - h.GetBinContent(j), 2);
      }
      double variance = sqrt(sumv / N);
      
      prediction->SetBinContent(i, mean);
      prediction->SetBinError(i, variance);

   }
   
   return prediction;
   
}

TH1F RestoreBinning(TH1F const& src, TH1F const& ref) {
  TH1F res = ref;
  res.Reset();
  for (int x = 1; x <= res.GetNbinsX(); ++x) {
    res.SetBinContent(x, src.GetBinContent(x));
  }
  return res;
}

TH1F* eff_bg(TH1F* h1, TH1F* h2, TH1F* h3, TH1F* h4, const char* name="eff"){

  // first, verify that all histograms have same binning
  // nx is the number of visible bins
  // nxtot = nx+2 includes underflow and overflow
  Int_t nx = h1->GetNbinsX();
  Int_t nxtot = nx + 2;
  if (h2->GetNbinsX() != nx) {
    //    cout << "Histograms must have same number of bins" << endl;
    return 0;
  }
  if (h3->GetNbinsX() != nx) {
    //    cout << "Histograms must have same number of bins" << endl;
    return 0;
  }
  if (h3->GetNbinsX() != nx) {
    //    cout << "Histograms must have same number of bins" << endl;
    return 0;
  }

  // get the new histogram
  TH1F* temp = (TH1F*) h1->Clone(name);
  temp->SetTitle(name);
  temp->Reset();
  temp->Sumw2();

  // Loop over bins, calculate efficiency and error, put it in histogram
  for (Int_t i=0; i<nxtot; i++) {
    Double_t x1 = h1->GetBinContent(i);
    Double_t x2 = h2->GetBinContent(i);
    Double_t x3 = h3->GetBinContent(i);
    Double_t x4 = h4->GetBinContent(i);
    Double_t denom = x1 - x3;
    Double_t eff;
    if (denom == 0.) {
      eff = 0;
    } else {
      eff   = (x2-x4)/denom;
    }
    Double_t failSig = x1 - x2;
    Double_t failBg  = x3 - x4;
    Double_t blah    = (1-eff)*(1-eff)*(x2+x4) + eff*eff*(failSig+failBg);
    if (blah <= 0.) blah=0.0;
    Double_t err;
    if (denom == 0) {
      err = 0.;
    } else {
      err = sqrt(blah)/denom;
    }
    temp->SetBinContent(i,eff);
    temp->SetBinError(i,err);
  }

  // Done
  return temp;
}

TH1F * MakeOnePSE(TH1F * sig, int nsig, TH1F * back, int nback){

  //threw one data PSE  and fill it
  TH1F * data  = (TH1F*) sig ->Clone("data");
  data->Reset();
  data->FillRandom(sig,nsig);
  data->FillRandom(back,nback);

  return data;

}//MakeOnePSE

TH1F * getErrHist(TH1F * h, bool rel = true)
{
	TH1F * res = (TH1F *) h->Clone("errhist");
	res->Reset();
	
	for(int j = 1; j <= res->GetNbinsX(); j++)
		if(rel) res->SetBinContent(j,h->GetBinError(j)/h->GetBinContent(j));
		else res->SetBinContent(j,h->GetBinError(j));

	return res;
}

TH1F * sqSum(TH1F * h1, TH1F * h2)
{
	TH1F * res = (TH1F *) h2->Clone("sum");
	res->Reset();
	
	for(int j = 1; j <= res->GetNbinsX(); j++)
	{
		if(h1) res->SetBinContent(j,TMath::Sqrt( TMath::Power(h1->GetBinContent(j),2) + TMath::Power(h2->GetBinContent(j),2)));
		else res->SetBinContent(j,h2->GetBinContent(j));
	}
	
	return res;
}

double Chi2F(const double *xx ){
 const Double_t scale = xx[0];
 fileInMC->cd();
 TTree* MyTreeMC = (TTree*) fileInMC->Get(treeNameMC.c_str());
 outFile->cd();
 TString nameDATA = Form("hDATA_%d_%d_%.5f",Data_or_MC,nIter,ScaleTrue);
 TH1F* hDATA = (TH1F*) outFile->Get(nameDATA);
 TString NameMC = Form("hMC_Chi2_%.5f",scale);
 
 //std::cerr << " NameMC = "  << NameMC.Data() << " => " << scale << std::endl;
 
 TH1F* hMC;
 if (!gROOT->FindObject(NameMC.Data())){
  hMC = new TH1F(NameMC,NameMC,numBINS,minBINS,maxBINS);
  hMC->Reset();

//  std::cerr << " AdditionalCut.Data() = " << AdditionalCut.Data() << std::endl;
  
  MyTreeMC->SetEntryList(0); 
  MyTreeMC->Draw(">> myListMCTot",(AdditionalCut + Form(" && (ET * (1+(%f)))>%f",scale,minET)).Data(),"entrylist");
  TEntryList *mylistMCTot = (TEntryList*)gDirectory->Get("myListMCTot");
  MyTreeMC->SetEntryList(mylistMCTot);
  
  TString DrawMC = Form("(%s * (1+(%f)))>>%s",variableName.c_str(),scale,NameMC.Data());
//  std::cerr << " DrawMC = " << DrawMC.Data() << std::endl;
  MyTreeMC->Draw(DrawMC);
//   MyTreeMC->Draw(DrawMC,(AdditionalCut+Form("&& (ET * (1+(%f)))>%f",scale,minET)).Data());
  hMC->Sumw2();
  hMC->Scale(hDATA->GetEffectiveEntries()/hMC->GetEffectiveEntries());
  outFile->cd();
  hMC->Write();
 }
 else {
//   std::cerr << " KM old "  << NameMC.Data() << std::endl;
  hMC = (TH1F*) outFile->Get(NameMC.Data());
 }
 
 outFile->cd();
//  hDATA.Write();
//  hMC.Write();
 
 double result = hMC->KolmogorovTest(hDATA,"M");
// double result = - hMC->KolmogorovTest(hDATA,"X");
//  double result = hMC->Chi2Test(&hDATA,"CHI2/NDF");
//=========> E' QUESTO! ==> double result = hMC->Chi2Test(hDATA,"CHI2");
 //double result = - hMC.Chi2Test(&hDATA,""); ///==== http://root.cern.ch/root/html/TH1.html#TH1:Chi2Test
 return result;
}

void make(TDirectory & out, TObject * o) {
  TDirectory * dir;
  TH1F * th1f;
  TH1D * th1d;
  TH2F * th2f;
  TH2D * th2d;
  out.cd();
  if((dir = dynamic_cast<TDirectory*>(o)) != 0) {
    TDirectory * outDir = out.mkdir(dir->GetName(), dir->GetTitle());
    TIter next(dir->GetListOfKeys());
    TKey *key;
    while( (key = dynamic_cast<TKey*>(next())) ) {
      string className(key->GetClassName());
      string name(key->GetName());
      TObject * obj = dir->Get(name.c_str());
      if(obj == 0) {
	cerr <<"error: key " << name << " not found in directory " << dir->GetName() << endl;
	exit(-1);
      }
      make(*outDir, obj);
    }
  } else if((th1f = dynamic_cast<TH1F*>(o)) != 0) {
    TH1F *h = (TH1F*) th1f->Clone();
    h->Reset();
    h->Sumw2();
    h->SetDirectory(&out);
  } else if((th1d = dynamic_cast<TH1D*>(o)) != 0) {
    TH1D *h = (TH1D*) th1d->Clone();
    h->Reset();
    h->Sumw2();
    h->SetDirectory(&out);
  } else if((th2f = dynamic_cast<TH2F*>(o)) != 0) {
    TH2F *h = (TH2F*) th2f->Clone();
    h->Reset();   
    h->Sumw2();
    h->SetDirectory(&out);
  } else if((th2d = dynamic_cast<TH2D*>(o)) != 0) {
    TH2D *h = (TH2D*) th2d->Clone();
    h->Reset();   
    h->Sumw2();
    h->SetDirectory(&out);
  }
}

// Method by pointer
TH1F* eff(TH1F* h1, TH1F* h2, const char* name="eff"){

  // first, verify that all histograms have same binning
  // nx is the number of visible bins
  // nxtot = nx+2 includes underflow and overflow
  Int_t nx = h1->GetNbinsX();
  if (h2->GetNbinsX() != nx) {
    cout << "Histograms must have same number of bins" << endl;
    return 0;
  }

  // get the new histogram
  TH1F* temp = (TH1F*) h1->Clone(name);
  temp->SetTitle(name);
  temp->Reset();
  temp->Sumw2();

  // Do the calculation
  temp->Divide(h2,h1,1.,1.,"B");

  // Done
  return temp;
}

void efficiency( TString denfile="60-120.root",
                 TString numfile="orig.root",
                 TString outfile="res.root" ) {

  gROOT ->Reset();
  
  TFile * nume = new TFile(numfile);
  TFile * deno = new TFile(denfile);
  
  nume->cd("/");
  deno->cd("/");

  TString theName = "d02-x01-y01";

  const TH1F *myNume;
  nume->GetObject(theName,myNume);
  myNume->Sumw2();
  myNume->Print("all");

  const TH1F *myDeno;
  deno->GetObject(theName,myDeno);
  myDeno->Sumw2();
  myDeno->Print("all");

  TH1F *Acce = (TH1F*)myNume->Clone("efficiency");
  Acce->Reset();
  Double_t f1(1.),f2(1.);
  Option_t* opt("b");
  Acce->Divide(myNume,myDeno,f1,f2,opt);
  Acce->Print("all");

  TFile * out = new TFile(outfile,"NEW");
  Acce->Write();
  out->Close();

}

void cetaflatHFP12(int nIterN=1, double Ethr1=10, double Ethr2=150) {
  
  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");

  sprintf(ctit,"/home/vodib/beam12/intercal/%s.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 HFP;i#eta;i#phi",13,28.5,41.5,36,0,72);
  TH2F* hSmapP = new TH2F("hSmapP","E S HFP;i#eta;i#phi",13,28.5,41.5,36,0,72);
  TH2F* hLmapP0 = new TH2F("hLmapP0","E0 L HFP;i#eta;i#phi",13,28.5,41.5,36,0,72);
  TH2F* hSmapP0 = new TH2F("hSmapP0","E0 S HFP;i#eta;i#phi",13,28.5,41.5,36,0,72);
  TH2F* hLmapPc = new TH2F("hLmapPc","corr L HFP;i#eta;i#phi",13,28.5,41.5,36,0,72);
  TH2F* hSmapPc = new TH2F("hSmapPc","corr S HFP;i#eta;i#phi",13,28.5,41.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",ii+29,2*jj+1);
    hLdatP[ii][jj] = new TH1F(ctit,ctit,8000,0,250);
    sprintf(ctit,"hS%d_%d",ii+29,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);

  for (int ii=0;ii<13;ii++) {
  //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) {
      int iphi=ll;
      if (abs(ieta)>39 && (iphi-1)%4==0) continue;
      hSmapPc->SetBinContent(ii+1,ll/2+1,1);
      hLmapPc->SetBinContent(ii+1,ll/2+1,1);
      hSmapPc->SetBinError(ii+1,ll/2+1,1.e-6);
      hLmapPc->SetBinError(ii+1,ll/2+1,1.e-6);
      sprintf(ctit,"phaseHF/espec/E_+%d_%d_1",ieta,iphi);
      hLdatPx[ii][ll/2]  =   new TH1F(*((TH1F*)fila->Get(ctit)));
      hLdatPx[ii][ll/2]->SetAxisRange(Ethr1,Ethr2);
      rLP = hLdatPx[ii][ll/2]->Integral()*hLdatPx[ii][ll/2]->GetMean();
      hLmapP0->SetBinContent(ii+1,ll/2+1,rLP);
      sprintf(ctit,"phaseHF/espec/E_+%d_%d_2",ieta,iphi);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
   if(jetPtHatHi.size()!= nSize_jet){cout << "error 4"<< endl; return;}


   // last fill mix trees
   fillTrees(mixFile,mixTree,treeName);
   const unsigned int nSize_mix = mixFile.size();
   if(mixTree.size()!= nSize_mix){cout << "error 1"<< endl; return;}
   if(mixWeight.size()!= nSize_mix){cout << "error 2"<< endl; return;}
   if(mixPtHatLo.size()!= nSize_mix){cout << "error 3"<< endl; return;}
   if(mixPtHatHi.size()!= nSize_mix){cout << "error 4"<< endl; return;}
   
   TH1F *hTemplate = new TH1F("hTemplate","",70,-5,2);
   TH1F* hEcalIsoPho = (TH1F*)hTemplate->Clone();
   std::string histoName = outputName + "Pho";
   hEcalIsoPho->SetName(histoName.data());

   TH1F* hEcalIsoJet = (TH1F*)hTemplate->Clone();
   histoName = outputName + "Jet";
   hEcalIsoJet->SetName(histoName.data());

   TH1F* hEcalIsoMixSig = (TH1F*)hTemplate->Clone();
   histoName = outputName + "MixSig";
   hEcalIsoMixSig->SetName(histoName.data());

   TH1F* hEcalIsoMixBkg = (TH1F*)hTemplate->Clone();
   histoName = outputName + "MixBkg";
   hEcalIsoMixBkg->SetName(histoName.data());


   TH1F* hEcalIsoMixData = (TH1F*)hTemplate->Clone();
   histoName = outputName + "MixData";
   hEcalIsoMixData->SetName(histoName.data());
   
   cout << "making histograms from photon+jet MC samples" << endl;
   TCut allCut = basicCut + sigCut;
//    makePlot(phoTree,phoWeight,Form("%s",var.data()),hardScatterCut,hEcalIsoPho,normalize);
   makePlot(phoTree,phoWeight,phoPtHatLo,phoPtHatHi,Form("%s",var.data()),allCut,hEcalIsoPho,normalize);

   cout << "making histograms from dijet MC samples" << endl;
   allCut = basicCut + decayCut;
   makePlot(jetTree,jetWeight,jetPtHatLo,jetPtHatHi,Form("%s",var.data()),allCut,hEcalIsoJet,normalize);

   cout << "making histograms from mixed MC signal samples" << endl;     
   allCut = basicCut + sigCut;
   makePlot(mixTree,mixWeight,mixPtHatLo,mixPtHatHi,Form("%s",var.data()),allCut,hEcalIsoMixSig,normalize);

   cout << "making histograms from mixed MC background samples" << endl;     
   allCut = basicCut + bkgCut;
   makePlot(mixTree,mixWeight,mixPtHatLo,mixPtHatHi,Form("%s",var.data()),allCut,hEcalIsoMixBkg,normalize);
   
   hEcalIsoPho->SetMarkerColor(2);
   hEcalIsoPho->SetLineColor(2);
   hEcalIsoJet->SetMarkerColor(1);
   hEcalIsoJet->SetLineColor(1);
   hEcalIsoMixSig->SetMarkerColor(5);
   hEcalIsoMixSig->SetLineColor(5);
   hEcalIsoMixBkg->SetMarkerColor(4);
   hEcalIsoMixBkg->SetLineColor(4);


   cout << "hEcalIsoPho->Integral()  = " << hEcalIsoPho->Integral() << endl;
   cout << "hEcalIsoJet->Integral()  = " << hEcalIsoJet->Integral() << endl;
   cout << "hEcalIsoMixSig->Integral()  = " << hEcalIsoMixSig->Integral() << endl;
   cout << "hEcalIsoMixBkg->Integral()  = " << hEcalIsoMixBkg->Integral() << endl;

   hEcalIsoMixData->Reset();
   hEcalIsoMixData->Sumw2();
   hEcalIsoMixData->Add(hEcalIsoMixSig,hEcalIsoMixBkg,1.0,1.0);
   cout << "hEcalIsoMixData->Integral()  = " << hEcalIsoMixData->Integral() << endl;
   

   hEcalIsoPho->SetXTitle(var.data());
   hEcalIsoPho->Draw("hist");
   hEcalIsoJet->Draw("histesame");
   hEcalIsoMixSig->Draw("histesame");
   hEcalIsoMixBkg->Draw("histesame");
   TLegend* leg = new TLegend(0.5,0.6,0.7,0.9);
   leg->SetFillColor(0);
   leg->SetFillStyle(0);
   leg->SetTextSize(0.03);
   leg->SetBorderSize(0);
   leg->AddEntry(hEcalIsoPho,"#gamma+jet MC");
   leg->AddEntry(hEcalIsoJet,"Dijet MC");
   leg->AddEntry(hEcalIsoMixSig,"Mixed MC: signal");
   leg->AddEntry(hEcalIsoMixBkg,"Mixed MC: background");
   leg->Draw("same");

   // dump histogram to a root file
   std::string histoFile = outputName + "_histo.root";

   TFile* outFile = new TFile(histoFile.data(),"recreate");
   hEcalIsoPho->Write();
   hEcalIsoJet->Write();
   hEcalIsoMixSig->Write();
   hEcalIsoMixBkg->Write();
   hEcalIsoMixData->Write();
   outFile->Close();


}

// Function to draw EE, ME, MM channel and all channel merged plot
void cutStepPlots(const char* cutStep, const char* histNameTitle[2],
                  double minY, double maxY, bool doLogY)
{
 const char* histName = histNameTitle[0];
 const char* histTitle = histNameTitle[1];

  //cout << histName << " : " << histTitle << endl;

  TH1F* hDataEE = (TH1F*)fEE->Get(Form("%s/hData_%s_%s", cutStep, cutStep, histName));
  TH1F* hDataME = (TH1F*)fME->Get(Form("%s/hData_%s_%s", cutStep, cutStep, histName));
  TH1F* hDataMM = (TH1F*)fMM->Get(Form("%s/hData_%s_%s", cutStep, cutStep, histName));

  if ( !hDataEE ) { cout << Form("%s/hData_%s_%s", cutStep, cutStep, histName) << " for EE " << "\n"; return; }
  if ( !hDataME ) { cout << Form("%s/hData_%s_%s", cutStep, cutStep, histName) << " for ME " << "\n"; return; }
  if ( !hDataMM ) { cout << Form("%s/hData_%s_%s", cutStep, cutStep, histName) << " for MM " << "\n"; return; }


  TH1F* hDataLL = (TH1F*)hDataEE->Clone(Form("hData_%s_%s", cutStep, histName));
  hDataLL->Reset();
  hDataLL->Add(hDataEE);
  hDataLL->Add(hDataME);
  hDataLL->Add(hDataMM);

  THStack* hStackEE = new THStack(TString("hEE_")+cutStep+"_"+histName, histTitle);
  THStack* hStackME = new THStack(TString("hME_")+cutStep+"_"+histName, histTitle);
  THStack* hStackMM = new THStack(TString("hMM_")+cutStep+"_"+histName, histTitle);
  THStack* hStackLL = new THStack(TString("hLL_")+cutStep+"_"+histName, histTitle);

  TH1F* hSigEE = (TH1F*)fEE->Get(Form("%s/%s_%s_%s", cutStep, sigNames[0], cutStep, histName));
  TH1F* hSigME = (TH1F*)fME->Get(Form("%s/%s_%s_%s", cutStep, sigNames[0], cutStep, histName));
  TH1F* hSigMM = (TH1F*)fMM->Get(Form("%s/%s_%s_%s", cutStep, sigNames[0], cutStep, histName));

  hSigEE->SetFillColor(color_sig);
  hSigME->SetFillColor(color_sig);
  hSigMM->SetFillColor(color_sig);

  hSigEE->SetFillStyle(style_sig);
  hSigME->SetFillStyle(style_sig);
  hSigMM->SetFillStyle(style_sig);

  //hSigEE->SetLineColor(color_sig);
  //hSigME->SetLineColor(color_sig);
  //hSigMM->SetLineColor(color_sig);




  if ( !hSigEE || !hSigME || !hSigMM ) { cout << "No signal hist for " << histName << "\n"; return; }

  TH1F* hSigLL = (TH1F*)hSigEE->Clone(Form("%s_%s_%s", sigNames[0], cutStep, histName));
  hSigLL->Reset();
  hSigLL->Add(hSigEE);
  hSigLL->Add(hSigME);
  hSigLL->Add(hSigMM);

  if( stackSig ){
    hStackEE->Add(hSigEE);
    hStackME->Add(hSigME);
    hStackMM->Add(hSigMM);
    hStackLL->Add(hSigLL);
  }

  // Build legends
  TLegend* legEE = buildLegend();
  TLegend* legME = buildLegend();
  TLegend* legMM = buildLegend();
  TLegend* legLL = buildLegend();

  if ( hDataEE->GetEntries() > 0 ) legEE->AddEntry(hDataEE, "Data", "p");
  if ( hDataME->GetEntries() > 0 ) legME->AddEntry(hDataME, "Data", "p");
  if ( hDataMM->GetEntries() > 0 ) legMM->AddEntry(hDataMM, "Data", "p");
  if ( hDataLL->GetEntries() > 0 ) legLL->AddEntry(hDataLL, "Data", "p");

  TH1F* hEEs[nBkg];
  TH1F* hMEs[nBkg];
  TH1F* hMMs[nBkg];
  TH1F* hLLs[nBkg];

  for ( int i=0; i<nBkg; ++i )
  {
    TH1F* hEE = (TH1F*)fEE->Get(Form("%s/%s_%s_%s", cutStep, bkgNames[i], cutStep, histName));
    TH1F* hME = (TH1F*)fME->Get(Form("%s/%s_%s_%s", cutStep, bkgNames[i], cutStep, histName));
    TH1F* hMM = (TH1F*)fMM->Get(Form("%s/%s_%s_%s", cutStep, bkgNames[i], cutStep, histName));

    hEE->SetFillColor(color[i]);
    hME->SetFillColor(color[i]);
    hMM->SetFillColor(color[i]);

    hEE->SetFillStyle(style[i]);
    hME->SetFillStyle(style[i]);
    hMM->SetFillStyle(style[i]);

   // hEE->SetLineColor(color[i]);
   // hME->SetLineColor(color[i]);
   // hMM->SetLineColor(color[i]);


    if ( !hEE || !hME || !hMM ) { cout << "No bkg hist " << bkgNames[i] << endl; continue; }

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