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

//**********MAIN**************************************************************************
int main(int argc, char* argv[])
{
    if(argc < 2)
    {
        cout << argv[0] << " cfg file " << "[run]" << endl; 
        return -1;
    }

    //---memory consumption tracking---
    float cpu[2]{0}, mem[2]={0}, vsz[2]={0}, rss[2]={0};

    //---load options---    
    CfgManager opts;
    opts.ParseConfigFile(argv[1]);

    //-----input setup-----    
    if(argc > 2)
    {
        vector<string> run(1, argv[2]);
        opts.SetOpt("h4reco.run", run);
    }
    string outSuffix = opts.GetOpt<string>("h4reco.outNameSuffix");
    string run = opts.GetOpt<string>("h4reco.run");
    TChain* inTree = new TChain("H4tree");
    ReadInputFiles(opts, inTree);
    H4Tree h4Tree(inTree);

    //-----output setup-----
    uint64 index=stoul(run)*1e9;
    TFile* outROOT = new TFile(outSuffix+TString(run)+".root", "RECREATE");
    outROOT->cd();
    RecoTree mainTree(&index);

    //---Get plugin sequence---
    PluginLoader<PluginBase>* loader;
    vector<PluginLoader<PluginBase>* > pluginLoaders;    
    map<string, PluginBase*> pluginMap;
    vector<PluginBase*> pluginSequence;
    vector<string> pluginList = opts.GetOpt<vector<string> >("h4reco.pluginList");    
    //---plugin creation
    pluginLoaders.reserve(pluginList.size());
    for(auto& plugin : pluginList)
    {
        cout << ">>> Loading plugin <" << plugin << ">" << endl;
        //---create loader 
        loader = new PluginLoader<PluginBase>(opts.GetOpt<string>(plugin+".pluginType"));
        pluginLoaders.push_back(loader);
        pluginLoaders.back()->Create();
        //---get instance and put it in the plugin sequence   
        PluginBase* newPlugin = pluginLoaders.back()->CreateInstance();
        if(newPlugin)
        {
            pluginSequence.push_back(newPlugin);
            pluginSequence.back()->SetInstanceName(plugin);
            pluginMap[plugin] = pluginSequence.back();
        }
        else
        {
            cout << ">>> ERROR: plugin type " << opts.GetOpt<string>(plugin+".pluginType") << " is not defined." << endl;
            return 0;
        }
    }

    //---begin
    for(auto& plugin : pluginSequence)
    {
        //---call Begin() methods and check the return status
        bool r_status = plugin->Begin(opts, &index);
        if(!r_status)
        {
            cout << ">>> ERROR: plugin returned bad flag from Begin() call: " << plugin->GetInstanceName() << endl;
            exit(-1);
        }
        //---Get plugin shared data
        for(auto& shared : plugin->GetSharedData("", "TTree", true))
        {
            TTree* tree = (TTree*)shared.obj;
            tree->SetMaxVirtualSize(10000);
            tree->SetDirectory(outROOT);
        }
    }
            
    //---events loop
    int maxEvents = opts.OptExist("h4reco.maxEvents") ? opts.GetOpt<int>("h4reco.maxEvents") : -1;
    cout << ">>> Processing H4DAQ run #" << run << " <<<" << endl;
    while(h4Tree.NextEntry() && (index-stoul(run)*1e9<maxEvents || maxEvents==-1))
    {
        if(index % 1000 == 0)
        {
            cout << ">>> Processed events: " << index-stoul(run)*1e9 << "/"
                 << (maxEvents<0 ? h4Tree.GetEntries() : min(h4Tree.GetEntries(), (uint64)maxEvents))
                 << endl;
            TrackProcess(cpu, mem, vsz, rss);
        }

        //---call ProcessEvent for each plugin and check the return status
        bool status=true;
        for(auto& plugin : pluginSequence)
            if(status)
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      TF1*  frelrsp = (TF1*)hrsp[i][j]->GetListOfFunctions()->Last();
      peak    = (frelrsp==0) ? hrsp[i][j]->GetMean()      :  frelrsp->GetParameter(1);
      epeak   = (frelrsp==0) ? hrsp[i][j]->GetMeanError() :  frelrsp->GetParError(1);

      //cout << " peak : " << peak << "  mean : " << hrsp[i][j]->GetMean()  << "  " << frelrsp->GetParameter(1) << endl;

      double cor  = 1.0/peak;
      double ecor = cor*cor*epeak;

      //if(i==0)cout <<ptbins[j] << " ratio of  " << refpt/jetpt << " " << cor << "  " <<  ftemp->Eval(jetpt) << " " << 1/ ftemp->Eval(jetpt) << endl;
      grsp[i]->SetPoint     (j, refpt,  peak);
      grsp[i]->SetPointError(j, erefpt, epeak);
      gcor[i]->SetPoint     (j, jetpt,  cor);
      gcor[i]->SetPointError(j, ejetpt, ecor);
    }
  }
  if( gPad )gPad->Close();
  //hrsp[0][0]->Draw("p");
  //return 0;
  
  int ipad=0;
  //! 0 - 20 GeV
  TCanvas *c98[knj];
  for(int nj=njmin;nj<njmax;nj++){
    c98[nj] = new TCanvas(Form("c99_%d",nj),Form("Fine %s Fitting plots",calgo[nj]),100,102,1399,942);
    c98[nj]->Divide(6,3,0,0);
    //c98[nj]->Divide(10,5,0,0);
    //c98[nj]->Divide(6,5,0,0);
    //c98[nj]->Divide(7,4,0,0);
    //c98[nj]->Divide(6,4,0,0);
    ipad=0;
    for(int ip=0;ip<nbins;ip++){      
      if(ptbins[ip]<minx)continue;
      c98[nj]->cd(++ipad);
      gPad->SetLeftMargin(0.15);
      gPad->SetRightMargin(0.01);
      gPad->SetBottomMargin(0.15);
      gPad->SetLogy();
      //if(ipad%10==0)gPad->SetRightMargin(0.02);
      //if(ipad%(ip+1)==0)gPad->SetLeftMargin(0.15);
    
      //hrsp[nj][ip]->SetMaximum(hrsp[nj][ip]->GetMaximum() + 0.25*hrsp[nj][ip]->GetMaximum());
      hrsp[nj][ip]->SetMaximum(hrsp[nj][ip]->Integral()*2e-00);
      hrsp[nj][ip]->SetMinimum(hrsp[nj][ip]->Integral()*1e-07);
      hrsp[nj][ip]->SetTitle(0);
      hrsp[nj][ip]->GetXaxis()->SetTitle("<reco jet p_{T} / gen jet p_{T}>");
      hrsp[nj][ip]->GetXaxis()->SetTitleFont(42);
      hrsp[nj][ip]->GetXaxis()->SetLabelFont(42);
      hrsp[nj][ip]->GetXaxis()->SetLabelSize(0.08);
      hrsp[nj][ip]->GetXaxis()->SetTitleSize(0.07);
      hrsp[nj][ip]->GetXaxis()->SetNdivisions(507);
      hrsp[nj][ip]->GetYaxis()->SetTitle("");
      hrsp[nj][ip]->GetYaxis()->SetTitleFont(42);
      hrsp[nj][ip]->GetYaxis()->SetLabelFont(42);
      hrsp[nj][ip]->GetYaxis()->SetLabelSize(0.08);
      hrsp[nj][ip]->GetYaxis()->SetNdivisions(507);
      
      hrsp[nj][ip]->SetMarkerStyle(20);
      hrsp[nj][ip]->SetMarkerColor(1);
      hrsp[nj][ip]->SetLineColor(1);
      hrsp[nj][ip]->SetMarkerSize(1.1);
      hrsp[nj][ip]->Draw("hist");  
      
      TF1 *fdscb = (TF1*)hrsp[nj][ip]->GetFunction("fdscb");
      if(fdscb){
	fdscb->SetLineWidth(1);

//.........这里部分代码省略.........
            std::string in(*(argv + f));
            inputs[inputs.size() - 1]->SetCont(atof(in.c_str()));
            cout<<"cont: "<<inputs[inputs.size() - 1]->GetCont()<<endl;

        } else if (arg_fth == "age") {
            f++;
            std::string in(*(argv + f));
            inputs[inputs.size() - 1]->SetAge(atof(in.c_str()));
            cout<<"age: "<<inputs[inputs.size() - 1]->GetAge()<<endl;

        } else if (arg_fth == "group") {
            f++;
            std::string in(*(argv + f));
            inputs[inputs.size() - 1]->SetGroup(in);
        }
    }
    int uSum = 0;
    int lSum = 0;
    int nStaff = 0;
    int nStudent = 0;
    int nPostDoc = 0;
    TH1D * price = new TH1D("price", "price", 120, 0, 600);
    price->GetXaxis()->SetTitle("maximum payement");

    TH2D * price_feed = new TH2D("price_feed", "price_feed", 120, 0, 600, 4, -1.5, 2.5);
    price_feed->GetYaxis()->SetBinLabel(1, "negative");
    price_feed->GetYaxis()->SetBinLabel(2, "positive");
    price_feed->GetYaxis()->SetBinLabel(3, "v. positive");
    price_feed->GetYaxis()->SetBinLabel(4, "occ. user");
    price_feed->GetXaxis()->SetTitle("maximum payement");

    TH2D * price_position = new TH2D("price_position", "price_position", 120, 0, 600, 3, 0, 3);
    price_position->GetYaxis()->SetBinLabel(1, "Staff");
    price_position->GetYaxis()->SetBinLabel(2, "Postdoc");
    price_position->GetYaxis()->SetBinLabel(3, "Student");
    price_position->GetXaxis()->SetTitle("maximum payement");

    TH2D * price_age = new TH2D("price_age", "price_age", 120, 0, 600, 2, -0.5, 1.5);
    price_age->GetYaxis()->SetBinLabel(1, "[30,35]");
    price_age->GetYaxis()->SetBinLabel(2, "> 35");
    price_age->GetXaxis()->SetTitle("maximum payement");

    TH2D * cont_feed = new TH2D("cont_feed", "cont_feed", 4, 0, 4, 4, -1.5, 2.5);
    cont_feed->GetXaxis()->SetBinLabel(1, "no contribution");
    cont_feed->GetXaxis()->SetBinLabel(1, "Coffee");
    cont_feed->GetXaxis()->SetBinLabel(3, "System");
    cont_feed->GetXaxis()->SetBinLabel(4, "Both");
    cont_feed->GetYaxis()->SetBinLabel(1, "negative");
    cont_feed->GetYaxis()->SetBinLabel(2, "positive");
    cont_feed->GetYaxis()->SetBinLabel(3, "v. positive");
    cont_feed->GetYaxis()->SetBinLabel(4, "occ. user");

    TH2D * feed_age = new TH2D("feed_age", "feed_age", 4, -1.5, 2.5, 2, -0.5, 1.5);
    feed_age->GetXaxis()->SetBinLabel(1, "negative");
    feed_age->GetXaxis()->SetBinLabel(2, "positive");
    feed_age->GetXaxis()->SetBinLabel(3, "v. positive");
    feed_age->GetXaxis()->SetBinLabel(4, "occ. user");
    feed_age->GetYaxis()->SetBinLabel(1, "[30,35]");
    feed_age->GetYaxis()->SetBinLabel(2, "> 35");
    for (int i = 0; i < inputs.size(); i++) {
        uSum += inputs.at(i)->GetUpLim();
        lSum += inputs.at(i)->GetLowLim();
        price->Fill(inputs.at(i)->GetUpLim());
        price_age->Fill(inputs.at(i)->GetUpLim(), inputs.at(i)->GetAge());
        price_feed->Fill(inputs.at(i)->GetUpLim(), inputs.at(i)->GetFeed());
        cont_feed->Fill(inputs.at(i)->GetCont(), inputs.at(i)->GetFeed());
        feed_age->Fill(inputs.at(i)->GetFeed(), inputs.at(i)->GetAge());
        if (inputs.at(i)->GetPosition() == "Staff") {
            price_position->Fill(inputs.at(i)->GetUpLim(), 0);
            nStaff += 1;
        } else if (inputs.at(i)->GetPosition() == "PostDoc") {
            price_position->Fill(inputs.at(i)->GetUpLim(), 1);
            nPostDoc += 1;
        } else if (inputs.at(i)->GetPosition() == "Student") {
            price_position->Fill(inputs.at(i)->GetUpLim(), 2);
            nStudent += 1;
        }
    }
    cout << "Number of participants: " << inputs.size() << endl;
    cout << "Staffs: " << nStaff << ", PostDocs: " << nPostDoc << ", Students: " << nStudent << endl;
    cout << "Based on the current numbers we have:" << endl;
    cout << "\tMaximum expected contribution in IRR:" << uSum * 10 << endl;
    cout << "\tMinimum expected contribution in IRR:" << lSum * 10 << endl;
    cout << "Note that the maximum value could be enhanced since some of the\n"
            << "contributors did not specify the value or gave only the lower\nlimit of their contribution" << endl;

    TFile * fout = new TFile("coffee.root", "recreate");
    fout->cd();
    price->Write();
    price_age->Write();
    price_feed->Write();
    price_position->Write();
    cont_feed->Write();
    feed_age->Write();
    fout->Write();
    fout->Close();


    return 0;
}

void R_scin_eff(){

  TString gDir = "/w/halla-2/e08014/Rootfiles/";
  Int_t  run_number;    TString filename;

  ofstream myfile; myfile.open("R_Scin_Eff.txt");
  ifstream runfile("../../runsample.txt"); const int NUMR = 1000;
 
  TString GeneralCut = "abs(R.tr.x)<0.75 && abs(R.tr.y)<0.55 && abs(R.tr.th)<0.15 && abs(R.tr.ph)<0.045";
  TString TriggerCut1 = "((DBB.evtypebits>>1)&1)";
  TString TriggerCut2 = "((DBB.evtypebits>>2)&1)";
  TString TriggerCut6 = "((DBB.evtypebits>>6)&1)";
  
  TFile *geocut = new TFile(Form("R_Geo_Cut_%d.root",run_number), "UPDATE");
  
  TCutG* phi_y_cut = (TCutG*) gROOT->FindObject("phi_y_cut"); //looking for old cut definition
  TCutG* theta_dp_cut = (TCutG*) gROOT->FindObject("theta_dp_cut"); //looking for old cut definition
  
  if(!phi_y_cut || !theta_dp_cut){
    
    TCanvas *c1= new TCanvas("c1","c1",1000,600);
    c1->Divide(2,1); 
    TPad *bpad = gPad;

    //    runfile >> run_number;
    TString file1 = gDir + Form("e08014_less_%d.root",3590);
    TChain *T1 = new TChain("T");
    T1->Add(file1);

    cerr <<"========= Create Acceptance CutG ..."<<endl;
    c1->cd(1);
    T1->Draw("R.tr.tg_ph:R.tr.tg_y>>gc1(300,-0.05,0.05,300,-0.05,0.05)",GeneralCut,"colz");
    gc1->SetXTitle("R.tr.tg_y"); gc1->SetYTitle("R.tr.tg_ph");  c1->Update();
    phi_y_cut = (TCutG*) (bpad->WaitPrimitive("CUTG", "CutG"));c1->Update();
    phi_y_cut->SetName("phi_y_cut"); phi_y_cut->SetVarX("R.tr.tg_y"); phi_y_cut->SetVarY("R.tr.tg_ph");
    cerr <<"Phi vs Y Done!"<<endl;

    c1->cd(2);
    T1->Draw("R.tr.tg_th:R.tr.tg_dp>>gc2(300,-0.05,0.05,300,-0.15,0.15)",GeneralCut,"colz");
    gc2->SetXTitle("R.tr.tg_dp"); gc2->SetYTitle("R.tr.tg_th");  c1->Update()
    theta_dp_cut = (TCutG*) (bpad->WaitPrimitive("CUTG", "CutG"));c1->Update();
    theta_dp_cut->SetName("theta_dp_cut"); theta_dp_cut->SetVarX("R.tr.tg_dp"); theta_dp_cut->SetVarY("R.tr.tg_th");
    cerr <<"Theta vs Dp Done!"<<endl;

    geocut->cd();
    phi_y_cut->Write("", TObject::kOverwrite); // Overwrite old cut
    theta_dp_cut->Write("", TObject::kOverwrite); // Overwrite old cut

    T1->Delete();
  }
  
  TCutG* electron_cut = (TCutG*) gROOT->FindObject("electron_cut"); //looking for old cut definition
   
  if(!electron_cut ){
    
    //    runfile >> run_number;
    TString file2 = gDir + Form("e08014_less_%d.root",3590);
    TChain *T2 = new TChain("T");
    T2->Add(file2);

    TCanvas *c2= new TCanvas("c2","c2",800,800);
    //  c2->Divide(2,1); 
    TPad *cpad = gPad;

    cerr <<"========= Create Electron and Pion CutG On Calorimeter ..."<<endl;
    c2->cd(1);
    T2->Draw("R.sh.asum_c:R.ps.asum_c>>gc3(300,-100.,6000.0,300,-100.0,6000.0)",GeneralCut,"colz");c2->Update();
    electron_cut = (TCutG*) (cpad->WaitPrimitive("CUTG", "CutG"));c2->Update();
    electron_cut->SetName("electron_cut"); electron_cut->SetVarX("R.ps.asum_c"); electron_cut->SetVarY("R.sh.asum_c");
    cerr <<"Electron Done!"<<endl;

     geocut->cd();
     electron_cut->Write("", TObject::kOverwrite); // Overwrite old cut

     T2->Delete();
  }
 
   TString cut_t1 = TriggerCut1+"&&" + GeneralCut + "&&phi_y_cut && theta_dp_cut "+ "&&electron_cut";
   TString cut_t2 = TriggerCut2+"&&" + GeneralCut + "&&phi_y_cut && theta_dp_cut "+ "&&electron_cut";  
   TString cut_t6 = TriggerCut6+"&&" + GeneralCut + "&&phi_y_cut && theta_dp_cut "+ "&&electron_cut";
  
  // TString cut_t1 = TriggerCut1+"&&" + GeneralCut +  "&&electron_cut";
  // TString cut_t2 = TriggerCut2+"&&" + GeneralCut +  "&&electron_cut";
  // TString cut_t6 = TriggerCut6+"&&" + GeneralCut +  "&&electron_cut";
  
//   TString cut_t1 = TriggerCut1+"&&" + GeneralCut;
//   TString cut_t2 = TriggerCut2+"&&" + GeneralCut;
//   TString cut_t6 = TriggerCut6+"&&" + GeneralCut;

  float rnum[NUMR]; 
  float T1eff[NUMR], T6eff[NUMR];

  int jj;  float evt,evt0,evt1,evt2,evt6;
  int ps1=1,ps2=1, ps3=1,ps4=1,ps6=1,ps7=1;

  while(!(runfile.eof())){
    runfile >> run_number; rnum[jj] = run_number;
    filename = gDir + Form("e08014_less_%d.root",run_number);
    //   TChain *T = new TChain("T");
    //   T->Add(filename);
//.........这里部分代码省略.........

void asymptoticLimit(const char* outputfile, const char* inputfiles, unsigned int verbosity=2, bool upper_exclusion=true)
{
  // prepare input names from inputfiles
  std::vector<std::string> tanb_files;
  string2Vector(cleanupWhitespaces(inputfiles), tanb_files);
  // prepare tanb values from inputfiles
  std::map<double, std::string> tanb_values;  
  for(std::vector<std::string>::const_iterator tanb_file = tanb_files.begin(); tanb_file!=tanb_files.end(); ++tanb_file){
    tanb_values[atof(tanb_file->substr(tanb_file->rfind("_")+1, tanb_file->find(".root")-tanb_file->rfind("_")-1).c_str())] = *tanb_file;
  }
  
  double limit, lowlimit;
  TGraph* graph = new TGraph();
  TFile* file = TFile::Open(outputfile, "update");
  TTree* tree = new TTree("limit", "limit");
  tree->Branch("limit", &limit, "limit/D");
  tree->Branch("lowlimit", &lowlimit, "lowlimit/D");
  for(unsigned int itype=0; itype<all_types; ++itype){
    graph->Clear();
    if( std::string(outputfile).find("quant0.027")!=std::string::npos ){
      if( itype == minus_2sigma ){
	fillTree(tree, graph, limit, lowlimit, itype, tanb_values, upper_exclusion, verbosity);
	break;
      }
    }
    else if( std::string(outputfile).find("quant0.160")!=std::string::npos ){
      if( itype == minus_1sigma ){
	fillTree(tree, graph, limit, lowlimit, itype, tanb_values, upper_exclusion, verbosity);
	break;
      }
    }
    else if( std::string(outputfile).find("quant0.500")!=std::string::npos ){
      if( itype == expected ){
	fillTree(tree, graph, limit, lowlimit, itype, tanb_values, upper_exclusion, verbosity);
	break;
      }
    }
    else if( std::string(outputfile).find("quant0.840")!=std::string::npos ){
      if( itype == plus_1sigma ){
	fillTree(tree, graph, limit, lowlimit, itype, tanb_values, upper_exclusion, verbosity);
	break;
      }
    }
    else if( std::string(outputfile).find("quant0.975")!=std::string::npos ){
      if( itype == plus_2sigma ){
	fillTree(tree, graph, limit, lowlimit, itype, tanb_values, upper_exclusion, verbosity);
	break;
      }
    }
    else{
      if( itype == observed ){
	fillTree(tree, graph, limit, lowlimit, itype, tanb_values, upper_exclusion, verbosity);
      }
    }
  }
  file->cd();
  tree->Write();
  file->Close();
  delete graph;

  return;
}

//.........这里部分代码省略.........
      for (int i= 0; i < statprecision_vs_ring[k]-> GetN(); i++){
	double spread, espread;
	double stat, estat;
	double residual, eresidual;
	double xdummy,ex;
	sigma_vs_ring[k]-> GetPoint(i, xdummy, spread );
	espread = sigma_vs_ring[k]-> GetErrorY(i);
	statprecision_vs_ring[k]-> GetPoint(i, xdummy, stat );
	estat = statprecision_vs_ring[k]-> GetErrorY(i);
	ex = statprecision_vs_ring[k]-> GetErrorX(i);
	if (spread > stat ){
	  residual  = sqrt( spread*spread - stat*stat );
	  eresidual = sqrt( pow(spread*espread,2) + pow(stat*estat,2))/residual;
	}
	else {
	  residual = 0;
	  eresidual = 0;
	}
	residual_vs_ring[k]->SetPoint(i,xdummy, residual);
	residual_vs_ring[k]->SetPointError(i,ex,eresidual);
      }
    }
 
 }
  
  ///-----------------------------------------------------------------
  ///--- Draw plots
  ///-----------------------------------------------------------------
  TCanvas *cEEP[12];
  TCanvas *cEEM[12];

  /// --- plot 0 : map of coefficients 
  cEEP[0] = new TCanvas("cmapEEP","cmapEEP");
  cEEP[0] -> cd();
  cEEP[0]->SetLeftMargin(0.1); 
  cEEP[0]->SetRightMargin(0.13); 
  cEEP[0]->SetGridx();
  cEEP[0]->SetGridy();
  //  hcmap[1]->GetXaxis()->SetNdivisions(1020);
  hcmap[1]->GetXaxis() -> SetLabelSize(0.03);
  hcmap[1]->Draw("COLZ");
  hcmap[1]->GetXaxis() ->SetTitle("ix");
  hcmap[1]->GetYaxis() ->SetTitle("iy");
  hcmap[1]->GetZaxis() ->SetRangeUser(0.8,1.2);

  cEEM[0] = new TCanvas("cmapEEM","cmapEEM");
  cEEM[0] -> cd();
  cEEM[0]->SetLeftMargin(0.1); 
  cEEM[0]->SetRightMargin(0.13); 
  cEEM[0]->SetGridx();
  cEEM[0]->SetGridy();
  //hcmap[0]->GetXaxis()->SetNdivisions(1020);
  hcmap[0]->GetXaxis() -> SetLabelSize(0.03);
  hcmap[0]->Draw("COLZ");
  hcmap[0]->GetXaxis() ->SetTitle("ix");
  hcmap[0]->GetYaxis() ->SetTitle("iy");
  hcmap[0]->GetZaxis() ->SetRangeUser(0.8,1.2);

  /// --- plot 1 : ring precision vs ieta
  cEEP[1] = new TCanvas("csigmaEEP","csigmaEEP");
  cEEP[1]->SetGridx();
  cEEP[1]->SetGridy();
  sigma_vs_ring[1]->GetHistogram()->GetYaxis()-> SetRangeUser(0.00,0.20);
  sigma_vs_ring[1]->GetHistogram()->GetXaxis()-> SetRangeUser(-85,85);
  sigma_vs_ring[1]->GetHistogram()->GetYaxis()-> SetTitle("#sigma_{c}");
  sigma_vs_ring[1]->GetHistogram()->GetXaxis()-> SetTitle("ring");

//.........这里部分代码省略.........
  Double_t pTPC = 0.; // Momentum at TPC inner wall
  UShort_t tpcSignalN = 0; // Number of clusters used for dEdx
  UChar_t  pidType = 0;
  Int_t    fMultiplicity = 0;
  //Double_t phiPrime = 0;

  // Only activate the branches of interest to save processing time
  tree->SetBranchStatus("*", 0); // Disable all branches
  tree->SetBranchStatus("pTPC", 1);
  tree->SetBranchStatus("dEdx", 1);
  tree->SetBranchStatus("dEdxExpected", 1);
  tree->SetBranchStatus("tanTheta", 1);
  tree->SetBranchStatus("tpcSignalN", 1);
  tree->SetBranchStatus("pidType", 1);
  //tree->SetBranchStatus("phiPrime", 1);
  if (isNonPP)
    tree->SetBranchStatus("fMultiplicity", 1);

  
  tree->SetBranchAddress("dEdx", &dEdx);
  tree->SetBranchAddress("dEdxExpected", &dEdxExpected);
  tree->SetBranchAddress("tanTheta", &tanTheta);
  tree->SetBranchAddress("tpcSignalN", &tpcSignalN);
  tree->SetBranchAddress("pTPC", &pTPC);
  tree->SetBranchAddress("pidType", &pidType);
  //tree->SetBranchAddress("phiPrime", &phiPrime);
  if (isNonPP)
    tree->SetBranchAddress("fMultiplicity", &fMultiplicity);

  
  // Output file
  TDatime daTime;
  TString savefileName = Form("%s%s_checkPullSigma_%04d_%02d_%02d__%02d_%02d.root", fileNameTree.ReplaceAll(".root", "").Data(),
                              recalculateExpecteddEdx ? "_recalcdEdx" : "",
                              daTime.GetYear(), daTime.GetMonth(), daTime.GetDay(), daTime.GetHour(), daTime.GetMinute());

  TFile* fSave = TFile::Open(Form("%s/%s", pathTree.Data(), savefileName.Data()), "recreate");
  if (!fSave) {
    std::cout << "Failed to open save file \"" << Form("%s/%s", pathTree.Data(), savefileName.Data()) << "\"!" << std::endl;
    return -1;
  }
  
  const Double_t pBoundLow = 0.1;
  const Double_t pBoundUp = 5;

  const Int_t nBins1 = TMath::Ceil(180 / downScaleFactor);
  const Int_t nBins2 = TMath::Ceil(100 / downScaleFactor);
  const Int_t nBins3 = TMath::Ceil(60 / downScaleFactor);
  
  const Int_t nPbinsForMap = nBins1 + nBins2 + nBins3;
  Double_t binsPforMap[nPbinsForMap + 1];
  
  Double_t binWidth1 = (1.0 - pBoundLow) / nBins1;
  Double_t binWidth2 = (2.0 - 1.0 ) / nBins2;
  Double_t binWidth3 = (pBoundUp - 2.0) / nBins3;
  
  for (Int_t i = 0; i < nBins1; i++)  {
    binsPforMap[i] = pBoundLow + i * binWidth1;
  }
  for (Int_t i = nBins1, j = 0; i < nBins1 + nBins2; i++, j++)  {
    binsPforMap[i] = 1.0 + j * binWidth2;
  }
  for (Int_t i = nBins1 + nBins2, j = 0; i < nBins1 + nBins2 + nBins3; i++, j++)  {
    binsPforMap[i] = 2.0 + j * binWidth3;
  }
  binsPforMap[nPbinsForMap] = pBoundUp;

double TrimEventContent(Int_t iRapBin = 1,
		      Int_t iPTBin = 1,
		      Double_t fracL = 0.5, Double_t nSigma = 2., 
		      Int_t nUpsState=0,//[0]... 1S, [1]... 2S, [2]... 3S
		      bool UpsMC=false,
		      bool f_BG_zero=false,
		      bool ProjectLSBdata=false,
		      bool ProjectRSBdata=false,
		      bool CombineSignalPeaks=false,
		      bool Y1Sto2S_SB=false,
		      bool LeftSided=false,
		      bool RightSided=false,
		      bool MassScan=false,
		      bool adjustOverlapBorders=true
		      ){

  printf("\n\n\nfracL = %1.3f, nSigma = %1.1f, iState = %d, rap %d, pT %d\n", fracL, nSigma, nUpsState, iRapBin, iPTBin);

  Char_t name[100], title[100];
  Char_t fileNameIn[100];
  sprintf(fileNameIn, "tmpFiles/data_Ups_rap%d_pT%d.root", iRapBin, iPTBin);
  //==============================
  //read inputs from input file:
  TFile *fIn = new TFile(fileNameIn);
  TLorentzVector *lepP;
  TLorentzVector *lepN;
  TTree *treeIn = (TTree *) gDirectory->Get("selectedData");
  if(gDirectory->Get("selectedData")==NULL){
    printf("\n\n\nskip processing this bin.\n\n\n");
    return -999.;
  }

  TH2D *hBG_cosThetaPhiLR[onia::kNbFrames][2];
  for(int iFrame = 0; iFrame < onia::kNbFrames; iFrame++){
    sprintf(name, "hBG_cosThetaPhi_%s_L", onia::frameLabel[iFrame]);
    hBG_cosThetaPhiLR[iFrame][L] = (TH2D *) gDirectory->Get(name);
    sprintf(name, "hBG_cosThetaPhi_%s_R", onia::frameLabel[iFrame]);
    hBG_cosThetaPhiLR[iFrame][R] = (TH2D *) gDirectory->Get(name);
  }
  //==============================

  //definition of output variables 
  Char_t fileNameOut[100];
  sprintf(fileNameOut, "AllStates_%1.2fSigma_FracLSB%dPercent/data_%dSUps_rap%d_pT%d.root", nSigma, int(fracL*100), nUpsState+1, iRapBin, iPTBin);
  TFile *fOut = new TFile(fileNameOut, "RECREATE");
  gStyle->SetPadRightMargin(0.2);
  TTree *treeOut = treeIn->CloneTree(0);
  // treeOut->SetName("data");
  TH2D *hBG_cosThetaPhi[onia::kNbFrames];
  // TH2D *hBG_cosThetaPhiSignal[onia::kNbFrames];
  for(int iFrame = 0; iFrame < onia::kNbFrames; iFrame++){
    // //book the histo for the signal
    // sprintf(name, "total_%s", onia::frameLabel[iFrame]);
    // sprintf(title, ";cos#theta_{%s};#phi_{%s} [deg]", onia::frameLabel[iFrame], onia::frameLabel[iFrame]);
    // hBG_cosThetaPhiSignal[iFrame] = new TH2D(name, title, onia::kNbBinsCosT, onia::cosTMin, onia::cosTMax, 
    // 					  onia::kNbBinsPhiPol, onia::phiPolMin, onia::phiPolMax);
    // hBG_cosThetaPhiSignal[iFrame]->Sumw2();
    //copy the L and R sideband histos into one output BG histogram
    hBG_cosThetaPhiLR[iFrame][L]->Scale(fracL/hBG_cosThetaPhiLR[iFrame][L]->Integral());
    hBG_cosThetaPhiLR[iFrame][R]->Scale((1.-fracL)/hBG_cosThetaPhiLR[iFrame][R]->Integral());
    sprintf(name, "background_costhphi%s", onia::frameLabel[iFrame]);
    hBG_cosThetaPhi[iFrame] = (TH2D *) hBG_cosThetaPhiLR[iFrame][L]->Clone(name);
    hBG_cosThetaPhi[iFrame]->Add(hBG_cosThetaPhiLR[iFrame][R]);
  }

  //==========================================================
  //reading fit parameters to establish signal mass window
  //as well as the L and R sideband window for the 3D BG histo
  //==========================================================
  fIn->cd();
  TTree *treeFitPar = (TTree *) gDirectory->Get("massFitParameters");
  TF1 *fUps[kNbSpecies], *fBG = 0;
  fUps[0] = 0, fUps[1] = 0, fUps[2] = 0;
  treeFitPar->SetBranchAddress("fUps1S", &fUps[0]);
  treeFitPar->SetBranchAddress("fUps2S", &fUps[1]);
  treeFitPar->SetBranchAddress("fUps3S", &fUps[2]);
  treeFitPar->SetBranchAddress("fBG", &fBG);
  treeFitPar->LoadTree(0);
  treeFitPar->GetEntry(0);


  Double_t mass[kNbSpecies], sigma[kNbSpecies];
  for(int iState = 0; iState < kNbSpecies; iState++){
    mass[iState] = fUps[iState]->GetParameter(1);
    sigma[iState] = fUps[iState]->GetParameter(2);
  }
  printf("1S: mass = %1.3f GeV, sigma = %1.3f GeV\n", mass[UPS1S], sigma[UPS1S]);
  printf("2S: mass = %1.3f GeV, sigma = %1.3f GeV\n", mass[UPS2S], sigma[UPS2S]);
  printf("3S: mass = %1.3f GeV, sigma = %1.3f GeV\n", mass[UPS3S], sigma[UPS3S]);
  Double_t poleMass = mass[nUpsState], massMin, massMax;
  massMin = poleMass - nSigma*sigma[nUpsState];
  massMax = poleMass + nSigma*sigma[nUpsState];

  if(LeftSided){
	  massMin = poleMass - nSigma*sigma[nUpsState];
	  massMax = poleMass;
  }

  if(RightSided){
	  massMin = poleMass;
//.........这里部分代码省略.........

void hadd(std::string& targetName, std::vector<std::string>& sources)
{
    using namespace std;
    TFile* target = 0;
    if(isForce) target = TFile::Open(targetName.c_str(), "RECREATE");
    else        target = TFile::Open(targetName.c_str(),   "CREATE");
    if(!target)
    {
        printf("target file already exists! (use -f to force recreation)\n");
        exit(0);
    }

    map<pair<string, string>, TObject*> outputMap;
    vector<pair<string, TObject*> > outputVec;

    for(vector<string>::const_iterator iF = sources.begin(); iF != sources.end(); ++iF)
    {
        if(verbosity >= 2) 
        {
            printf("Processing source file: %s\n", iF->c_str());
            fflush(stdout);
        }
        TFile * f = new TFile(iF->c_str());
        MergeRootfile(outputMap, outputVec, f, f);
        f->Close();
    }

    target->cd();
    map<string, TDirectory*> paths;
    for(std::vector<std::pair<std::string, TObject*> >::const_iterator iO = outputVec.begin(); iO != outputVec.end(); ++iO)
    {
        if(iO->second == 0)
        {
            if(paths.find(iO->first) == paths.end())
            {
                size_t pos = iO->first.rfind('/');
                if(pos == size_t(-1))
                {
                    target->cd();
                    pos = 0;
                }
                else
                {
                    target->cd(iO->first.substr(0, pos).c_str());
                    pos++;
                }
                paths[iO->first] = gDirectory->mkdir(iO->first.substr(pos).c_str());
            }
        }
        else
        {
            paths[iO->first]->cd();
            if(iO->second->IsA()->InheritsFrom(TH1::Class())) iO->second->Write();
            else
            {
                if(iO->second && ((TTree*)iO->second)->GetTreeIndex()) ((TTree*)iO->second)->GetTreeIndex()->Append(0, kFALSE); // Force the sorting
                TTree* tree = ((TTree*)iO->second)->CloneTree();
                ((TTree*)iO->second)->GetListOfClones()->Remove(tree);
                ((TTree*)iO->second)->ResetBranchAddresses();
                tree->ResetBranchAddresses();
                tree->Write();
            }
        }
    }
    
    if(verbosity >= 3)
    {
        printf("Results written to target file: %s\n", target->GetName());
        fflush(stdout);
    }

    target->Close();
    for(vector<pair<string, TFile*> >::const_iterator iF = tmpFiles.begin(); iF != tmpFiles.end(); ++iF)
    {
        if(iF->second)
        {
            iF->second->Close();
            system(("rm " + iF->first).c_str());
        }

    }
}

//.........这里部分代码省略.........
  
    //inputdata="/data/bmeson/data//nt_20141022_PPMuon_Run2013A_PromptReco_v1_resub20141126.root";
    //inputmc="/afs/cern.ch/work/w/wangj/public/nt_20140801_mixed_fromQMBFinder_Kp.root";
    inputdata="/data/bmeson/data/nt_20141022_PPMuon_Run2013A_PromptReco_v1_resub20141126_HLT_PAL1DoubleMu0_HighQ_v1.root";
    inputmc="/data/bmeson/MC/nt_2015251_PYTHIA6_BuJpsiK_pp_PAL1DoubleMu0_HighQ_v1.root";

    luminosity=5400*1e-3;
    outputname="../ResultsBplus_pp/SigmaBplus_extend.root";
    cut="abs(y)<2.4&&(HLT_PAL1DoubleMu0_HighQ_v1)&&abs(mumumass-3.096916)<0.15&&mass>5&&mass<6&& isbestchi2&&trk1Pt>0.9&&chi2cl>1.32e-02&&(d0/d0Err)>3.41&&cos(dtheta)>-3.46e01&&mu1pt>1.5&&mu2pt>1.5";
  }
  
  seldata_2y=Form("%s",cut.Data());
  selmc=Form("abs(y)<2.4&&gen==23333&&%s",cut.Data());
  selmcgen="abs(y)<2.4&&abs(pdgId)==521&&isSignal==1";

  if (doweight==0) weight="1";
  if (infname=="") infname=inputdata.Data();
  TFile *inf = new TFile(infname.Data());
  TTree *nt = (TTree*) inf->Get("ntKp");

  TFile *infMC = new TFile(inputmc.Data());
  TTree *ntGen = (TTree*)infMC->Get("ntGen");
  TTree *ntGen2 = (TTree*)inf->Get("ntGen");
  TTree *ntMC = (TTree*)infMC->Get("ntKp");
    
  ntGen->AddFriend(ntMC);
  ntGen2->AddFriend(ntMC);
    
  const int nBins = 5;
  double ptBins[nBins+1] = {10,15,20,25,30,60};
  //const int nBins = 1;
  //double ptBins[nBins+1] = {10,60};
  TH1D *hPt = new TH1D("hPt","",nBins,ptBins);
  TH1D *hPtRecoTruth = new TH1D("hPtRecoTruth","",nBins,ptBins);
  TH1D *hGenPtSelected = new TH1D("hGenPtSelected","",nBins,ptBins);
  TH1D *hPtMC = new TH1D("hPtMC","",nBins,ptBins);
  TH1D *hPtGen = new TH1D("hPtGen","",nBins,ptBins);
  TH1D *hPtGen2 = new TH1D("hPtGen2","",nBins,ptBins);

  TFile *outf = new TFile(outputname.Data(),"recreate");

  for (int i=0;i<nBins+1;i++)
    {
      if (i==nBins) {fit(nt,ntMC,10,60,ispPb,i);continue;}
      TF1 *f = fit(nt,ntMC,ptBins[i],ptBins[i+1],ispPb,i);
      double yield = f->Integral(5,6)/0.02;
      double yieldErr = f->Integral(5,6)/0.02*f->GetParError(0)/f->GetParameter(0);
      hPt->SetBinContent(i+1,yield/(ptBins[i+1]-ptBins[i]));
      hPt->SetBinError(i+1,yieldErr/(ptBins[i+1]-ptBins[i]));
    }  

  TCanvas *c=  new TCanvas("cResult","",600,600);
  hPt->SetXTitle("B^{+} p_{T} (GeV/c)");
  hPt->SetYTitle("Uncorrected B^{+} dN/dp_{T}");
  hPt->Sumw2();
  hPt->Draw();
  
  ntMC->Project("hPtMC","pt",TCut(weight)*(TCut(selmc.Data())&&"gen==23333"));
  nt->Project("hPtRecoTruth","pt",TCut(seldata_2y.Data())&&"gen==23333");
  ntGen->Project("hPtGen","pt",TCut(weight)*(TCut(selmcgen.Data())));
  ntGen2->Project("hPtGen2","pt",(TCut(selmcgen.Data())));
  divideBinWidth(hPtRecoTruth);
  
  hPtRecoTruth->Draw("same hist");
  divideBinWidth(hPtMC);
  divideBinWidth(hPtGen);
  divideBinWidth(hPtGen2);
  
  hPtMC->Sumw2();
  TH1D *hEff = (TH1D*)hPtMC->Clone("hEff");
  hPtMC->Sumw2();
  hEff->Divide(hPtGen);
  
  TH1D *hPtCor = (TH1D*)hPt->Clone("hPtCor");
  hPtCor->Divide(hEff);
  TCanvas *cCor=  new TCanvas("cCorResult","",600,600);
  hPtCor->SetYTitle("Corrected B^{+} dN/dp_{T}");
  hPtCor->Draw();
  hPtGen->Draw("same hist");
  hPtGen2->Draw("same hist");

  TH1D *hPtSigma= (TH1D*)hPtCor->Clone("hPtSigma");
  double BRchain=6.09604e-5;

  hPtSigma->Scale(1./(2*luminosity*BRchain));
  hPtSigma->SetYTitle("d#sigma (B^{+})/dp_{T}");

  TCanvas *cSigma=  new TCanvas("cSigma","",600,600);

  hPtSigma->Draw();
  
  outf->cd();
  hPt->Write();
  hEff->Write();
  hPtGen->Write();
  hPtCor->Write();
  hPtSigma->Write();
  outf->Close();
  delete outf;
}

//.........这里部分代码省略.........
  TH2D *h2D_Acc_pt_y_Pbp = new TH2D("h2D_Acc_pt_y_Pbp","",nRapBins,rapBinsArr,nPtBins,ptBinsArr);
  h2D_Den_pt_y_Pbp->Sumw2();
  h2D_Num_pt_y_Pbp->Sumw2();
  h2D_Acc_pt_y_Pbp->Sumw2();
  //fine bins for plots in AN (for 8rap9pt only)
  TH2D *h2D_Den_pt_y_Pbp_fine = new TH2D("h2D_Den_pt_y_Pbp_fine","",100,-2.4,2.4,100,0,30);
  TH2D *h2D_Num_pt_y_Pbp_fine = new TH2D("h2D_Num_pt_y_Pbp_fine","",100,-2.4,2.4,100,0,30);
  TH2D *h2D_Acc_pt_y_Pbp_fine = new TH2D("h2D_Acc_pt_y_Pbp_fine","",100,-2.4,2.4,100,0,30);
  h2D_Den_pt_y_Pbp_fine->Sumw2();
  h2D_Num_pt_y_Pbp_fine->Sumw2();
  h2D_Acc_pt_y_Pbp_fine->Sumw2();

  ////////////////////////////////////////////////////////////////////////////////////////////////////  
  //// event loop Start!
  for(int i=0; i<ana1->GetEntries(); i++){
    if(i%100000==0) cout << ">>>>> EVENT " << i << " / " << ana1->GetEntries() <<  endl;
    ana1->GetEntry(i);
    //// for denominator
    if (dimuCut(dmom0_Id,dkid0_ch,dkid1_ch) && minpt<=dimu_pt && dimu_pt<maxpt && minylab<=dimu_y && dimu_y<maxylab) { 
      h2D_Den_pt_y_Pbp->Fill(dimu_y,dimu_pt);
      if (mrapnpt==89) h2D_Den_pt_y_Pbp_fine->Fill(dimu_y,dimu_pt);
      //// for numerator 
      if (massCut(dimu_mass) && kineCut(dkid0_pt,dkid0_eta,dkid0_p) && kineCut(dkid1_pt,dkid1_eta,dkid1_p)) { 
        h2D_Num_pt_y_Pbp->Fill(dimu_y,dimu_pt);
        if (mrapnpt==89) h2D_Num_pt_y_Pbp_fine->Fill(dimu_y,dimu_pt);
      }
    }
  }//// end of event loop

  //// (Num/Den) to get acceptance (B : binomial error)
  h2D_Acc_pt_y_Pbp->Divide(h2D_Num_pt_y_Pbp,h2D_Den_pt_y_Pbp,1,1,"B");
  if (mrapnpt==89) h2D_Acc_pt_y_Pbp_fine->Divide(h2D_Num_pt_y_Pbp_fine,h2D_Den_pt_y_Pbp_fine,1,1,"B");

  ////////////////////////////////////////////////////////////////////////////////////////////////////  
  //// rapidity reverse, for 2nd run hist.
  TH2D *h2D_Den_pt_y_pPb = new TH2D("h2D_Den_pt_y_pPb","",nRapBins2nd,rapBinsArr2nd,nPtBins,ptBinsArr);
  TH2D *h2D_Num_pt_y_pPb = new TH2D("h2D_Num_pt_y_pPb","",nRapBins2nd,rapBinsArr2nd,nPtBins,ptBinsArr);
  TH2D *h2D_Acc_pt_y_pPb = new TH2D("h2D_Acc_pt_y_pPb","",nRapBins2nd,rapBinsArr2nd,nPtBins,ptBinsArr);
  h2D_Den_pt_y_pPb->Sumw2();
  h2D_Num_pt_y_pPb->Sumw2();
  h2D_Acc_pt_y_pPb->Sumw2();

  int tmpbinDenPbp, tmpbinNumPbp, tmpbinAccPbp;
  int tmpbinDenpPb, tmpbinNumpPb, tmpbinAccpPb;
  float tmpDenVal, tmpDenErr, tmpNumVal, tmpNumErr, tmpAccVal, tmpAccErr;
  
  for (Int_t iy=0; iy<nRapBins; iy++){
    for (Int_t ipt=0; ipt<nPtBins; ipt++) {
      tmpbinDenPbp = h2D_Den_pt_y_Pbp->FindBin((rapBinsArr[iy]+rapBinsArr[iy+1])/2, (ptBinsArr[ipt]+ptBinsArr[ipt+1])/2);
      tmpDenVal = h2D_Den_pt_y_Pbp->GetBinContent(tmpbinDenPbp);
      tmpDenErr = h2D_Den_pt_y_Pbp->GetBinError(tmpbinDenPbp);
      tmpbinDenpPb = h2D_Den_pt_y_pPb->FindBin((rapBinsArr2nd[nRapBins-iy-1]+rapBinsArr2nd[nRapBins-iy])/2,(ptBinsArr[ipt]+ptBinsArr[ipt+1])/2);
      h2D_Den_pt_y_pPb->SetBinContent(tmpbinDenpPb,tmpDenVal);
      h2D_Den_pt_y_pPb->SetBinError(tmpbinDenpPb,tmpDenErr);
      
      tmpbinNumPbp = h2D_Num_pt_y_Pbp->FindBin((rapBinsArr[iy]+rapBinsArr[iy+1])/2, (ptBinsArr[ipt]+ptBinsArr[ipt+1])/2);
      tmpNumVal = h2D_Num_pt_y_Pbp->GetBinContent(tmpbinNumPbp);
      tmpNumErr = h2D_Num_pt_y_Pbp->GetBinError(tmpbinNumPbp);
      tmpbinNumpPb = h2D_Num_pt_y_pPb->FindBin((rapBinsArr2nd[nRapBins-iy-1]+rapBinsArr2nd[nRapBins-iy])/2,(ptBinsArr[ipt]+ptBinsArr[ipt+1])/2);
      h2D_Num_pt_y_pPb->SetBinContent(tmpbinNumpPb,tmpNumVal);
      h2D_Num_pt_y_pPb->SetBinError(tmpbinNumpPb,tmpNumErr);
      
      tmpbinAccPbp = h2D_Acc_pt_y_Pbp->FindBin((rapBinsArr[iy]+rapBinsArr[iy+1])/2, (ptBinsArr[ipt]+ptBinsArr[ipt+1])/2);
      tmpAccVal = h2D_Acc_pt_y_Pbp->GetBinContent(tmpbinAccPbp);
      tmpAccErr = h2D_Acc_pt_y_Pbp->GetBinError(tmpbinAccPbp);
      tmpbinAccpPb = h2D_Acc_pt_y_pPb->FindBin((rapBinsArr2nd[nRapBins-iy-1]+rapBinsArr2nd[nRapBins-iy])/2,(ptBinsArr[ipt]+ptBinsArr[ipt+1])/2);
      h2D_Acc_pt_y_pPb->SetBinContent(tmpbinAccpPb,tmpAccVal);
      h2D_Acc_pt_y_pPb->SetBinError(tmpbinAccpPb,tmpAccErr);
      /*
      cout << "*** "<<iy<<"th rap, "<<ipt<<"th pt ***" <<endl;
      cout << "tmpbinAccPbp = " <<tmpbinAccPbp << ", tmpbinAccpPb = " << tmpbinAccpPb << endl;
      cout << "tmpAccVal = " <<tmpAccVal << ", tmpAccErr = " << tmpAccErr << endl;
      cout << "rapBinsArr :" <<(rapBinsArr[iy]+rapBinsArr[iy+1])/2 << endl;
      cout << "      from " <<rapBinsArr[iy]<<" to "<<rapBinsArr[iy+1] << endl;
      cout << "rapBinsArr2nd : " << (rapBinsArr2nd[nRapBins-iy-1]+rapBinsArr2nd[nRapBins-iy])/2 << endl;
      cout << "      from " << rapBinsArr2nd[nRapBins-iy-1] << " to " << rapBinsArr2nd[nRapBins-iy] << endl;  
      cout << "" << endl;
      */
    }
  }

  ////////////////////////////////////////////////////////////////////////////////////////////////////  
  //// Save the data as an root file
  TFile *outFile = new TFile(Form("AccAna_%s.root",strName.Data()),"RECREATE");
  std::cout << "strName: " << strName << std::endl;
  outFile->cd();
  h2D_Den_pt_y_Pbp->Write();
  h2D_Num_pt_y_Pbp->Write();
  h2D_Acc_pt_y_Pbp->Write();
  h2D_Den_pt_y_pPb->Write();
  h2D_Num_pt_y_pPb->Write();
  h2D_Acc_pt_y_pPb->Write();
  if (mrapnpt==89) {
    h2D_Den_pt_y_Pbp_fine->Write();
    h2D_Num_pt_y_Pbp_fine->Write();
    h2D_Acc_pt_y_Pbp_fine->Write();
  }
  outFile->Close();

} // end of main func

void GetHiggsBR (int maxEvents) {

  TChain* chain = new TChain(inputTree);

  EvtInfo.Register(chain);
  GenInfo.Register(chain);

  for (int ii = 0; ii < 10; ++ii) {
    chain -> Add(fnames[ii]);
    TFile *fin = TFile::Open(fnames[ii],"READ");
    fin->Close();
  }

  TFile* fout = new TFile("HiggsBR.root", "RECREATE") ; 
  fout->cd() ; 
  TH1D* h_NHiggs = new TH1D("h_NHiggs", "Higgs boson decay channels;Decay channel;Events;", 12, 0., 12.) ; 
  h_NHiggs -> GetXaxis() -> SetBinLabel(1 , "H->bb") ; 
  h_NHiggs -> GetXaxis() -> SetBinLabel(2 , "H->#tau#tau") ; 
  h_NHiggs -> GetXaxis() -> SetBinLabel(3 , "H->#mu#mu") ; 
  h_NHiggs -> GetXaxis() -> SetBinLabel(4 , "H->cc") ; 
  h_NHiggs -> GetXaxis() -> SetBinLabel(5 , "H->ss") ; 
  h_NHiggs -> GetXaxis() -> SetBinLabel(6 , "H->tt") ; 
  h_NHiggs -> GetXaxis() -> SetBinLabel(7 , "H->gg") ; 
  h_NHiggs -> GetXaxis() -> SetBinLabel(8 , "H->#gamma#gamma") ; 
  h_NHiggs -> GetXaxis() -> SetBinLabel(9 , "H->Z#gamma") ; 
  h_NHiggs -> GetXaxis() -> SetBinLabel(10, "H->WW") ; 
  h_NHiggs -> GetXaxis() -> SetBinLabel(11, "H->ZZ") ; 
  h_NHiggs -> GetXaxis() -> SetBinLabel(12, "H->everything") ; 
  TH1D* h_BRHiggs = new TH1D("h_BRHiggs", "Higgs boson BRs;Decay channels;BR", 12, 0., 12.) ; 
  h_BRHiggs -> GetXaxis() -> SetBinLabel(1 , "H->bb") ; 
  h_BRHiggs -> GetXaxis() -> SetBinLabel(2 , "H->#tau#tau") ; 
  h_BRHiggs -> GetXaxis() -> SetBinLabel(3 , "H->#mu#mu") ; 
  h_BRHiggs -> GetXaxis() -> SetBinLabel(4 , "H->cc") ; 
  h_BRHiggs -> GetXaxis() -> SetBinLabel(5 , "H->ss") ; 
  h_BRHiggs -> GetXaxis() -> SetBinLabel(6 , "H->tt") ; 
  h_BRHiggs -> GetXaxis() -> SetBinLabel(7 , "H->gg") ; 
  h_BRHiggs -> GetXaxis() -> SetBinLabel(8 , "H->#gamma#gamma") ; 
  h_BRHiggs -> GetXaxis() -> SetBinLabel(9 , "H->Z#gamma") ; 
  h_BRHiggs -> GetXaxis() -> SetBinLabel(10, "H->WW") ; 
  h_BRHiggs -> GetXaxis() -> SetBinLabel(11, "H->ZZ") ; 
  h_BRHiggs -> GetXaxis() -> SetBinLabel(12, "H->everything") ; 
  TH1D* h_BRSFH125 = new TH1D("h_BRSFH125", "Higgs boson SF(BR);Decay channels;SF(BR)", 12, 0., 12.) ; 
  h_BRSFH125 -> GetXaxis() -> SetBinLabel(1 , "H->bb") ; 
  h_BRSFH125 -> GetXaxis() -> SetBinLabel(2 , "H->#tau#tau") ; 
  h_BRSFH125 -> GetXaxis() -> SetBinLabel(3 , "H->#mu#mu") ; 
  h_BRSFH125 -> GetXaxis() -> SetBinLabel(4 , "H->cc") ; 
  h_BRSFH125 -> GetXaxis() -> SetBinLabel(5 , "H->ss") ; 
  h_BRSFH125 -> GetXaxis() -> SetBinLabel(6 , "H->tt") ; 
  h_BRSFH125 -> GetXaxis() -> SetBinLabel(7 , "H->gg") ; 
  h_BRSFH125 -> GetXaxis() -> SetBinLabel(8 , "H->#gamma#gamma") ; 
  h_BRSFH125 -> GetXaxis() -> SetBinLabel(9 , "H->Z#gamma") ; 
  h_BRSFH125 -> GetXaxis() -> SetBinLabel(10, "H->WW") ; 
  h_BRSFH125 -> GetXaxis() -> SetBinLabel(11, "H->ZZ") ; 
  h_BRSFH125 -> GetXaxis() -> SetBinLabel(12, "H->everything") ; 

  if(maxEvents<0 || maxEvents>chain->GetEntries()) maxEvents = chain->GetEntries();

  for(int entry = 0; entry < maxEvents; entry++) { 
    if( (entry%1000) == 0 ) std::cout << " entry = " << entry << " of " << maxEvents << std::endl ;
    chain -> GetEntry(entry);

    for (int igen=0; igen < GenInfo.Size; ++igen) {

      if ( GenInfo.Status[igen] == 3 && 
          TMath::Abs(GenInfo.PdgID[igen])==25 && 
          GenInfo.nDa[igen] >= 2 ) { 

        int higgsDau = TMath::Abs(GenInfo.Da0PdgID[igen]) ; 
        if( TMath::Abs(GenInfo.Da0PdgID[igen]) == TMath::Abs(GenInfo.Da1PdgID[igen]) ) { 

          h_BRHiggs -> Fill ("H->everything" ,1.) ; h_NHiggs -> Fill ("H->everything" ,1.) ; 
          switch (higgsDau) {
            case 5 :
              h_BRHiggs -> Fill ("H->bb" ,1.) ; h_NHiggs -> Fill ("H->bb" ,1.) ; 
              break ; 
            case 15 : 
              h_BRHiggs -> Fill ("H->#tau#tau" ,1.) ; h_NHiggs -> Fill ("H->#tau#tau" ,1.) ; 
              break ; 
            case 13 :
              h_BRHiggs -> Fill ("H->#mu#mu" ,1.) ; h_NHiggs -> Fill ("H->#mu#mu" ,1.) ; 
              break ;
            case 4 :
              h_BRHiggs -> Fill ("H->cc" ,1.) ; h_NHiggs -> Fill ("H->cc" ,1.) ; 
              break ;
            case 3 :
              h_BRHiggs -> Fill ("H->ss" ,1.) ; h_NHiggs -> Fill ("H->ss" ,1.) ; 
              break ;
            case 6 :
              h_BRHiggs -> Fill ("H->tt" ,1.) ; h_NHiggs -> Fill ("H->tt" ,1.) ; 
              break ;
            case 21 :
              h_BRHiggs -> Fill ("H->gg" ,1.) ; h_NHiggs -> Fill ("H->gg" ,1.) ; 
              break ;
            case 22 : 
              h_BRHiggs -> Fill ("H->#gamma#gamma" ,1.) ; h_NHiggs -> Fill ("H->#gamma#gamma" ,1.) ; 
              break ;
            case 24 :
              h_BRHiggs -> Fill ("H->WW" ,1.) ; h_NHiggs -> Fill ("H->WW" ,1.) ; 
              break ;
            case 23 :
//.........这里部分代码省略.........

//.........这里部分代码省略.........
        treeHiForestInfo->GetEntry(j_entry);
        if(isMC) treeGen->GetEntry(j_entry);

        bool eventAdded = em->addEvent(run,lumis,event,j_entry);
        //std::cout << run << " " << lumis << " " << event << " " << j_entry << std::endl;
        if(!eventAdded) // this event is duplicate, skip this one.
        {
            duplicateEntries++;
            continue;
        }

        if(isMC) {
            if((pthat>=pthatCut[0]) && (pthat<pthatCut[1])) pthatWeight = tmpWeight[0];
            else if((pthat>=pthatCut[1]) && (pthat<pthatCut[2])) pthatWeight = tmpWeight[1];
            else if((pthat>=pthatCut[2]) && (pthat<pthatCut[3])) pthatWeight = tmpWeight[2];
            else if((pthat>=pthatCut[3]) && (pthat<pthatCut[4])) pthatWeight = tmpWeight[3];
            else if((pthat>=pthatCut[4]) && (pthat<pthatCut[5])) pthatWeight = tmpWeight[4];
            else continue;
        }
        // event selection
        if (!(TMath::Abs(vz) < cut_vz))  continue;
        if (isHI) {
            if ((pcollisionEventSelection < cut_pcollisionEventSelection))  continue;
        }
        else {
            if (pPAprimaryVertexFilter < cut_pPAprimaryVertexFilter || pBeamScrapingFilter < cut_pBeamScrapingFilter)  continue;
        }
        entriesPassedEventSelection++;

        // photon block
        // find leading photon
        int phoIdx = -1;     // index of the leading photon
        double maxPhoEt = -1;
        for(int i=0; i<ggHi.nPho; ++i)
        {
            bool failedEtCut  = (ggHi.phoEt->at(i) < cutPhoEt) ;
            bool failedEtaCut = (TMath::Abs(ggHi.phoEta->at(i)) > cutPhoEta) ;
            bool failedSpikeRejection;
            bool failedHotSpotRejection;
            //if (isHI) {
            failedSpikeRejection =( (ggHi.phoEta->at(i)<1.44) && 
                    (ggHi.phoSigmaIEtaIEta->at(i) < 0.002 ||
                     ggHi.pho_swissCrx->at(i)     > 0.9   ||
                     TMath::Abs(ggHi.pho_seedTime->at(i)) > 3) );
            // }
            // else {
            //     failedSpikeRejection = (ggHi.phoSigmaIEtaIEta->at(i) < 0.002);
            // }

            failedHotSpotRejection = (
                (ggHi.phoE3x3->at(i)/ggHi.phoE5x5->at(i) > 2./3.-0.03 && ggHi.phoE3x3->at(i)/ggHi.phoE5x5->at(i) < 2./3.+0.03) &&
                (ggHi.phoE1x5->at(i)/ggHi.phoE5x5->at(i) > 1./3.-0.03 && ggHi.phoE1x5->at(i)/ggHi.phoE5x5->at(i) < 1./3.+0.03) &&
                (ggHi.phoE2x5->at(i)/ggHi.phoE5x5->at(i) > 2./3.-0.03 && ggHi.phoE2x5->at(i)/ggHi.phoE5x5->at(i) < 2./3.+0.03) );
            
            bool failedHoverE = (ggHi.phoHoverE->at(i) > 0.2);      // <0.1 cut is applied after corrections
            //               bool failedEnergyRatio = ((float)ggHi.phoSCRawE->at(i)/ggHi.phoE->at(i) < 0.5);

            if (failedEtCut)          continue;
            if (failedEtaCut)         continue;
            if (failedSpikeRejection) continue;
            if (failedHotSpotRejection) {entriesSpikeRejected++; continue;}
            if (failedHoverE)         continue;
            //               if (failedEnergyRatio)    continue;    // actually applied after corrections

            if (ggHi.phoEt->at(i) > maxPhoEt)
            {
                maxPhoEt = ggHi.phoEt->at(i);
                phoIdx = i;
            }
        }
        if (phoIdx == -1) continue;
        entriesAnalyzed++;

        outputTreeHLT->Fill();
        outputTreeggHiNtuplizer->Fill();
        outputTreeHiEvt->Fill();
        outputTreeSkim->Fill();
        outputTreeHiForestInfo->Fill();
        if(isMC) outputTreeGen->Fill();
    }// event loop closed here

    std::cout<<  "Loop ENDED : ggHiNtuplizer/EventTree" <<std::endl;
    std::cout << "entries            = " << entries << std::endl;
    std::cout << "duplicateEntries   = " << duplicateEntries << std::endl;
    std::cout << "entriesPassedEventSelection   = " << entriesPassedEventSelection << std::endl;
    std::cout << "entriesAnalyzed               = " << entriesAnalyzed << std::endl;
    std::cout << "entriesSpikeRejected          = " << entriesSpikeRejected << std::endl;
    std::cout << "outputTreeHLT->GetEntries()   = " << outputTreeHLT->GetEntries() << std::endl;
    std::cout << "outputTreeggHiNtuplizer->GetEntries()   = " << outputTreeggHiNtuplizer->GetEntries() << std::endl;
    std::cout << "outputTreeSkim->GetEntries()  = " << outputTreeSkim->GetEntries() << std::endl;
    std::cout << "outputTreeHiEvt->GetEntries() = " << outputTreeHiEvt->GetEntries() << std::endl;
    std::cout << "outputTreeHiForestInfo->GetEntries() = " << outputTreeHiForestInfo->GetEntries() << std::endl;
    if(isMC) std::cout << "outputTreeGen->GetEntries() = " << outputTreeGen->GetEntries() << std::endl;


    output->cd();
    configTree->Write("",TObject::kOverwrite);
    output->Write("",TObject::kOverwrite);
    output->Close();
}

//.........这里部分代码省略.........
  } //-- End cent loop

  //-- Make sweet, sweet TGraphErrors
  double cErr[NCENT];
  for(int i = 0; i < NCENT; i++) cErr[i] = 0;
  grVn2DoSys_RatioToNominal       = new TGraphErrors(NCENT, centBinCenter, vn2DoSys_RatioToNominal, cErr, vn2DoSys_RatioToNominal_staterr);
  formatGraph(grVn2DoSys_RatioToNominal, "Centrality %", ratioMin, ratioMax, Form("v_{%i}{2} [RespErr] / [Nominal]", norder_), 1, 24, "grVn2DoSys_RatioToNominal");
  grVn4DoSys_RatioToNominal       = new TGraphErrors(NCENT, centBinCenter, vn4DoSys_RatioToNominal, cErr, vn4DoSys_RatioToNominal_staterr);
  formatGraph(grVn4DoSys_RatioToNominal, "Centrality %", ratioMin, ratioMax, Form("v_{%i}{4} [RespErr] / [Nominal]", norder_), kSpring+4, 25, "grVn4DoSys_RatioToNominal");
  grVn6DoSys_RatioToNominal       = new TGraphErrors(NCENT, centBinCenter, vn6DoSys_RatioToNominal, cErr, vn6DoSys_RatioToNominal_staterr);
  formatGraph(grVn6DoSys_RatioToNominal, "Centrality %", ratioMin, ratioMax, Form("v_{%i}{6} [RespErr] / [Nominal]", norder_), 6, 28, "grVn6DoSys_RatioToNominal");
  grVn8DoSys_RatioToNominal       = new TGraphErrors(NCENT, centBinCenter, vn8DoSys_RatioToNominal, cErr, vn8DoSys_RatioToNominal_staterr);
  formatGraph(grVn8DoSys_RatioToNominal, "Centrality %", ratioMin, ratioMax, Form("v_{%i}{8} [RespErr] / [Nominal]", norder_), kOrange+7, 27, "grVn8DoSys_RatioToNominal");
  grGamma1ExpDoSys_RatioToNominal = new TGraphErrors(NCENT, centBinCenter, gamma1expDoSys_RatioToNominal, cErr, gamma1expDoSys_RatioToNominal_staterr);
  formatGraph(grGamma1ExpDoSys_RatioToNominal, "Centrality %", ratioG1Min, ratioG1Max, "|1-#gamma_{1}^{exp} [RespErr] / [Nominal]|", 2, 20, "grGamma1ExpDoSys_RatioToNominal");
  grVn6Vn4DoSys_RatioToNominal    = new TGraphErrors(NCENT, centBinCenter, vn6vn4DoSys_RatioToNominal, cErr, vn6vn4DoSys_RatioToNominal_staterr);
  formatGraph(grVn6Vn4DoSys_RatioToNominal, "Centrality %", ratioCumuRatioMin, ratioCumuRatioMax, Form("v_{%i}{6}/v_{%i}{4} [RespErr] / [Nominal]", norder_, norder_), 4, 21, "grVn6Vn4DoSys_RatioToNominal");
  grVn8Vn4DoSys_RatioToNominal    = new TGraphErrors(NCENT, centBinCenter, vn8vn4DoSys_RatioToNominal, cErr, vn8vn4DoSys_RatioToNominal_staterr);
  formatGraph(grVn8Vn4DoSys_RatioToNominal, "Centrality %", ratioCumuRatioMin, ratioCumuRatioMax, Form("v_{%i}{8}/v_{%i}{4} [RespErr] / [Nominal]", norder_, norder_), kGreen+2, 34, "grVn8Vn4DoSys_RatioToNominal");
  grVn8Vn6DoSys_RatioToNominal    = new TGraphErrors(NCENT, centBinCenter, vn8vn6DoSys_RatioToNominal, cErr, vn8vn6DoSys_RatioToNominal_staterr);
  formatGraph(grVn8Vn6DoSys_RatioToNominal, "Centrality %", ratioCumuRatioMin, ratioCumuRatioMax, Form("v_{%i}{8}/v_{%i}{6} [RespErr] / [Nominal]", norder_, norder_), kViolet-1, 33, "grVn8Vn6DoSys_RatioToNominal");
  grVn46_Vn68DoSys_RatioToNominal    = new TGraphErrors(NCENT, centBinCenter, vn46_vn68DoSys_RatioToNominal, cErr, vn46_vn68DoSys_RatioToNominal_staterr);
  formatGraph(grVn46_Vn68DoSys_RatioToNominal, "Centrality %", 0.75, 1.25, Form("(v_{%i}{4} - v_{%i}{6})/(v_{%i}{6} - v_{%i}{8}) [RespErr] / [Nominal]", norder_, norder_, norder_, norder_), kGray+2, 22, "grVn46_Vn68DoSys_RatioToNominal");


  TLine * line = new TLine(grVn2DoSys_RatioToNominal->GetXaxis()->GetXmin(), 1.0, grVn2DoSys_RatioToNominal->GetXaxis()->GetXmax(), 1.0);
  line->SetLineColor(1);
  line->SetLineStyle(2);
  line->SetLineWidth(2);


  TCanvas * cCumuSys = new TCanvas("cCumuSys", "cCumuSys", 1000, 1000);
  cCumuSys->Divide(2,2);
  cCumuSys->cd(1);
  grVn2DoSys_RatioToNominal->Draw("ap");
  line->Draw("same");
  cCumuSys->cd(2);
  grVn4DoSys_RatioToNominal->Draw("ap");
  line->Draw("same");
  cCumuSys->cd(3);
  grVn6DoSys_RatioToNominal->Draw("ap");
  line->Draw("same");
  cCumuSys->cd(4);
  grVn8DoSys_RatioToNominal->Draw("ap");
  line->Draw("same");
  cCumuSys->SaveAs("../../plots/systematicStudies/cSysResp_CumuCent.pdf");

  TCanvas * cGamma1ExpSys = new TCanvas("cGamma1ExpSys", "cGamma1ExpSys", 500, 500);
  cGamma1ExpSys->cd();
  grGamma1ExpDoSys_RatioToNominal->Draw("ap");
  //line->Draw("same");
  cGamma1ExpSys->SaveAs("../../plots/systematicStudies/cSysResp_Gamma1ExpCent.pdf");

  TCanvas * cCumuRatioSys = new TCanvas("cCumuRatioSys", "cCumuRatioSys", 1500, 500);
  cCumuRatioSys->Divide(3,1);
  cCumuRatioSys->cd(1);
  cCumuRatioSys->cd(1)->SetLeftMargin(0.2);
  grVn6Vn4DoSys_RatioToNominal->GetYaxis()->SetTitleOffset(1.6);
  grVn6Vn4DoSys_RatioToNominal->Draw("ap");
  line->Draw("same");
  cCumuRatioSys->cd(2);
  grVn8Vn4DoSys_RatioToNominal->Draw("ap");
  line->Draw("same");
  cCumuRatioSys->cd(3);
  grVn8Vn6DoSys_RatioToNominal->Draw("ap");
  line->Draw("same");

void Draw2DPlot(TString name, TString histotitle, TString titleX, TString titleY, bool log){



 TH2F* hist_b;
 TH2F* hist_c;
 TH2F* hist_gsplit;
 TH2F* hist_l;
 TH2F* hist_data;
 
 
 TFile *myFile     = new TFile(filename);
 
 myFile->cd();
 hist_b         = (TH2F*)gROOT->FindObject(name+"_b");
 hist_c         = (TH2F*)gROOT->FindObject(name+"_c");
 hist_gsplit    = (TH2F*)gROOT->FindObject(name+"_bfromg");
 hist_l         = (TH2F*)gROOT->FindObject(name+"_l");
 hist_data      = (TH2F*)gROOT->FindObject(name+"_data");
 

 TH2F* histo_tot = (TH2F*) hist_b->Clone();
 histo_tot ->Add(hist_c);
 histo_tot ->Add(hist_gsplit);
 histo_tot ->Add(hist_l); 
 
 float scale_f = (hist_data->Integral())/(hist_b->Integral() + hist_c ->Integral()+ hist_gsplit->Integral() + hist_l->Integral());

 hist_b       ->Scale(scale_f);
 hist_c       ->Scale(scale_f);
 hist_gsplit  ->Scale(scale_f);
 hist_l       ->Scale(scale_f);
 histo_tot    ->Scale(scale_f);

  
 TProfile * pro_mc = histo_tot->ProfileX(name+"_tot");
 TProfile * pro_mc_b = hist_b->ProfileX();
 TProfile * pro_mc_c = hist_c->ProfileX();
 TProfile * pro_mc_udsg = hist_l->ProfileX();
 TProfile * pro_mc_gspl = hist_gsplit->ProfileX();
 TProfile * pro_data = hist_data->ProfileX();
 
   // SET COLORS
  pro_mc->SetLineColor(1);
  pro_mc_b->SetLineColor(2);
  pro_mc_c->SetLineColor(8);
  pro_mc_udsg->SetLineColor(4);
  pro_mc_gspl->SetLineColor(7);

  pro_data->SetMarkerStyle(20);
  pro_data->SetMarkerSize(0.75);
  
  pro_mc_gspl->GetXaxis()->SetTitle(titleX);  
  pro_mc_gspl->GetYaxis()->SetTitle(titleY);  
  
  
  // SET COSMETICS
  pro_data->SetMarkerStyle(20);
  pro_data->SetMarkerSize(0.75);
  //pro_mc_gspl->GetXaxis()->SetTitle();  
  //pro_mc_gspl->GetYaxis()->SetTitle();  


  gStyle->SetOptTitle(0);


  Double_t titleoffsetx=0.8;
  Double_t titleoffsety=0.8;
  Double_t titlesizex=0.05;
  Double_t titlesizey=0.05;
  Double_t labelsizex=0.035;
  Double_t labelsizey=0.035;

  pro_data->GetYaxis()->SetLabelSize(labelsizey);
  pro_data->GetYaxis()->SetTitleSize(titlesizey);
  pro_data->GetYaxis()->SetTitleOffset(titleoffsety);
  pro_mc->GetYaxis()->SetLabelSize(labelsizey);
  pro_mc->GetYaxis()->SetTitleSize(titlesizey);
  pro_mc->GetYaxis()->SetTitleOffset(titleoffsety);
  pro_mc_b->GetYaxis()->SetLabelSize(labelsizey);
  pro_mc_b->GetYaxis()->SetTitleSize(titlesizey);
  pro_mc_b->GetYaxis()->SetTitleOffset(titleoffsety);
  pro_mc_c->GetYaxis()->SetLabelSize(labelsizey);
  pro_mc_c->GetYaxis()->SetTitleSize(titlesizey);
  pro_mc_c->GetYaxis()->SetTitleOffset(titleoffsety);

  pro_mc_gspl->GetYaxis()->SetLabelSize(labelsizey);
  pro_mc_gspl->GetYaxis()->SetTitleSize(titlesizey);
  pro_mc_gspl->GetYaxis()->SetTitleOffset(titleoffsety);

  pro_mc_udsg->GetYaxis()->SetLabelSize(labelsizey);
  pro_mc_udsg->GetYaxis()->SetTitleSize(titlesizey);
  pro_mc_udsg->GetYaxis()->SetTitleOffset(titleoffsety);



  pro_data->GetXaxis()->SetLabelSize(labelsizex);
  pro_data->GetXaxis()->SetTitleSize(titlesizex);
  pro_data->GetXaxis()->SetTitleOffset(titleoffsetx);
  pro_mc->GetXaxis()->SetLabelSize(labelsizex);
//.........这里部分代码省略.........