C++ TH1D::Clone方法代码示例

TH1D* plot (std::string num, std::string den,std::vector<TFile*>& tfiles,std::vector<double>& weights) {

    double weightTot = 0.0;

    for ( int i = 0; i<weights.size(); i++) {
        weightTot+=weights[i];
    }

    std::string name = num;
    name+="NEW";
    TH1D* hNum = (TH1D*)tfiles[0]->Get(num.c_str());
    TH1D* HNum = (TH1D*)hNum->Clone(name.c_str());
    HNum->Sumw2();



    name=den;
    name+="NEW";
    TH1D* hDen = (TH1D*)tfiles[0]->Get(den.c_str());
    TH1D* HDen = (TH1D*)hDen->Clone(name.c_str());
    HDen->Sumw2();

    for (int i=1; i<tfiles.size(); i++) {
        TH1D* htempNum = (TH1D*)tfiles[i]->Get(num.c_str());
        TH1D* htempDen = (TH1D*)tfiles[i]->Get(den.c_str());
        HNum->Add(htempNum,weights[i]/weightTot);
        HDen->Add(htempDen,weights[i]/weightTot);
    }
    name = num;
    name+="Div";
    TH1D* hDiv = (TH1D*)HNum->Clone(name.c_str());
    //hDiv->Sumw2();
    hDiv->Divide(HDen);
    return hDiv;
}

void MCtoSTATerrors() {
  TLatex *tplus = labelLatex(0.1364943,0.8114407,"MC sample");

  TFile *file0 = TFile::Open("results/StandardCutflow/allmc.root");

  TCanvas *c1 = new TCanvas();
  //TCanvas *c2 = new TCanvas();
  TH1D * LPplus = (TH1D*)file0->Get("RECO_PolPlots_50toinf/RECO_ICVarPFPlus");
  TH1D * LPminus = (TH1D*)file0->Get("RECO_PolPlots_50toinf/RECO_ICVarPFMinus");

  LPplus->Rebin(20);
  LPminus->Rebin(20);

  LPplus->Scale(0.3);
  LPminus->Scale(0.3);

  TH1D * newLPplus = (TH1D*)LPplus->Clone();
  TH1D * newLPminus = (TH1D*)LPminus->Clone();

  for(unsigned int i=1; i < LPplus->GetXaxis()->GetNbins(); i++) {
    newLPplus->SetBinError(i, TMath::Sqrt(LPplus->GetBinContent(i)));
    newLPminus->SetBinError(i, TMath::Sqrt(LPminus->GetBinContent(i)));
  }

  c1->cd();
  newLPplus->GetYaxis()->SetTitle("Events / 300 nb^{-1}");
  //newLPplus->GetXaxis()->SetTitle("LP(#mu^{+})");
  newLPplus->GetXaxis()->SetTitle("LP(#mu)");
  newLPplus->GetXaxis()->SetRangeUser(-1.0,2.0);
  newLPplus->SetLineWidth(3);
  newLPplus->SetLineStyle(2);
  newLPplus->SetLineColor(kRed);
  newLPplus->DrawCopy();
  tplus->DrawClone("same");

//   c2->cd();
//   newLPminus->GetYaxis()->SetTitle("Events / 100 nb^{-1}");
//   newLPminus->GetXaxis()->SetTitle("LP(#mu^{-})");
  newLPminus->SetLineWidth(3);
  newLPminus->SetLineStyle(2);
  newLPminus->SetLineColor(kBlack);
  newLPminus->DrawCopy("same");
//   tplus->DrawClone("same");

  leg = new TLegend(0.762931,0.720339,0.8649425,0.8622881,NULL,"brNDC");
  leg->AddEntry(newLPplus,"#mu^{+}","l");
  leg->AddEntry(newLPminus,"#mu^{-}","l");

  leg->SetFillColor(kWhite);
  leg->SetBorderSize(0);
  leg->SetTextFont(62);
  leg->DrawClone();

  file0->Close();

  return;
}

void Drawbbc(){
	TFile *f = new TFile("merged.root");
	TH1D* hbbcs;
	TH1D* hbbcn;
	TCanvas *c1;
	TCanvas *c2;
	TLegend *leg = new TLegend(0.5,0.7,0.7,0.85);
	leg->SetBorderSize(0);
	leg->SetFillColor(0);
	leg->SetTextSize(0.048);
	TH2F* bbcsbbcn = (TH2F*)f->Get(Form("bbcsbbcn_0"));
	for(int xcent=1; xcent < 6; xcent++){
		TH1F* bbcsbbcndis_t = (TH1F*)f->Get(Form("bbcsbbcn_%d",xcent));
		bbcsbbcn->Add(bbcsbbcndis_t);
	}
	hbbcs = (TH1D*)bbcsbbcn->ProjectionX("bbcs",0,-1);
	hbbcn = (TH1D*)bbcsbbcn->ProjectionY("bbcn",0,-1);
	hbbcs->Rebin(4);
	hbbcn->Rebin(4);
	c1= new TCanvas();
	c1->cd();
	c1->SetLogy();
	SetTitle(*hbbcs,"bbc charge","# of events","");
	SetRange(*hbbcs,0,1e-1,100,hbbcn->GetMaximum()*5);
	SetStyle(*hbbcs,1.2,1,20,0,0);
	SetStyle(*hbbcn,1.2,2,24,0,0);
	hbbcs->Draw("P");
	hbbcn->Draw("Psame");
	leg->AddEntry(hbbcs,"bbc south");
	leg->AddEntry(hbbcn,"bbc north");
	leg->Draw("same");
	c1->Print(Form("bbcdis.png"));
	c2= new TCanvas();
	c2->cd();
	c2->SetLogy();
	hbbcs_norm = (TH1D*)hbbcs->Clone("hbbcs_norm");
	hbbcn_norm = (TH1D*)hbbcn->Clone("hbbcn_norm");
	hbbcs_norm->Scale(1./hbbcs_norm->Integral());
	hbbcn_norm->Scale(1./hbbcn_norm->Integral());
	SetRange(*hbbcs_norm,0,1e-9,100,1);
	hbbcs_norm->Draw("P");
	hbbcn_norm->Draw("Psame");
	leg->Draw("same");
	c2->Print(Form("bbcdis_norm.png"));	
	c3= new TCanvas();
	c3->cd();
	c3->SetLogz();
	SetTitle(*bbcsbbcn,"bbc charge south","bbc charge north","# of events");
	SetRange(*bbcsbbcn,0,0,100,100);
	bbcsbbcn->Draw("colz");
	c3->Print("bbcsbbcn.png");
}

Double_t MakeHistsForSlides(TString goodFile, TString badFile, TString rootFile)
{
  //Extracting the histogram
  TFile* f = TFile::Open(rootFile);
  TDirectoryFile* dirFile = (TDirectoryFile*)(f->Get("AliAnalysisTaskCalibEmcal"));
  TList* list = (TList*)(dirFile->Get("histogram"));
  TH2D* cellVAmplitude = (TH2D*)(list->FindObject("_histogram_cell_id_amplitude"));

    
  //Getting a Good Cell
  TH1D* goodC = new TH1D();
  ifstream myFile(goodFile);
  TString title = "";

  const int num = 500;

  if(myFile.is_open())
    {
      string line = "";
      std::getline(myFile,line);

      int index = 0;
      double chi2 = 0.0;
      
      while(myFile >> index >> chi2 && index < num)
	{
	  //cout << index << "\t" << chi2 << endl;
	  title = Form("Energy distribution of Cell %d",index);
	  TH1D* goodCell = cellVAmplitude->ProjectionY("goodCell", index+1,index+1);
	  goodC = (TH1D*)goodCell->Clone();
	}
    }

void NTupleSS() {
    double Lumi = 19700.;
    double signal_weight = 2.39163;//Lumi/(3263180./1177.3);
    double ttbar_weight  = Lumi / (4246440. / 23.64);
    double tautau_weight = Lumi / (3297032. / 1966.7);
    double tbarw_weight  = Lumi / (493460. / 11.1);
    double tw_weight     = Lumi / (497658. / 11.1);
    double ww_weight     = Lumi / (10000430. / 54.94);
    double wz_weight     = Lumi / (10000280. / 33.21);
    double zz_weight     = Lumi / (9799908. / 17.7);

    TH1D* Data = GetDataPhiStar();
    TH1D* h_ee    = GetBGPhiStar(File_Signal_reco,    signal_weight);
    TH1D* h_zz    = GetBGPhiStar(File_zz,    zz_weight);
    TH1D* h_tt    = GetBGPhiStar(File_tt,    ttbar_weight);
    TH1D* h_tautau = GetBGPhiStar(File_tautau, tautau_weight);
    TH1D* h_tbarw = GetBGPhiStar(File_tbarw, tbarw_weight);
    TH1D* h_tw    = GetBGPhiStar(File_tw,    tw_weight);
    TH1D* h_ww    = GetBGPhiStar(File_ww,    ww_weight);
    TH1D* h_wz    = GetBGPhiStar(File_wz,    wz_weight);
    //  TH1D* h_zz    =GetBGPhiStar(File_zz,    zz_weight);

    double data_sel = Data->GetSumOfWeights();
    double ee_sel = h_ee->GetSumOfWeights();
    double tt_sel = h_tt->GetSumOfWeights();
    double tautau_sel = h_tautau->GetSumOfWeights();
    double tbarw_sel = h_tbarw->GetSumOfWeights();
    double tw_sel = h_tw->GetSumOfWeights();
    double ww_sel = h_ww->GetSumOfWeights();
    double wz_sel = h_wz->GetSumOfWeights();
    double zz_sel = h_zz->GetSumOfWeights();

    // double t_bg=ee_sel+tt_sel+tautau_sel+tbarw_sel+tw_sel+ww_sel+wz_sel+zz_sel;
    double t_bg = ee_sel + tt_sel + tautau_sel + tbarw_sel + tw_sel + ww_sel +
                  wz_sel + zz_sel;
    cout << "total left:" << data_sel - t_bg;
    cout << "Data: " << data_sel << "  ee: " << ee_sel << "  tt: " << tt_sel <<
         "  tautau: " << tautau_sel << "  tbarw: " << tbarw_sel << "  tw: " << tw_sel <<
         " singletop: " << tbarw_sel + tw_sel << "  ww: " << ww_sel << "  wz: " << wz_sel
         << "  zz: " << zz_sel << endl;
    cout << "ratio:"          << "  ee: " << ee_sel / data_sel << "  tt: " << tt_sel
         / data_sel << "  tautau: " << tautau_sel / data_sel << "  tbarw: " << tbarw_sel
         / data_sel << "  tw: " << tw_sel / data_sel << " singletop: " <<
         (tbarw_sel + tw_sel) / data_sel << "  ww: " << ww_sel / data_sel << "  wz: " <<
         wz_sel / data_sel << "  zz: " << zz_sel / data_sel << "data: " <<
         (data_sel - t_bg) / data_sel << endl;

    TH1D* h_qcd  = (TH1D*)Data->Clone();
    h_qcd->Add(h_ee, -1.0);
    h_qcd->Add(h_tt, -1.0);
    h_qcd->Add(h_tautau, -1.0);
    h_qcd->Add(h_tbarw, -1.0);
    h_qcd->Add(h_tw, -1.0);
    h_qcd->Add(h_ww, -1.0);
    h_qcd->Add(h_wz, -1.0);
    h_qcd->Add(h_zz, -1.0);

    h_qcd->Draw();
    cout << "norm=" << h_qcd->GetSumOfWeights() << endl;;
}

void CheckTaggEffButVsH() {
	Char_t Name[256];
	Char_t *VarName;

	//Connect to all existing histograms
	RootFileBut = new TFile(Str_RootFilesBut);
	TH2D *hBut2D = (TH2D*)gROOT->FindObject("hTaggerTime");
	RootFileH = new TFile(Str_RootFilesH);
	TH2D *hH2D = (TH2D*)gROOT->FindObject("hTaggerTime");

	gROOT->cd();

	TH1D *hButBg = hBut2D->ProjectionY("TempBut1D", 800, 950);
	TH1D *hHBg = hH2D->ProjectionY("TempH1D", 800, 950);


	TCanvas *c1 = new TCanvas();
	gStyle->SetPalette(1);
	c1->Divide(2,2);
	c1->cd(1);
	hButBg->Draw("");
	c1->cd(2);
	hHBg ->Draw("");
	c1->cd(3);
	delete gROOT->FindObject("Diff");
	TH1D *hDiff = (TH1D*) hButBg->Clone("Diff");
	hDiff->Divide(hHBg);
	hDiff->Draw();
	

	RootFileBut->Close();
	RootFileH->Close();

	printf("Finished.\n");
}

void MCefficiency(TString inputmc, TString outputfile, TString tfend, TString selmcgen, TString cut, TString weight, Float_t centmin, Float_t centmax)
{
  Float_t hiBinMin,hiBinMax;
  hiBinMin = centmin*2;
  hiBinMax = centmax*2;

  selmcgen = Form("%s&&hiBin>=%f&&hiBin<%f",selmcgen.Data(),hiBinMin,hiBinMax);
  cut = Form("%s&&hiBin>=%f&&hiBin<%f",cut.Data(),hiBinMin,hiBinMax);

  TFile* infMC = new TFile(inputmc.Data());
  TTree* ntMC = (TTree*)infMC->Get("ntDkpi");
  ntMC->AddFriend("ntHi");
  ntMC->AddFriend("ntSkim");
  TTree* ntGen = (TTree*)infMC->Get("ntGen");
  ntGen->AddFriend("ntHi");
  ntGen->AddFriend("ntSkim");

  TH1D* hPtMC = new TH1D("hPtMC","",nPtBins,ptBins);
  TH1D* hPtGen = new TH1D("hPtGen","",nPtBins,ptBins);

  ntMC->Project("hPtMC","Dpt",TCut(weight)*(TCut(cut.Data())&&"(Dgen==23333)"));
  divideBinWidth(hPtMC);
  ntGen->Project("hPtGen","Gpt",TCut(weight)*(TCut(selmcgen.Data())));
  divideBinWidth(hPtGen);
  TH1D* hEff = (TH1D*)hPtMC->Clone("hEff");
  hEff->Divide(hPtGen);

  TFile* fout=new TFile(Form("%s_cent_%.0f_%.0f_%s.root",outputfile.Data(),centmin,centmax,tfend.Data()),"recreate");
  fout->cd();
  hEff->Write();
  fout->Close();
}

TH1D * scaleHistsFromPtHatRootFilesReturnHist( string histname, vector <TFile* > &rootFiles )
{
    const double lumi = 40.03 ; // pb-1
    vector <double> mcscales;
    // mcscales.push_back(  117276.          / 3468514. * lumi);  // xs / nevents * lumi
    // mcscales.push_back(    7823.          / 2936644. * lumi);
    mcscales.push_back(     648.2         / 1971800. * lumi);
    mcscales.push_back(     186.9         / 1981608. * lumi);
    mcscales.push_back(      32.293       / 1990208. * lumi);
    mcscales.push_back(       9.4183      / 1487712. * lumi);
    mcscales.push_back(       0.84265     /  197959. * lumi);
    mcscales.push_back(       0.114943    /  194924. * lumi);
    mcscales.push_back(       0.00682981  /  198383. * lumi);
    mcscales.push_back(       0.000165445 /  194528. * lumi);

    TH1D * htemp = (TH1D*)  rootFiles[0] -> Get( histname.c_str() );
    TH1D * h_sum = (TH1D*)htemp->Clone();
    h_sum->Reset();
    h_sum->Sumw2();

    for (unsigned int i=0; i<rootFiles.size(); i++) {
        TH1D * h1 = (TH1D*)  rootFiles[i] -> Get( histname.c_str() );
        // h1->Sumw2();
        h1->Scale(mcscales[i]);
        h_sum->Add(h1);
        // cout<<mcscales[i]<<"  "<<h1->Integral()<<endl;
    }
    return h_sum;
}

// Convert plots to paper format - no title, bigger fonts etc
void combineHists( TFile* fSig, TFile* fBg, TFile* fBg2, std::string histName, std::string plotOpt, std::string outputName, std::vector<double> scalingFactors, std::string label, std::string yTitle="DEFAULT"){
    TH1::SetDefaultSumw2();

    // for 3 hists - sig and 2 bg
    TCanvas c1;
    TH1D* hSig = fSig->Get(histName.c_str());
    hSig->SetLineColor(kRed);
    hSig->SetMarkerSize(0);
    doSignalHist(hSig);
    
    // Make combined BG hist
    // Need to rescale carefully
    std::vector<TFile*> files;
    files.push_back(fBg);
    files.push_back(fBg2);
    // TH1D* hBg = combine(files, histName, scalingFactors);
    TH1D* hBgA = fBg->Get(histName.c_str());
    TH1D* hBgB =  fBg2->Get(histName.c_str());
    TH1D* hBg = (TH1D*) hBgA->Clone();
    double total = scalingFactors[0]+scalingFactors[1];
    hBg->Scale(scalingFactors[0]/total);
    hBg->Add(hBgB, scalingFactors[1]/total);
    hBg->SetMarkerSize(0);
    doAltBGHist(hBg);
    
    THStack st("h","");
    st.Add(hSig);
    st.Add(hBg);
    st.Draw((plotOpt+"NOSTACK").c_str());
    st.GetXaxis()->SetTitle(hSig->GetXaxis()->GetTitle());
    if (yTitle == "DEFAULT") {
        st.GetYaxis()->SetTitle(hSig->GetYaxis()->GetTitle());
    } else {
        st.GetYaxis()->SetTitle(yTitle.c_str());
    }
    setAltTitleLabelSizes(&st.GetHistogram());
    st.SetTitle("");
    st.Draw((plotOpt+"NOSTACK").c_str());

    TLegend* l_all = new TLegend(0.65,0.6,0.89,0.89);
    l_all->AddEntry(hBg,"Gen. level QCD MC","lp");
    l_all->AddEntry((TObject*)0,"(b#bar{b} + q-g scatter,",""); //null pointers for blank entries
    l_all->AddEntry((TObject*)0,"q = b, #bar{b}, c, #bar{c})","");
    l_all->AddEntry(hSig, "Signal MC", "lp");
    l_all->AddEntry((TObject*)0,"m_{#phi} = 8 GeV", "");
    doStandardLegend(l_all);
    l_all->Draw();

    TPaveText t(0.15, 0.75, 0.5, 0.85, "NDC");
    t.AddText(label.c_str());
    doStandardText(&t);
    if (label != "") {
        t.Draw();
    }
    c1.SaveAs(outputName.c_str());

    if (!hSig) delete hSig;
    if (!hBg) delete hBg;
}

void  makeGaussianSignals(SigData_t& m_sigdata)
{
  //std::vector<TH1D *> vgsh(NUMCHAN);
  std::vector<TH1D *> vcdfh(NUMCHAN);

  if( m_sigdata.find("gs") == m_sigdata.end() ) {
    cerr << "Gaussian signal data not found, not making CDF signal!" << endl;
    return;
  }

  for (int ichan=0; ichan<NUMCHAN; ichan++) {
    TH1D *cdfh;

    TString channame(channames[ichan]);
    TString name;

    TH1D * gsh = m_sigdata["gs"].at(ichan);

    assert(gsh) ;

#if 0
    name = "Signalgs_"+channame;
    gsh = (TH1D *)tch->Clone(name.Data());

    assert(gsh);

    gsh->SetTitle("Gaussian signal");
    
    gsh->Reset();

    TF1 *g = (TF1 *)gROOT->GetFunction("gaus");
    g->SetParameters(1,gaussian_mean_mass_gev,gaussian_mass_sigma_gev);
    gsh->FillRandom("gaus",100000);

    // norm to 1pb signal with 1/fb integrated luminosity
    double norm = 1000 * gseffxacc[ichan]/gsh->Integral(0,gsh->GetNbinsX()+1,"width");

    //gsh->Scale(norm/eff_fudgefactor); // kludge: pre-undo the fudge in the next module
    gsh->Scale(norm);

    vgsh[ichan] = gsh;
#endif

    // New CDF bump, same as Gauss but set to CDF (obs/theor)*(LHC theor) = 3.43pb
    cdfh = (TH1D *)gsh->Clone("CDFbump");

    cdfh->Scale(3.43);

    vcdfh[ichan] = cdfh;

    cdfh->Draw();

    gsh->Draw("same");

  } // channel loop

  //m_sigdata["gs"]  = vgsh;
  m_sigdata["cdf"] = vcdfh;
}                                                           // makeGaussianSignals

void loglikdistrib(Int_t ntrials = 10000, Bool_t print = kFALSE)
{
  // compute distribution of log likelihood value
  TH1D * hmc   = gStack[gPadNr][gOrder[gPadNr][0]];
  TH1D * hdata = gStack[gPadNr][gMaxProcess-1];
  Int_t nbins = hmc->GetNbinsX();
  Double_t loglik = loglikelihood(hmc, hdata, 1, nbins);
  TH1D * htest = new TH1D(*hdata);
  TH1D * lldistrib = new TH1D("lldistrib", "log(Likelihood) distribution", 
			      1000, loglik-200, loglik+200);
  setopt(lldistrib);
  for (Int_t n = 0; n < ntrials; n++) {
    // generate poisson around theorie
    for (Int_t i = 1; i <= nbins; i++) {
      htest->SetBinContent(i, gRandom->Poisson(hmc->GetBinContent(i)));
    }
    lldistrib->Fill(loglikelihood(hmc, htest, 1, nbins));
  }
  TCanvas * llcanvas = new TCanvas("llcanvas", "Log(Likelihood) distribution", 
				   40, 40, 800, 600);
  setopt(llcanvas);
  lldistrib->SetFillColor(kYellow);
  lldistrib->Draw();
  lldistrib->GetYaxis()->SetTitle("Anzahl Ereignisse");
  lldistrib->GetXaxis()->SetTitle("-ln L");
  // autozoom
  Int_t lowbin = 1;
  while (lldistrib->GetBinContent(lowbin) == 0)
    lowbin++;
  Int_t highbin = lldistrib->GetNbinsX();
  while (lldistrib->GetBinContent(highbin) == 0)
    highbin--;
  lldistrib->SetAxisRange(lldistrib->GetBinLowEdge(lowbin), 
			  lldistrib->GetBinLowEdge(highbin));
  TH1D * hworse = (TH1D *) lldistrib->Clone();
  for (Int_t nbin = 1; nbin < 501; nbin++) {
    hworse->SetBinContent(nbin, 0);
  }
  hworse->SetFillColor(95);
  hworse->Draw("same");
  Double_t pvalue = lldistrib->Integral(501,1000) / lldistrib->Integral();
  TLatex * tex = new TLatex(0.18, 0.96, Form("-ln L_{obs} = %5.2f", loglik));
  tex->SetNDC();
  tex->SetTextAlign(13);
  tex->Draw();
  tex = new TLatex(0.18, 0.86, Form("CL_{obs} = %.3f", pvalue));
  tex->SetNDC();
  tex->SetTextAlign(13);
  tex->Draw();
  TLine * l = new TLine(loglik, 0, loglik, lldistrib->GetMaximum());
  l->SetLineWidth(3);
  l->SetLineColor(kBlue);
  l->Draw();
  llcanvas->Modified();
  llcanvas->Update();
  if (print)
    llcanvas->Print("lldistrib.pdf");
  cd(gPadNr+1);
}

void getComparisons_ABCD(TString INPUTDIR_PREFIX, TString varType, TString whichHist, TString SCEN_TRIG, TString shapeComb, TString RR) {

   //Add categories according to the factorization ansats. Fixed
   const int nconfs = 3;
   TString overallConfs[nconfs] = {"Glb_Tight2012_IsolPFRelCombNoEGammaR03_"+SCEN_TRIG,"Glb_Tight2012_IsolPFRelCombNoEGammaR03","Glb_Tight2012"};

   //get TP efficiencies for above categories
   TH1D* htrigAB;
   TH1D* hisoAB;
   TH1D* hrecoidAB;
   TH1D* htrig_mc;
   TH1D* hiso_mc;
   TH1D* hrecoid_mc; 

   if (whichHist == "trig")     htrigAB   = getSimplePlot(INPUTDIR_PREFIX,overallConfs[0],"data_all",varType,SCEN_TRIG,"RunAB");
   if (whichHist == "iso")      hisoAB    = getSimplePlot(INPUTDIR_PREFIX,overallConfs[1],"data_all",varType,SCEN_TRIG,"RunAB");
   if (whichHist == "recoid")   hrecoidAB = getSimplePlot(INPUTDIR_PREFIX,overallConfs[2],"data_all",varType,SCEN_TRIG,"RunAB");

   if (whichHist == "trig")     htrigC   = getSimplePlot(INPUTDIR_PREFIX,overallConfs[0],"data_all",varType,SCEN_TRIG,"RunC");
   if (whichHist == "iso")      hisoC    = getSimplePlot(INPUTDIR_PREFIX,overallConfs[1],"data_all",varType,SCEN_TRIG,"RunC");
   if (whichHist == "recoid")   hrecoidC = getSimplePlot(INPUTDIR_PREFIX,overallConfs[2],"data_all",varType,SCEN_TRIG,"RunC");

   if (whichHist == "trig")     htrigD   = getSimplePlot(INPUTDIR_PREFIX,overallConfs[0],"data_all",varType,SCEN_TRIG,"RunD");
   if (whichHist == "iso")      hisoD    = getSimplePlot(INPUTDIR_PREFIX,overallConfs[1],"data_all",varType,SCEN_TRIG,"RunD");
   if (whichHist == "recoid")   hrecoidD = getSimplePlot(INPUTDIR_PREFIX,overallConfs[2],"data_all",varType,SCEN_TRIG,"RunD");

   if (whichHist == "trig")     htrig_mc   = (TH1D*)htrigAB->Clone();
   if (whichHist == "iso")      hiso_mc    = (TH1D*)hisoAB->Clone(); 
   if (whichHist == "recoid")   hrecoid_mc = (TH1D*)hrecoidAB->Clone();

   if (varType != "rrr" && varType != "rrr2" && varType != "rrr3") {    
     if (whichHist == "trig")   { 
        htrig_mc   = getSimplePlot(INPUTDIR_PREFIX,overallConfs[0],"datalike_mc",varType,SCEN_TRIG,RR);
        getComparisonABCDratio(htrig_mc, htrigAB, htrigC, htrigD, varType, true, INPUTDIR_PREFIX,RR);
     } 
     if (whichHist == "iso")    {
        hiso_mc    = getSimplePlot(INPUTDIR_PREFIX,overallConfs[1],"datalike_mc",varType,SCEN_TRIG,RR);
        getComparisonABCD(hiso_mc, hisoAB, hisoC, hisoD, varType, true, INPUTDIR_PREFIX,RR);
     }
     if (whichHist == "recoid") { 
        hrecoid_mc = getSimplePlot(INPUTDIR_PREFIX,overallConfs[2],"datalike_mc",varType,SCEN_TRIG,RR);
        getComparisonABCD(hrecoid_mc, hrecoidAB, hrecoidC, hrecoidD, varType,true, INPUTDIR_PREFIX,RR);
     }
   }
}

	// 1 file argument -> add e+m from same files e.g. mc
	// this is considered to be MC -> added to sum histo
	TH1D* get(const TString& File){ //
		TFile file(path+File);
		TH1D* mH = (TH1D*)file.Get(mName);
		TH1D* eH = (TH1D*)file.Get(eName);
		gROOT->cd();
		TH1D*  H = mH->Clone();// Sumw2() is copied as well
		H->Add(eH);
		return H;
	}

// -----------------------------------------------------------------------------
//
TH1* getHisto( std::string nameFile,
	       std::string nameHist,
	       std::string Dirname, 
	       int rebin ) {
  std::string name = nameFile;
  TFile* file =  new TFile(name.c_str());
  if (file) { std::cout << "Opened file: " << file->GetName() << std::endl; }
  else { 
    std::cout << "Could not find file: " << name << std::endl; 
    return 0; 
  }
  
  TDirectory* dir = (TDirectory*)file->Get(Dirname.c_str());
  if (dir) { std::cout << "Opened dir: " << dir->GetName() << std::endl; }
  else { 
    std::cout << "Could not find dir: " << Dirname << std::endl; 
    return 0; 
  }
  
  int low = 375;
  TH1* hist = 0;
  if ( false || nameHist.find("HtMultiplicity_HT375") == std::string::npos ) { 
    hist = (TH1*)dir->Get(nameHist.c_str());
  } else {
    
    for ( uint ii = low; ii <= 975; ii+=100 ) {
      std::stringstream tmp; tmp << "HtMultiplicity_HT" << ii << nameHist.substr(20);
      if ( !hist ) { 
	dir->cd();
	TH1D* temp = (TH1D*)dir->Get( "HtMultiplicity_HT375_aT0" );
	//TH1D* temp = (TH1D*)file->Get( tmp.str().c_str() );
	if (temp) { hist = (TH1D*)temp->Clone(); } 
	else { std::cout << "1 Unable to retrieve histo with name " << tmp.str() << std::endl; }
      } else { 
	dir->cd();
	TH1D* temp = (TH1D*)dir->Get( tmp.str().c_str() );
	if (temp) { hist->Add( (TH1D*)temp ); } 
	else { std::cout << "2 Unable to retrieve histo with name " << tmp.str() << std::endl; }
      }
    }

  }
  if (hist) { std::cout << "Opened histo: " << hist->GetName() << std::endl; }
  else { 
    std::cout << "Could not find histo: " << nameHist << std::endl; 
    return 0; 
  }

  hist->SetLineWidth(3);
  if ( rebin > 0 ) { hist->Rebin(rebin); }
  hist->GetXaxis()->SetTitleSize(0.055);
  hist->GetYaxis()->SetTitleSize(0.055);
  hist->GetXaxis()->SetLabelSize(0.05);
  hist->GetYaxis()->SetLabelSize(0.05);
  hist->SetStats(kFALSE);
  return hist;
}

int main(){
	TFile * data = TFile::Open("TreesMu_Data_plots.root");
	TH1D * dataMtop = (TH1D*)data->Get("antiEtaFwDTrue_allW/antiEtaFwDTrue_allWcosTheta");
	
	TFile * tt = TFile::Open("TreesMu_TTBar_RW.root");
	
	TH1D * ttMtop = (TH1D*)tt->Get("MtopOutWindowTrue_allW/MtopOutWindowTrue_allWcosTheta");
	
	TH2D * ttMtop2D = (TH2D*)tt->Get("MtopOutWindowTrue_allW/MtopOutWindowTrue_allWcosTheta2D");
	
  std::pair<TF1, WeightFunctionCreator*> WeightFuncUD(WeightFunctionCreator::getWeightFunction("WeightFuncUDF", F0, FL));
  std::pair<TF1, WeightFunctionCreator*> WeightFuncDU(WeightFunctionCreator::getWeightFunction("WeightFuncDUF", F0, FL));
  cout<<F0 + F0Sys<<"\t"<<FL - FLSys<<endl;
	WeightFuncUD.first.SetParameters(F0 + F0Sys, FL - FLSys);
	WeightFuncDU.first.SetParameters(F0 - F0Sys, FL + FLSys);
	
	TH1D * ttUD = myReweightor(ttMtop2D,WeightFuncUD,"WeightFuncUD_");
	cout<<"ratio: "<<ttMtop2D->Integral()/ttUD->Integral()<<endl;
	ttUD->Scale(ttMtop2D->Integral()/ttUD->Integral());
	ttUD->Add(ttMtop);
	TH1D * wUD = (TH1D*)ttUD->Clone("wUD");
	wUD->Scale(-1.);
	wUD->Add(dataMtop);
	
	TH1D * ttDU = myReweightor(ttMtop2D,WeightFuncDU,"WeightFuncDU_");
	cout<<"ratio: "<<ttMtop2D->Integral()/ttDU->Integral()<<endl;
	ttDU->Scale(ttMtop2D->Integral()/ttDU->Integral());
	ttDU->Add(ttMtop);
	TH1D * wDU = (TH1D*)ttDU->Clone("wDU");
	wDU->Scale(-1.);
	wDU->Add(dataMtop);
	
	TFile * out = new TFile("file.root","recreate");
	out->cd();
	WeightFuncUD.first.Write();
	WeightFuncDU.first.Write();
	ttUD->Write();
	ttDU->Write();	
	wUD->Write();
	wDU->Write();	
	out->Close();
	return 1;
}