C++ TFile::WriteTObject方法代码示例

bool TrGainDB::Save(const char* filename) {
  Init();
  TFile* rootfile = TFile::Open(filename,"recreate");
  if (rootfile==0) return false;
  rootfile->WriteTObject(GetHead());
  return true;
}

void MakeHQTHiggsPtDistribution( string inputFilename ,  const string Label) {

  //--------------------------------------------------------------------------------------------------------------
  // Settings 
  //==============================================================================================================
  
  Double_t lumi;              // luminosity (pb^-1)
  string label = "";
  if (Label != "") label = "_" + Label;
 


  //--------------------------------------------------------------------------------------------------------------
  // Main analysis code 
  //==============================================================================================================  
  
  //
  // Access samples and fill histograms
  //  
  ifstream infile(inputFilename.c_str());
  vector<double> lowEdges;
  vector<double> binContent;
  double x_min;
  double x_max;
  double y;
  Int_t binNumber = 1;
  while(infile >> x_min) {
    infile >> x_max >> y;
    lowEdges.push_back(x_min);
    binContent.push_back(y);
  }

  Double_t *xLowEdges = new Double_t[lowEdges.size()];
  for (UInt_t i=0; i < lowEdges.size();i++) {
    xLowEdges[i] = lowEdges[i];
  }


  //--------------------------------------------------------------------------------------------------------------
  // Histograms
  //==============================================================================================================  
  TH1D *bosonSystemPt = new TH1D((string("HiggsPt")+ label).c_str(), ";  Higgs p_{T} [GeV/c]; Number of Events", lowEdges.size()-1, xLowEdges);

  for (UInt_t i=1;i<lowEdges.size();i++) {
    bosonSystemPt->SetBinContent(i,binContent[i]);
  }

  TFile *file = new TFile("HwwAcceptanceSystematics.root", "UPDATE");
  file->WriteTObject(bosonSystemPt, bosonSystemPt->GetName(), "WriteDelete");
  file->Close();

}

//#include <typeinfo.h>
void addFlatNuisances(std::string fi){
  gSystem->Load("libHiggsAnalysisCombinedLimit.so");
  TFile *fin = TFile::Open(fi.c_str());
  RooWorkspace *wspace = (RooWorkspace*)fin->Get("w_hmumu");

  wspace->Print("");

  RooStats::ModelConfig *mc = (RooStats::ModelConfig*)wspace->genobj("ModelConfig");
  RooArgSet *nuis = (RooArgSet*) mc->GetNuisanceParameters();
  std::cout << "Before...." << std::endl;
  nuis->Print();
  
  RooRealVar *mgg = (RooRealVar*)wspace->var("mmm");
  // Get all of the "flat" nuisances to be added to the nusiances:
  RooArgSet pdfs = (RooArgSet) wspace->allVars();
  RooAbsReal *pdf;
  TIterator *it_pdf = pdfs.createIterator();
  

  while ( (pdf=(RooAbsReal*)it_pdf->Next()) ){
  	  if (!(std::string(pdf->GetName()).find("zmod") != std::string::npos )) {
  	   if (!(std::string(pdf->GetName()).find("__norm") != std::string::npos )) {
	   	continue;
	   }
	  }
	  pdf->Print();
	  RooArgSet* pdfpars = (RooArgSet*)pdf->getParameters(RooArgSet(*mgg));
	  pdfpars->Print();

	  std::string newname_pdf = (std::string("unconst_")+std::string(pdf->GetName()));
	  wspace->import(*pdf,RooFit::RenameVariable(pdf->GetName(),newname_pdf.c_str()));
	  pdf->SetName(newname_pdf.c_str());
	  nuis->add(*pdf);
  }
 
  wspace->var("MH")->setVal(125.0);
  std::cout << "After..." << std::endl;
  nuis->Print();
  mc->SetNuisanceParameters(*nuis);
  //RooWorkspace *wspace_new = wspace->Clone();
  //mc->SetWorkspace(*wspace_new);
  //wspace_new->import(*mc,true);

  TFile *finew = new TFile((std::string(fin->GetName())+std::string("_unconst.root")).c_str(),"RECREATE");
  //wspace_new->SetName("w");
  finew->WriteTObject(wspace);
  finew->Close();
}

void PileUpMaker() {
  // Access data file
  TFile *file = TFile::Open("ZmumuGammaNtuple_Full2012_MuCorr.root");
  TTree *Datatree = (TTree*)file->Get("ZmumuGammaEvent");

  // Set up histogram for data pile up
  TH1F *pileupraw = new TH1F("pileup","",40,0,80);
 
  // Fill pile up distribution
  Int_t nentries = Datatree->GetEntriesFast();
  Int_t nbytes = 0; 
  for (Int_t i=0; i<nentries;i++) {
    nbytes += Datatree->GetEvent(i);
    UInt_t NVertices;
    Datatree->SetBranchAddress("NVertices",&NVertices);
    pileupraw->Fill(float(NVertices)/0.7);
  }

  // Normalize the distribution
  TH1F* PU = (TH1F*)pileupraw->Clone("pileup");
  PU->Scale(1./PU->Integral());

  // Plot both the normalized and unnormalized pile up
  TCanvas *cv = 0;
  cv= new TCanvas("cv","cv",800, 600);
  pileupraw->Draw("");
  cv->SaveAs("DataPU.pdf");

  cv= new TCanvas("cv","cv",800,600);
  PU->Draw("");
  cv->SaveAs("DataPU_normalized.pdf");

  // Save the normalized distribution to a root file for weighting MC
  TFile *file = TFile::Open("2012_PhosphorData_PileUp.root","UPDATE");
  file->cd();
  file->WriteTObject(PU,PU->GetName(),"WriteDelete");
  file->Close();
  delete file;
}

//--------------------------------------------------------------------------------------------------
// Main of the macro
//--------------------------------------------------------------------------------------------------
void computeFakeRates(char *FakeRateFilename ) {

  string fakeTypeString = "";
  //Samples
  vector<string> madgraphDatasetFiles;
  vector<string> madgraphDatasetNames;
  vector<string> pythiaDatasetFiles;
  vector<string> pythiaDatasetNames;

  //*****************************************************************************************
  //Do FakeRate Computation 
  //*****************************************************************************************

  //Create root file to store fake rates
  TFile *file = new TFile(FakeRateFilename, "RECREATE");
  file->cd();

  //do both muon and electrons
  Double_t maxY = 1.0;
  for (int faketype=11 ; faketype <= 13; faketype+=2) {
    if (faketype == 11) {
      fakeTypeString = "Electron";
      maxY = 0.5;
    } else if (faketype == 13) {
      fakeTypeString = "Muon";
      maxY = 0.01;
    } else {
      cerr << "Error: faketype = " << faketype << " not recognized.\n";
      return;
    }
    
    //*****************************************************************************************
    //Load input histogram files
    //*****************************************************************************************
    madgraphDatasetFiles.clear(); 
    madgraphDatasetNames.clear();
    pythiaDatasetFiles.clear();
    pythiaDatasetNames.clear();
    
    //Madgraph WJets
    madgraphDatasetFiles.push_back("/home/sixie/hist/Compute" + fakeTypeString + "FakeRate/filler/006/Compute" + fakeTypeString + "FakeRate_f8-wjets-mg-id9.root");
    madgraphDatasetNames.push_back("f8-wjets-mg-id9");
    //Pythia WJets
    pythiaDatasetFiles.push_back("/home/sixie/hist/Compute" + fakeTypeString + "FakeRate/filler/006/Compute" + fakeTypeString + "FakeRate_s8-we-id9.root");
    pythiaDatasetNames.push_back("s8-we-id9");
    pythiaDatasetFiles.push_back("/home/sixie/hist/Compute" + fakeTypeString + "FakeRate/filler/006/Compute" + fakeTypeString + "FakeRate_s8-wm-id9.root");
    pythiaDatasetNames.push_back("s8-wm-id9");
    pythiaDatasetFiles.push_back("/home/sixie/hist/Compute" + fakeTypeString + "FakeRate/filler/006/Compute" + fakeTypeString + "FakeRate_s8-wtau-id9.root");
    pythiaDatasetNames.push_back("s8-wtau-id9");

    //*****************************************************************************************
    //Construct the fake rate Histgrams. Fit with function.
    //*****************************************************************************************
    vector<TH1F*> fakeRateHistograms; fakeRateHistograms.clear();
    vector<TF1*> fakeRateFits; fakeRateFits.clear();
    vector<string> fakeRateLegendNames; fakeRateLegendNames.clear();

    //pythia version
    TH1F *FakeRatePt_pythia = createFakeRatePtHist(pythiaDatasetFiles, pythiaDatasetNames, faketype,
                                                   (fakeTypeString + "FakeRate_Pt_pythia").c_str());
    FakeRatePt_pythia->GetYaxis()->SetTitle((fakeTypeString + " Fake Rate").c_str());    
    TF1 *FakeRateFitFunction_pythia = fitFakeRatePtHist(FakeRatePt_pythia);
    FakeRateFitFunction_pythia->SetName((fakeTypeString + "FakeRateFunction_Pt_pythia").c_str());
    
    fakeRateHistograms.push_back(FakeRatePt_pythia);
    fakeRateFits.push_back(FakeRateFitFunction_pythia);
    fakeRateLegendNames.push_back("Pythia W+Jets");

    //madgraph version
    TH1F *FakeRatePt_madgraph = createFakeRatePtHist(
      madgraphDatasetFiles, madgraphDatasetNames,faketype,
      (fakeTypeString + "FakeRate_Pt_madgraph").c_str());
    FakeRatePt_madgraph->GetYaxis()->SetTitle((fakeTypeString + " Fake Rate").c_str());    
    TF1 *FakeRateFitFunction_madgraph = fitFakeRatePtHist(FakeRatePt_madgraph);
    FakeRateFitFunction_madgraph->SetName(
      (fakeTypeString + "FakeRateFunction_Pt_madgraph").c_str());

    fakeRateHistograms.push_back(FakeRatePt_madgraph);
    fakeRateFits.push_back(FakeRateFitFunction_madgraph);
    fakeRateLegendNames.push_back("Madgraph W+Jets");

    cout << "pythia at x = 30 : " << FakeRateFitFunction_pythia->Eval(30.0) << endl;
    cout << "pythia at x = 30 : " << FakeRateFitFunction_madgraph->Eval(30.0) << endl;

    //*****************************************************************************************
    //Draw the plots
    //*****************************************************************************************
    drawFakeRatePlots(fakeRateHistograms,fakeRateFits,fakeRateLegendNames,
                      (fakeTypeString + "FakeRateEt").c_str(),maxY);

    //*****************************************************************************************
    //save them to the output file
    //*****************************************************************************************
    for (UInt_t i=0; i<fakeRateHistograms.size();i++) {
      fakeRateHistograms[i]->SetDirectory(file);
      file->WriteTObject(fakeRateFits[i]);
    }
//.........这里部分代码省略.........

//.........这里部分代码省略.........
                    }
                    if (OSSFweight[1]>0){
                      totalOS *= OSSFweight[1];
                      totalOSerror += pow(OSSFweighterror[1]/OSSFweight[1], 2);
                    }
                    totalOSerror = sqrt(totalOSerror) * totalOS;
                    if (isocut==2){
                      sumOS = OSSFweight[0] + OSSFweight[1];
                      sumOSerror = sqrt(pow(OSSFweighterror[0], 2) + pow(OSSFweighterror[1], 2));
                    }
                    if (iCR==2){
                      totalOS*=MC_weight;
                      sumOS*=MC_weight;
                      totalOSerror*=MC_weight;
                      sumOSerror*=MC_weight;
                    }
                    wgt = totalOS;
                    wgterror = totalOSerror;
                  }
                  else{
                    if (iCR==2){
                      totalOS*=MC_weight;
                      sumOS*=MC_weight;
                    }
                  }

                  hmZZTxyOS->Fill(strdraw_mZZ, strdraw_Txy, totalOS);
                  hmZZmZ1mZ2OS->Fill(strdraw_mZZ, strdraw_mZ1, strdraw_mZ2, totalOS);

                  if (iwgt>0){
/*
                    int binx = hmZZTxyOS->GetXaxis()->FindBin(strdraw_mZZ);
                    int biny = hmZZTxyOS->GetYaxis()->FindBin(strdraw_Txy);
                    hmZZTxyOS->SetBinError(binx, biny, sqrt(pow(hmZZTxyOS->GetBinError(binx, biny), 2)+pow(totalOSerror, 2)));

                    binx = hmZZmZ1mZ2OS->GetXaxis()->FindBin(strdraw_mZZ);
                    biny = hmZZmZ1mZ2OS->GetYaxis()->FindBin(strdraw_mZ1);
                    int binz = hmZZmZ1mZ2OS->GetZaxis()->FindBin(strdraw_mZ2);
                    hmZZmZ1mZ2OS->SetBinError(binx, biny, binz, sqrt(pow(hmZZmZ1mZ2OS->GetBinError(binx, biny, binz), 2)+pow(totalOSerror, 2)));
*/
                    if (isocut==2){
                      hmZZTxyOS_extra->Fill(strdraw_mZZ, strdraw_Txy, sumOS);
                      hmZZmZ1mZ2OS_extra->Fill(strdraw_mZZ, strdraw_mZ1, strdraw_mZ2, sumOS);
/*
                      int binx = hmZZTxyOS_extra->GetXaxis()->FindBin(strdraw_mZZ);
                      int biny = hmZZTxyOS_extra->GetYaxis()->FindBin(strdraw_Txy);
                      hmZZTxyOS_extra->SetBinError(binx, biny, sqrt(pow(hmZZTxyOS_extra->GetBinError(binx, biny), 2)+pow(sumOSerror, 2)));

                      binx = hmZZmZ1mZ2OS_extra->GetXaxis()->FindBin(strdraw_mZZ);
                      biny = hmZZmZ1mZ2OS_extra->GetYaxis()->FindBin(strdraw_mZ1);
                      int binz = hmZZmZ1mZ2OS_extra->GetZaxis()->FindBin(strdraw_mZ2);
                      hmZZmZ1mZ2OS_extra->SetBinError(binx, biny, binz, sqrt(pow(hmZZmZ1mZ2OS_extra->GetBinError(binx, biny, binz), 2)+pow(sumOSerror, 2)));
*/
                    }
                  }

                }

                hsip->Fill(strdraw, wgt);
                if (iwgt>0){
                  int binx = hsip->GetXaxis()->FindBin(strdraw);
                  hsip->SetBinError(binx, sqrt(pow(hsip->GetBinError(binx), 2)+pow(wgterror, 2)));
                }

              }

              if (iCR!=1 && isocut>1 && catZ2<3){
                frecord->WriteTObject(hmZZTxyOS);
                delete hmZZTxyOS;
                frecord->WriteTObject(hmZZmZ1mZ2OS);
                delete hmZZmZ1mZ2OS;
                if (isocut==2 && iwgt>0){
                  frecord->WriteTObject(hmZZTxyOS_extra);
                  delete hmZZTxyOS_extra;
                  frecord->WriteTObject(hmZZmZ1mZ2OS_extra);
                  delete hmZZmZ1mZ2OS_extra;
                }
              }

              double totalerror = 0;
              double totalintegral = hsip->IntegralAndError(1,hsip->GetNbinsX(),totalerror);

              cout << '\t' << totalintegral << '\t' << totalerror;

              delete hsip;
            }
            cout << endl;
          }
        }
        cout << endl;
      }
      cout << endl;
    }
    cout << endl;
  }
  delete tc;

  for (int ww=0; ww<2; ww++) fweight[ww]->Close();
  frecord->Close();
}

//.........这里部分代码省略.........
                               RooFit::Save());

      cout << "Fitted parameters\n";
      cout << mean->getVal() << endl;
      cout << sigma->getVal() << endl;
     
      if (i >= 1 && ( j >= 1 && j <= NEtaBins )) {
        GaussParamArray_Muons_mean->SetBinContent(i,j,mean->getVal());
        GaussParamArray_Muons_sigma->SetBinContent(i,j,sigma->getVal());
      } else {
        GaussParamArray_Muons_mean->SetBinContent(i,j,0);
        GaussParamArray_Muons_sigma->SetBinContent(i,j,0);        
      }




      //Save Workspace
      RooWorkspace *w = new RooWorkspace(Form("LeptonPtResolutionModel_Muons_PtBin%d_EtaBin%d",i,j),Form("LeptonPtResolutionModel_Muons_PtBin%d_EtaBin%d",i,j));
      w->import(*model);
      w->import(*data);
      //w->Print();

      //Make Plot
      RooPlot *frame = leptonPtRes.frame(RooFit::Bins(100));
      data->plotOn(frame,RooFit::MarkerStyle(kFullCircle),RooFit::MarkerSize(0.8),RooFit::DrawOption("ZP"));    
      model->plotOn(frame);

      TCanvas *cv = new TCanvas("cv","cv",800,600);
      
      frame->Draw();
      cv->SaveAs(Form("LeptonPtResolution_Muons%s_PtBin%d_EtaBin%d.gif",label.c_str(),i,j)); 

      fileOutput->WriteTObject(model, Form("LeptonPtResolutionModel_Muons_PtBin%d_EtaBin%d",i,j), "WriteDelete");
      fileOutput->WriteTObject(cv, Form("LeptonPtResolutionFit_Muons_PtBin%d_EtaBin%d",i,j), "WriteDelete");
      fileOutput->WriteTObject(w, w->GetName(), "WriteDelete");

      fileOutput->WriteTObject(GaussParamArray_Muons_mean, "GaussParamArray_Muons_mean", "WriteDelete");
      fileOutput->WriteTObject(GaussParamArray_Muons_sigma, "GaussParamArray_Muons_sigma", "WriteDelete");


    }
  }

  //********************************************************
  // Produce output lookup table
  //******************************************************** 
  ofstream outf_e("FakeMuonResponseMap.h");

  outf_e << "UInt_t FindFakeMuonResponseBin( double value, double bins[], UInt_t nbins) {" << endl;
  outf_e << "  UInt_t nbinboundaries = nbins+1;" << endl;
  outf_e << "  UInt_t bin = 0;" << endl;
  outf_e << "  for (uint i=0; i < nbinboundaries; ++i) {" << endl;
  outf_e << "    if (i < nbinboundaries-1) {" << endl;
  outf_e << "      if (value >= bins[i] && value < bins[i+1]) {" << endl;
  outf_e << "        bin = i+1;" << endl;
  outf_e << "        break;" << endl;
  outf_e << "      }" << endl;
  outf_e << "    } else if (i == nbinboundaries-1) {" << endl;
  outf_e << "      if (value >= bins[i]) {" << endl;
  outf_e << "        bin = nbinboundaries;" << endl;
  outf_e << "        break;" << endl;
  outf_e << "      }" << endl;
  outf_e << "    }    " << endl;
  outf_e << "  }" << endl;
  outf_e << "  return bin;" << endl;

//.........这里部分代码省略.........
  }
  for (int bineta=1; bineta<=nbins_eta; bineta++){
    double sum_invwsq=0;
    double sum_prob_invwsq=0;
    for (int bin=0; bin<nbinsx; bin++){
      double bincontent = hProb[bin]->GetBinContent(bineta);
      double binerror = hProb[bin]->GetBinError(bineta);

      double invwsq = 0, prob_invwsq = 0;
      if (binerror!=0){
        invwsq = 1./pow(binerror, 2);
        prob_invwsq = bincontent*invwsq;
      }
      sum_invwsq += invwsq;
      sum_prob_invwsq += prob_invwsq;
    }

    if (sum_invwsq>0){
      sum_prob_invwsq /= sum_invwsq;
      sum_invwsq = 1./sqrt(sum_invwsq);
    }
    hProb_all->SetBinContent(bineta, sum_prob_invwsq);
    hProb_all->SetBinError(bineta, sum_invwsq);
  }
  hProb_all->Scale(1./hProb_all->Integral("width"));


  double eta_center_all[nbins_eta];
  double prob_center_all[nbins_eta];
  double eta_center_all_err[nbins_eta];
  double prob_center_all_err[nbins_eta];
  int nNonEmpty_all = 0;
  for (int bineta=1; bineta<=nbins_eta; bineta++){
    if (hProb_all->GetBinError(bineta)>0){
      eta_center_all[nNonEmpty_all] = hEta[nbinsx]->GetBinContent(bineta);
      eta_center_all_err[nNonEmpty_all] = hEta[nbinsx]->GetBinError(bineta);
      prob_center_all[nNonEmpty_all] = hProb_all->GetBinContent(bineta);
      prob_center_all_err[nNonEmpty_all] = hProb_all->GetBinError(bineta);
      cout << nNonEmpty_all << '\t' << prob_center_all[nNonEmpty_all] << '\t' << prob_center_all_err[nNonEmpty_all] << endl;
      nNonEmpty_all++;
    }
  }

  TGraphErrors* tg_all = new TGraphErrors(nNonEmpty_all, eta_center_all, prob_center_all, eta_center_all_err, prob_center_all_err);
  tg_all->SetNameTitle("tg_njets1_pjvbf_all", Form("m_{4l}: [%.0f, %.0f] GeV", bins_mzz[0], bins_mzz[nbinsx]));
  tg_all->GetXaxis()->SetTitle("#eta");
  tg_all->GetYaxis()->SetTitle("<P_{VBF}>");
  
  TF1* fit_tg_all = new TF1("fit_tg_njets1_pjvbf_all", "gaus*(1+gaus(3))", bins_eta[0], bins_eta[nbins_eta]);
  const int npars = 6;
  double mypars[npars]={ 0.155, 5.193, 1.46, 0.362, 3.524, 0.715 };
  fit_tg_all->SetParameters(mypars);
  fit_tg_all->SetParLimits(0, 0, 1);
  fit_tg_all->SetParLimits(3, 0, 1);
  fit_tg_all->SetParLimits(1, 0, 20);
  fit_tg_all->SetParLimits(4, 0, 20);
  fit_tg_all->SetParLimits(2, 0, 20);
  fit_tg_all->SetParLimits(5, 0, 20);
  tg_all->Fit(fit_tg_all);

  double* par_postfit = fit_tg_all->GetParameters();
  double* parerr_postfit = fit_tg_all->GetParErrors();
  double parIndex[npars];
  double parIndexErr[npars]={ 0 };
  for (int ip=0; ip<npars; ip++) parIndex[ip] = (double)ip;

  double integral = 0.5*par_postfit[0]*(1.+TMath::Erf(par_postfit[1]/par_postfit[2]/TMath::Sqrt(2)));
  double integral2 = 1./(2.*TMath::Sqrt(2.*TMath::Pi()*(TMath::Power(par_postfit[2], 2)+TMath::Power(par_postfit[5], 2))));
  integral2 *= par_postfit[0]*par_postfit[3]*TMath::Exp(-0.5*pow(par_postfit[1]-par_postfit[4], 2)/(TMath::Power(par_postfit[2], 2)+TMath::Power(par_postfit[5], 2)))*
    (1.+TMath::Erf((TMath::Power(par_postfit[2], 2)*par_postfit[4]+TMath::Power(par_postfit[5], 2)*par_postfit[1]) / (par_postfit[2]*par_postfit[5]*TMath::Sqrt(2.*(TMath::Power(par_postfit[2], 2)+TMath::Power(par_postfit[5], 2))))));
  integral += integral2;

  cout << integral << endl;
  par_postfit[0] /= integral;
  parerr_postfit[0] /= integral;

  TGraphErrors* tg_pars = new TGraphErrors(npars, parIndex, par_postfit, parIndexErr, parerr_postfit);
  tg_pars->SetNameTitle("njets1_pjvbf_pars", "gaus*(1+gaus(3))");
  tg_pars->GetXaxis()->SetTitle("Parameter index");
  tg_pars->GetYaxis()->SetTitle("Parameter");
  foutput->WriteTObject(tg_pars);
  delete tg_pars;

  foutput->WriteTObject(tg_all);
  delete tg_all;

  for (int bin=0; bin<nbinsx+1; bin++){
    if (bin<nbinsx) foutput->WriteTObject(hProb[bin]);
    else foutput->WriteTObject(hProb_all);
    foutput->WriteTObject(hEta[bin]);
    foutput->WriteTObject(hDist[bin]);
    delete hDist[bin];
    delete hEta[bin];
    if (bin<nbinsx) delete hProb[bin];
    else delete hProb_all;
  }
  delete tree;

  foutput->Close();
}

//.........这里部分代码省略.........
//
//fit->SetParLimits(6, 0.1,  0.2);
//fit->SetParLimits(7, 1.9,  2.05);
//fit->SetParLimits(8, 0.04, 0.08);
//graph1->Fit(fit,"R+","",1.8,5);

//fit = new TF1("f4",DGausExp,2.,5,8);
//fit->SetNpx(1000);
//fit->SetParameters(0.2,   2.25, 0.05,
//                   0.1,   2.0,  0.06,
//                   0.03,  1.8);
//
//fit->SetParLimits(0, 0.1,  10);
//fit->SetParLimits(1, 2.2,  2.4);
//fit->SetParLimits(2, 0.03, 0.052);
//
//fit->SetParLimits(3, 0.1,  1.0);
//fit->SetParLimits(4, 1.9,    2.0);
//fit->SetParLimits(5, 0.06,   0.1);
//
//fit->SetParLimits(6, 0.02,  9);
//fit->SetParLimits(7, 0,    1.8);
//
//graph1->Fit(fit,"R+","",1.75,5);
 
//fit = new TF1("f4",TGaus,2.,5,9);
//fit->SetNpx(1000);
//fit->SetParameters(0.2,   2.25, 0.05,
//                   0.1,   2.0,  0.06,
//                   0.03,  1.8,  1.0);
//
//fit->SetParLimits(0, 0.1,  10);
//fit->SetParLimits(1, 2.2,  2.4);
//fit->SetParLimits(2, 0.03, 0.052);
//
//fit->SetParLimits(3, 0.1,  1.0);
//fit->SetParLimits(4, 1.9,    2.0);
//fit->SetParLimits(5, 0.06,   0.1);
//
//fit->SetParLimits(6, 0.02,  90);
////fit->SetParLimits(7, -10,    2);
//fit->SetParLimits(7, -10,    2.4);
//fit->SetParLimits(8, 0.2,    2);
  
//fit = new TF1("f4",DGausAExp,2.,5,9);
//fit->SetNpx(1000);
//fit->SetParameters(0.5,   2.2,  30,
//                   0.1,   2.0,  0.06,
//                   0.03,  1.8,  1.0);
//
//fit->SetParLimits(0, 0.1,  100);
//fit->SetParLimits(1, 2.1,  2.3);
//fit->SetParLimits(2, 10,   60);
//
//fit->SetParLimits(3, 0.1,  1.0);
//fit->SetParLimits(4, 1.9,    2.0);
//fit->SetParLimits(5, 0.06,   0.1);
//
//fit->SetParLimits(6, 0.02,  90);
//fit->SetParLimits(7, -10,    2.4);
//fit->SetParLimits(8, 0.2,    2);
 
  fit = new TF1("f4",AGausDExp,2.,5,9);
  fit->SetNpx(1000);
  fit->SetParameters(50,   2.2,  30,
                     0.1,  2.0,  0.06,
                     0.1,  2.2,  1.8);
  
  fit->SetParLimits(0, 0.1,  100);
  fit->SetParLimits(1, 2.1,  2.3);
  fit->SetParLimits(2, 10,   60);
  
  fit->SetParLimits(3, 0.1,  1.0);
  fit->SetParLimits(4, 1.9,    2.1);
  fit->SetParLimits(5, 0.06,   0.1);
  
  fit->SetParLimits(6, 0.02,  90);
  fit->SetParLimits(7, -10,    2.4);
  fit->SetParLimits(8, 0.2,    2);
  
  graph1->Fit(fit,"R+","",1.77,5);



  for (double x=1.77; x<5-0.0000001;x+=gap){
    double xs = fit->Eval(x);
    //fitxs<<x<<"\t"<<xs<<"\t"<<0<<endl;
    fitxs << "xx.push_back(" << x << ");\tyy.push_back("<< xs << ");\t er.push_back(0);" <<endl;
  }

  gbes->Draw("P");


  file->WriteTObject(c2);
  graph1->GetXaxis()->SetRangeUser(1.6,3.5);
  graph1->GetYaxis()->SetRangeUser(0,0.5);
  c2->Print("cross.pdf");

  return 0;
}

//.........这里部分代码省略.........
      double chi2result =histos.at(i)->Chi2Test(histos.at(i+1),"WW");
      chi2result=floor(chi2result*1000+0.5)/1000;

      stringstream ss;
      //ss << "     KS: " <<std::setprecision(3) << ksresult << " chi2: " <<std::setprecision(3) << chi2result << " Private Work"; 
      ss << "            Private Work";
      const char * ch = & ss.str().c_str();;
      TLatex * ks = new TLatex(0.1, 0.9-0.03*i, ch );
      ks->SetTextColor(histos.at(i)->GetLineColor());
      ks->SetNDC();
      ks->Draw("");      

    }

    TLegend* l = new TLegend(0.40,0.9,0.69,0.99);
    // Options for legend
    l->SetBorderSize(0);
    l->SetLineStyle(0);
    l->SetTextSize(0.049);
    l->SetFillStyle(0);
    for(size_t i=0;i<names.size();i++){
      l->AddEntry(histos.at(i),names.at(i),"L");
    }
    l->Draw("same");

/*
// Clone histograms and draw ratio plot
c->cd(2);
 TH1F* ratioHisto = (TH1F*)histos.at(0)->Clone();
ratioHisto->Divide(histos.at(1));
ratioHisto->SetLineColor(kBlue);
ratioHisto->SetStats(false);
ratioHisto->GetYaxis()->SetTitle("Ratio #frac{noFxFx}{FxFx}");
// Same Size like in histogram
ratioHisto->SetLabelSize(histos.at(0)->GetLabelSize() * 0.7 / 0.3);
ratioHisto->SetTitleOffset((histos.at(0)->GetTitleOffset("Y") * 0.3 / 0.7), "Y");
ratioHisto->SetTitleSize((histos.at(0)->GetTitleSize("Y") * 0.7 / 0.3), "Y");
ratioHisto->SetTitleOffset((histos.at(0)->GetTitleOffset("X")), "X");
ratioHisto->SetTitleSize((histos.at(0)->GetTitleSize("X") * 0.7 / 0.3), "X");
// Use nicer range
ratioHisto->GetYaxis()->SetRangeUser(0, 2.2);
ratioHisto->GetYaxis()->SetNdivisions(503);
ratioHisto->GetYaxis()->SetLabelSize(0.06 * 0.7 / 0.3);
ratioHisto->Draw();*/
}
else
{
        TPad *mainPad = new TPad("","",0.0,0.0,1.0,1.0);
        mainPad->SetNumber(1);
        mainPad->SetBottomMargin(0.15);
        mainPad->SetRightMargin(0.04);
        mainPad->SetLeftMargin(0.13);
        mainPad->Draw();
        gStyle->SetOptTitle(0);
        //mainPad->SetLogx(1);
	c->cd(1);

    histos.at(0)->Draw("histo E");
    for(size_t i=0;i<histos.size();i++){
      histos.at(i)->Draw("histo same E");
    }


    for(size_t i=0;i<histos.size()-1;i=i+2){
      /*
      double ksresult = histos.at(i)->KolmogorovTest(histos.at(i+1));
      ksresult=floor(ksresult*1000+0.5)/1000;
      double chi2result =histos.at(i)->Chi2Test(histos.at(i+1),"WW");
      chi2result=floor(chi2result*1000+0.5)/1000;

      stringstream ss;
      ss << "KS: " <<std::setprecision(3) << ksresult << " chi2: " <<std::setprecision(3) << chi2result; 
      const char * ch = & ss.str().c_str();;
      TText * ks = new TText(0.1, 0.9-0.03*i, ch );
      ks->SetTextColor(histos.at(i)->GetLineColor());
      ks->SetNDC();
      ks->Draw("");      
      */
    }

    TLegend* l = new TLegend(0.65,0.5,0.9,0.7);
    l->SetBorderSize(0);
    l->SetLineStyle(0);
    //    l->SetTextSize(0.039);
    l->SetFillStyle(0);
    for(size_t i=0;i<names.size();i++){
      l->AddEntry(histos.at(i),names.at(i),"L");
    }
    l->Draw("same");
}

    c->Print("plots.pdf");
    c->SaveAs(pictureName);
    vergleich->WriteTObject(c);

}
  c->Print("plots.pdf]");


}

void buildPlotAxial_DD(){

	TH2F *dummyHist = new TH2F("dum","",1,1,900,1,6E-43,1E-37);
	dummyHist->GetYaxis()->SetTitle("#sigma_{SD} (cm^{2})");
	dummyHist->GetXaxis()->SetTitle("m_{DM} (GeV)");
	dummyHist->GetYaxis()->SetTitleOffset(1.5);
	dummyHist->GetYaxis()->SetTitleSize(0.046);
	dummyHist->GetXaxis()->SetTitleSize(0.046);
	dummyHist->GetYaxis()->SetTitleOffset(1.5);
	dummyHist->GetXaxis()->SetTitleOffset(1.25);
	dummyHist->GetYaxis()->SetTitleSize(0.045);
	dummyHist->GetXaxis()->SetTitleSize(0.045);
	dummyHist->GetYaxis()->SetLabelSize(0.04);
	dummyHist->GetXaxis()->SetLabelSize(0.04);


	gROOT->ProcessLine(".x paperStyle.C"); 
	gStyle->SetOptStat(0);
	gStyle->SetPalette(51);
	
	TFile *fiDD = TFile::Open("MassLimit_1_801_0_Both_updatePICO.root");
	TGraph *DDF = (TGraph*)fiDD->Get("DD");

	TFile *fi = TFile::Open("MassLimit_1_801_0_Both.root");

	TCanvas *can = new TCanvas();
	dummyHist->Draw("AXIS");
	
	TGraph *combinedD 	= makeOBV(combined_obs); combinedD->SetLineColor(1); combinedD->SetName("EXP");
	TGraph *combined_obsD 	= makeOBV(combined); combined_obsD->SetLineColor(1);combined_obsD->SetName("OBS");
	TGraph *monojetD 	= makeOBV(monojet); monojetD->SetLineColor(kBlue); monojetD->SetName("mjet");
	// Phil swapped boosted and resolved
	TGraph *boostedD 	= makeOBV(resolved); boostedD->SetLineColor(kMagenta); boostedD->SetName("bjet");
	TGraph *resolvedD 	= makeOBV(boosted); resolvedD->SetLineColor(kGreen); resolvedD->SetName("rjet");
	TGraph *monovD 		= makeOBV(monov); monovD->SetLineColor(kRed); monovD->SetName("MV");

	/* other lines */
	combinedD->Draw("lsame"); 
	combined_obsD->Draw("lsame"); 
	monojetD->Draw("lsame"); 
	boostedD->Draw("lsame"); 
	resolvedD->Draw("lsame"); 
	monovD->Draw("lsame"); 
	DDF->Draw("lsame");
	

	TLegend *leg = new TLegend(0.19,0.40,0.51,0.63,NULL,"brNDC");
	leg->SetFillStyle(0);
	leg->AddEntry(combinedD,"Median Expected","L");
	leg->AddEntry(combined_obsD,"Observed","L");
	leg->AddEntry(monojetD,"Monojet","L");
	leg->AddEntry(boostedD,"Boosted","L");
	leg->AddEntry(resolvedD,"Resolved","L");
	leg->AddEntry(monovD,"V-tagged","L");
	leg->AddEntry(DDF,"PICO-2L","L");
	leg->Draw();
	
   	TLatex *   texCMS = new TLatex(0.20,0.84,"#bf{CMS}");
	//TLegend *WHITEBOX = new TLegend(0.18,0.83,0.3,0.89);
   	//TLatex *   texCMS = new TLatex(0.20,0.84,"#bf{CMS}#it{Preliminary}");
	TLegend *WHITEBOX = new TLegend(0.18,0.83,0.42,0.89);
	WHITEBOX->SetFillColor(kWhite);
	//WHITEBOX->Draw();
	texCMS->SetNDC();
   	texCMS->SetTextFont(42);
   	texCMS->SetLineWidth(2);
   	texCMS->SetTextSize(0.042); texCMS->Draw();
	//tex->SetFillColor(kWhite);
	TLatex * tex = new TLatex();
	tex->SetNDC();
   	tex->SetTextFont(42);
   	tex->SetLineWidth(2);
   	tex->SetTextSize(0.035);
	tex->Draw();
   	tex->DrawLatex(0.69,0.94,"19.7 fb^{-1} (8 TeV)");
   	tex->DrawLatex(0.20,0.8,"g_{DM}=g_{SM}=1");
   	tex->DrawLatex(0.20,0.7,"Axial");
	
	can->SetRightMargin(0.11);
	can->SetLogx();
	can->SetLogy();
	can->RedrawAxis();	

	can->SaveAs("MassLimit_1_801_0_Both_DD.pdf");
	
        TFile *rout = new TFile("axialDD_out.root","RECREATE");
        rout->WriteTObject(combined_obsD,"combined");
        rout->WriteTObject(combinedD,"combined_expected"); 
        rout->WriteTObject(DDF,"DD_mass");
	rout->Close();
	
/*
 KEY: TCanvas	A;1	A
  KEY: TGraph	DD_mass;1	DD_mass
  KEY: TGraph	DD;1	DD
  KEY: TGraph	monojet;1	monojet
  KEY: TGraph	resolved;1	resolved
  KEY: TGraph	boosted;1	boosted
  KEY: TGraph	monov;1	monov
  KEY: TGraph	combined;1	combined
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    hfill[1][0]->SetLineColor(kBlack);
    hfill[1][0]->SetLineStyle(3);
    hfill[1][0]->Draw("histsame");
  }

  hfill[iDraw][0]->Draw("histsame");

  TLegendEntry* legendtext;
  if (SignalOnly==0){
    legendtext = ll->AddEntry(hfill[iDraw][0], "SM total", "f");
    legendtext = ll->AddEntry(hfill[1][0], "SM bkg.", "f");
    legendtext->SetFillStyle(1001);
    legendtext->SetFillColor(hfill[1][0]->GetFillColor());
  }
  else{
    legendtext = ll->AddEntry(hfill[iDraw][0], "SM signal", "f");
    legendtext->SetFillStyle(3001);
  }

  if (SignalOnly==0){
    legendtext = ll2->AddEntry(hfill[iDraw][4], Form("#Gamma_{H}=#Gamma^{SM}_{H}, %s", strACtitle[4].Data()), "f");
    legendtext->SetFillStyle(3001);
    legendtext->SetFillColor(hfill[iDraw][4]->GetFillColor());
    legendtext = ll2->AddEntry(hfill[iDraw][2], Form("#Gamma_{H}=#Gamma^{SM}_{H}, %s", strACtitle[2].Data()), "f");
    legendtext->SetFillStyle(3001);
    legendtext->SetFillColor(hfill[iDraw][2]->GetFillColor());
    legendtext = ll2->AddEntry(hfill[iDraw][1], Form("#Gamma_{H}=#Gamma^{SM}_{H}, %s", strACtitle[1].Data()), "f");
    legendtext->SetFillStyle(3001);
    legendtext->SetFillColor(hfill[iDraw][1]->GetFillColor());
    legendtext = ll2->AddEntry(hfill[iDraw][3], Form("#Gamma_{H}=#Gamma^{SM}_{H}, %s", strACtitle[3].Data()), "f");
    legendtext->SetFillStyle(3001);
    legendtext->SetFillColor(hfill[iDraw][3]->GetFillColor());
  }
  else{
    legendtext = ll->AddEntry(hfill[iDraw][4], strACtitle[4].Data(), "f");
    legendtext->SetFillStyle(3001);
    legendtext->SetFillColor(hfill[iDraw][4]->GetFillColor());
    legendtext = ll->AddEntry(hfill[iDraw][2], strACtitle[2].Data(), "f");
    legendtext->SetFillStyle(3001);
    legendtext->SetFillColor(hfill[iDraw][2]->GetFillColor());
    legendtext = ll->AddEntry(hfill[iDraw][1], strACtitle[1].Data(), "f");
    legendtext->SetFillStyle(3001);
    legendtext->SetFillColor(hfill[iDraw][1]->GetFillColor());
    legendtext = ll->AddEntry(hfill[iDraw][3], strACtitle[3].Data(), "f");
    legendtext->SetFillStyle(3001);
    legendtext->SetFillColor(hfill[iDraw][3]->GetFillColor());
  }
  if (SignalOnly==1){
    legendtext = ll2->AddEntry(hfill[iDraw][8], strACtitle[8].Data(), "f");
    legendtext->SetFillStyle(3001);
    legendtext->SetFillColor(hfill[iDraw][8]->GetFillColor());
    legendtext = ll2->AddEntry(hfill[iDraw][5], strACtitle[5].Data(), "f");
    legendtext->SetFillStyle(3001);
    legendtext->SetFillColor(hfill[iDraw][5]->GetFillColor());
    legendtext = ll2->AddEntry(hfill[iDraw][7], strACtitle[7].Data(), "f");
    legendtext->SetFillStyle(3001);
    legendtext->SetFillColor(hfill[iDraw][7]->GetFillColor());
    legendtext = ll2->AddEntry(hfill[iDraw][6], strACtitle[6].Data(), "f");
    legendtext->SetFillStyle(3001);
    legendtext->SetFillColor(hfill[iDraw][6]->GetFillColor());
  }

  ll->Draw("same");
  ll2->Draw("same");
  ptx->Draw();
  pt->Draw();
  cc->RedrawAxis();
  cc->Update();

  canvasname.Prepend(coutput_common);
  TString canvasname_pdf = canvasname;
  TString canvasname_eps = canvasname;
  TString canvasname_png = canvasname;
  TString canvasname_root = canvasname;
  TString canvasname_c = canvasname;
  canvasname_pdf.Append(".pdf");
  canvasname_eps.Append(".eps");
  canvasname_png.Append(".png");
  canvasname_root.Append(".root");
  canvasname_c.Append(".C");
  cc->SaveAs(canvasname_pdf);
  cc->SaveAs(canvasname_eps);
  cc->SaveAs(canvasname_png);
  cc->SaveAs(canvasname_root);
  cc->SaveAs(canvasname_c);

  foutput->WriteTObject(cc);
  delete ll2;
  delete ll;
  cc->Close();
  delete ptx;
  delete pt;
  for (int t=0; t<4; t++){
    for (int ac=0; ac<5; ac++){
      foutput->WriteTObject(hfill[t][ac]);
      delete hfill[t][ac];
    }
  }
  foutput->Close();
}

//.........这里部分代码省略.........
        mela.setMelaHiggsWidth(wPOLE);
        p0_g1prime2_VAJHU = getJHUGenMELAWeight(mela, lepIdOrdered, angularOrdered, selfDHvvcoupl);

        for (int ii = 0; ii < SIZE_HVV; ii++){ for (int jj = 0; jj < 2; jj++) selfDHvvcoupl[ii][jj] = 0; }
        selfDHvvcoupl[13][0]=-pow(10000.0/mPOLE, 2);
        mela.setMelaLeptonInterference(TVar::InterfOn);
        mela.setMelaHiggsWidth(wPOLE);
        p0_g1prime4_VAJHU = getJHUGenMELAWeight(mela, lepIdOrdered, angularOrdered, selfDHvvcoupl);


        for (int ii = 0; ii < SIZE_HVV; ii++){ for (int jj = 0; jj < 2; jj++) selfDHvvcoupl[ii][jj] = 0; }
        selfDHvvcoupl[0][0]=1;
        selfDHvvcoupl[1][0]=1.638;
        mela.setMelaLeptonInterference(TVar::InterfOn);
        mela.setMelaHiggsWidth(wPOLE);
        pg1g2_VAJHU = getJHUGenMELAWeight(mela, lepIdOrdered, angularOrdered, selfDHvvcoupl);

        for (int ii = 0; ii < SIZE_HVV; ii++){ for (int jj = 0; jj < 2; jj++) selfDHvvcoupl[ii][jj] = 0; }
        selfDHvvcoupl[0][0]=1;
        selfDHvvcoupl[3][0]=2.521;
        mela.setMelaLeptonInterference(TVar::InterfOn);
        mela.setMelaHiggsWidth(wPOLE);
        pg1g4_VAJHU = getJHUGenMELAWeight(mela, lepIdOrdered, angularOrdered, selfDHvvcoupl);

        for (int ii = 0; ii < SIZE_HVV; ii++){ for (int jj = 0; jj < 2; jj++) selfDHvvcoupl[ii][jj] = 0; }
        selfDHvvcoupl[0][0]=1;
        selfDHvvcoupl[11][0]=-12046.01;
        mela.setMelaLeptonInterference(TVar::InterfOn);
        mela.setMelaHiggsWidth(wPOLE);
        pg1g1prime2_VAJHU = getJHUGenMELAWeight(mela, lepIdOrdered, angularOrdered, selfDHvvcoupl);

        for (int ii = 0; ii < SIZE_HVV; ii++){ for (int jj = 0; jj < 2; jj++) selfDHvvcoupl[ii][jj] = 0; }
        selfDHvvcoupl[0][0]=1;
        selfDHvvcoupl[13][0]=-pow(10000.0/mPOLE, 2);
        mela.setMelaLeptonInterference(TVar::InterfOn);
        mela.setMelaHiggsWidth(wPOLE);
        pg1g1prime4_VAJHU = getJHUGenMELAWeight(mela, lepIdOrdered, angularOrdered, selfDHvvcoupl);

        if (p0plus_VAJHU!=0){
          p0hplus_VAJHU/=p0plus_VAJHU;
          p0minus_VAJHU/=p0plus_VAJHU;
          p0_g1prime2_VAJHU/=p0plus_VAJHU;
          p0_g1prime4_VAJHU/=p0plus_VAJHU;

          pg1g2_VAJHU/=p0plus_VAJHU;
          pg1g4_VAJHU/=p0plus_VAJHU;
          pg1g1prime2_VAJHU/=p0plus_VAJHU;
          pg1g1prime4_VAJHU/=p0plus_VAJHU;
          p0plus_VAJHU=1;
        }
        else{
          p0plus_VAJHU=0;
          p0hplus_VAJHU=0;
          p0minus_VAJHU=0;
          p0_g1prime2_VAJHU=0;
          p0_g1prime4_VAJHU=0;

          pg1g2_VAJHU=0;
          pg1g4_VAJHU=0;
          pg1g1prime2_VAJHU=0;
          pg1g1prime4_VAJHU=0;
        }
      }
      else{
        p0plus_VAJHU=1;
        p0hplus_VAJHU=1;
        p0minus_VAJHU=1;
        p0_g1prime2_VAJHU=1;
        p0_g1prime4_VAJHU=1;

        pg1g2_VAJHU=1;
        pg1g4_VAJHU=1;
        pg1g1prime2_VAJHU=1;
        pg1g1prime4_VAJHU=1;
      }

      p0plus_VAJHU*=templateWeight;
      p0hplus_VAJHU*=templateWeight;
      p0minus_VAJHU*=templateWeight;
      p0_g1prime2_VAJHU*=templateWeight;
      p0_g1prime4_VAJHU*=templateWeight;

      pg1g2_VAJHU*=templateWeight;
      pg1g4_VAJHU*=templateWeight;
      pg1g1prime2_VAJHU*=templateWeight;
      pg1g1prime4_VAJHU*=templateWeight;

      templateTree->Fill();
    }

    foutput->WriteTObject(templateTree);
    delete tree;
    delete templateTree;
    foutput->Close();
  }

  delete tg_interf;
  delete tgkf;
  finput_KDFactor->Close();
}

//.........这里部分代码省略.........
                                                    RooFit::Import("Pass",*pseudoData_pass),
                                                    RooFit::Import("Fail",*pseudoData_fail));

    pseudoDataCombined->write(Form("toy%d.txt",t));

    RooFitResult *fitResult=0;
    fitResult = totalPdf->fitTo(*pseudoDataCombined,
                                RooFit::Extended(),
                                RooFit::Strategy(2),
                                //RooFit::Minos(RooArgSet(eff)),
                                RooFit::Save());

    cout << "\n\n";
    cout << "Eff Fit: " << eff->getVal() << " -" << fabs(eff->getErrorLo()) << " +" << eff->getErrorHi() << endl;

    //Fill Tree
    varEff = eff->getVal();
    varEffErrL = fabs(eff->getErrorLo());
    varEffErrH = eff->getErrorHi();
    outTree->Fill();


//   //*************************************************************************************
//   //Plot Toys
//   //*************************************************************************************
//   TCanvas *cv = new TCanvas("cv","cv",800,600);
//   char pname[50];
//   char binlabelx[100];
//   char binlabely[100];
//   char yield[50];
//   char effstr[100];
//   char nsigstr[100];
//   char nbkgstr[100];
//   char chi2str[100];

//   //
//   // Plot passing probes
//   //

//   RooPlot *mframeFail_default = m.frame(Bins(Int_t(130-60)/2));
//   modelFail_default->plotOn(mframeFail_default);
//   modelFail_default->plotOn(mframeFail_default,Components("ebackgroundFail"),LineStyle(kDashed),LineColor(kRed));
//   modelFail_default->plotOn(mframeFail_default,Components("bkgexpFail"),LineStyle(kDashed),LineColor(kGreen+2));
//   mframeFail_default->Draw();
//   cv->SaveAs("DefaultModel.gif");





//   RooPlot *mframeFail = m.frame(Bins(Int_t(130-60)/2));
//   modelFail->plotOn(mframeFail);
//   modelFail->plotOn(mframeFail,Components("ebackgroundFail_shifted"),LineStyle(kDashed),LineColor(kRed));
//   modelFail->plotOn(mframeFail,Components("bkgexpFail_shifted"),LineStyle(kDashed),LineColor(kGreen+2));

//   sprintf(yield,"%u Events",(Int_t)passTree->GetEntries());
//   sprintf(nsigstr,"N_{sig} = %.1f #pm %.1f",NsigPass.getVal(),NsigPass.getPropagatedError(*fitResult));
//     plotPass.AddTextBox(yield,0.21,0.76,0.51,0.80,0,kBlack,-1);    
//   plotPass.AddTextBox(effstr,0.70,0.85,0.94,0.90,0,kBlack,-1);
//     plotPass.AddTextBox(0.70,0.73,0.94,0.83,0,kBlack,-1,1,nsigstr);//,chi2str);

//   mframeFail->Draw();
//   cv->SaveAs(Form("ShiftedModel_%d.gif",Option));


 
//   //
//   // Plot failing probes
//   //
//   sprintf(pname,"fail%s_%i",name.Data(),ibin);
//   sprintf(yield,"%u Events",(Int_t)failTree->GetEntries());
//   sprintf(nsigstr,"N_{sig} = %.1f #pm %.1f",NsigFail.getVal(),NsigFail.getPropagatedError(*fitResult));
//   sprintf(nbkgstr,"N_{bkg} = %.1f #pm %.1f",NbkgFail.getVal(),NbkgFail.getPropagatedError(*fitResult));
//   sprintf(chi2str,"#chi^{2}/DOF = %.3f",mframePass->chiSquare(nflfail));
//   CPlot plotFail(pname,mframeFail,"Failing probes","tag-probe mass [GeV/c^{2}]","Events / 2.0 GeV/c^{2}");
//   plotFail.AddTextBox(binlabelx,0.21,0.85,0.51,0.90,0,kBlack,-1);
//   if((name.CompareTo("etapt")==0) || (name.CompareTo("etaphi")==0)) {
//     plotFail.AddTextBox(binlabely,0.21,0.80,0.51,0.85,0,kBlack,-1);    
//     plotFail.AddTextBox(yield,0.21,0.76,0.51,0.80,0,kBlack,-1);    
//   } else {
//     plotFail.AddTextBox(yield,0.21,0.81,0.51,0.85,0,kBlack,-1);
//   }
//   plotFail.AddTextBox(effstr,0.70,0.85,0.94,0.90,0,kBlack,-1);  
//   plotFail.AddTextBox(0.70,0.68,0.94,0.83,0,kBlack,-1,2,nsigstr,nbkgstr);//,chi2str);
//   plotFail.Draw(cfail,kTRUE,format);  




  } //for loop over all toys
  



  //*************************************************************************************
  //Save To File
  //*************************************************************************************
  outputfile->WriteTObject(outTree, outTree->GetName(), "WriteDelete");

}

//.........这里部分代码省略.........
  t1.Draw();
  t2.Draw();
  t3.Draw();
  cv->SaveAs(Form("Unrolled%s.pdf",Label.c_str()));

  // Unrolled plots in percentages
  cv = new TCanvas("cv","cv", 800,600);
  legend = new TLegend(0.85,0.20,0.95,0.80);
  legend->SetTextSize(0.03);
  legend->SetBorderSize(0);

  for (Int_t i = histMRRsq.size()-1 ; i >= 0; --i) {
    if (hasSignal && i==0) {
      legend->AddEntry(histMRRsq[i],processLabels[i].c_str(), "L");
    } else {
      legend->AddEntry(histMRRsq[i],processLabels[i].c_str(), "F");
    }
  }
  stackUnrolledPercentage->Draw();
  stackUnrolledPercentage->GetHistogram()->GetXaxis()->SetTitle(((TH1F*)(stackUnrolledPercentage->GetHists()->At(0)))->GetXaxis()->GetTitle());
  // stackUnrolledPercentage->GetHistogram()->GetXaxis()->SetRangeUser(0, 35);
  stackUnrolledPercentage->GetHistogram()->GetYaxis()->SetTitle(((TH1F*)(stackUnrolledPercentage->GetHists()->At(0)))->GetYaxis()->GetTitle());
  if(hasSignal) histUnrolledPercentage[0]->Draw("same hist"); 
  legend->Draw();
  t1.Draw();
  t2.Draw();
  t3.Draw();
  cv->SaveAs(Form("UnrolledPercentage%s.pdf",Label.c_str()));

  // Unrolled plots in sideband vs signal box
  cv = new TCanvas("cv","cv", 800,600);
  legend = new TLegend(0.85,0.20,0.95,0.80);
  legend->SetTextSize(0.03);
  legend->SetBorderSize(0);

  for (Int_t i = histMRRsq.size()-1 ; i >= 0; --i) {
    if (hasSignal && i==0) {
      legend->AddEntry(histMRRsq[i],processLabels[i].c_str(), "L");
    } else {
      legend->AddEntry(histMRRsq[i],processLabels[i].c_str(), "F");
    }
  }
  stackUnrolled2bins->Draw();
  stackUnrolled2bins->GetHistogram()->GetXaxis()->SetTitle(((TH1F*)(stackUnrolled2bins->GetHists()->At(0)))->GetXaxis()->GetTitle());
  stackUnrolled2bins->GetHistogram()->GetYaxis()->SetTitle(((TH1F*)(stackUnrolled2bins->GetHists()->At(0)))->GetYaxis()->GetTitle());
  if(hasSignal) histUnrolled2bins[0]->Draw("same hist");
  legend->Draw();
  t1.Draw();
  t2.Draw();
  t3.Draw();
  cv->SaveAs(Form("Unrolled2bins%s.pdf",Label.c_str()));

  // Unrolled plots in sideband vs signal box in fractions
  cv = new TCanvas("cv","cv", 800,600);
  legend = new TLegend(0.7,0.23,0.90,0.88);
  legend->SetTextSize(0.03);
  legend->SetBorderSize(0);

  for (Int_t i = histMRRsq.size()-1 ; i >= 0; --i) {
    if (hasSignal && i==0) {
      legend->AddEntry(histMRRsq[i],processLabels[i].c_str(), "L");
    } else {
      legend->AddEntry(histMRRsq[i],processLabels[i].c_str(), "F");
    }
  }

  stackUnrolledPercentage2bins->Draw();
  stackUnrolledPercentage2bins->GetHistogram()->GetXaxis()->SetTitle(((TH1F*)(stackUnrolledPercentage2bins->GetHists()->At(0)))->GetXaxis()->GetTitle());
  stackUnrolledPercentage2bins->GetHistogram()->GetYaxis()->SetTitle(((TH1F*)(stackUnrolledPercentage2bins->GetHists()->At(0)))->GetYaxis()->GetTitle());
  stackUnrolledPercentage2bins->GetHistogram()->GetXaxis()->SetBinLabel(1, "Sideband");
  stackUnrolledPercentage2bins->GetHistogram()->GetXaxis()->SetBinLabel(2, "Signal Sensitive Region");
  if(hasSignal) histUnrolledPercentage2bins[0]->Draw("same hist");
  legend->Draw();
  t1.Draw();
  t2.Draw();
  t3.Draw();
  cv->SaveAs(Form("UnrolledPercentage2bins%s.pdf",Label.c_str()));


   //--------------------------------------------------------------------------------------------------------------
  // Output
  //==============================================================================================================
  TFile *file = TFile::Open(("RazorPlots"+Label+".root").c_str(), "RECREATE");
  file->cd();

  for(int i=0; i<int(inputfiles.size()); i++) {
    file->WriteTObject(histMR[i], Form("histMR_%s",processLabels[i].c_str()), "WriteDelete");
    file->WriteTObject(histRsq[i], Form("histRsq_%s",processLabels[i].c_str()), "WriteDelete");
    file->WriteTObject(histMRRsq[i], Form("histMRRsq_%s",processLabels[i].c_str()), "WriteDelete");
    histUnrolled[i]->Write();  
    histUnrolled2bins[i]->Write();  
    histUnrolledPercentage[i]->Write();  
    histUnrolledPercentage2bins[i]->Write();  
  }
  
  stackUnrolled->Write();
  stackUnrolled2bins->Write();
  stackUnrolledPercentage->Write();
  stackUnrolledPercentage2bins->Write();
 }