int InputForLimits(){
  TH1::SetDefaultSumw2(true);
  if(pcp)cout<<"going to set inputs"<<endl;

  Int_t NBR = 3;
  Float_t BR[] = { 1., 0.75, 0.5};

  TFile* bkgFile = new TFile( "../../BkgPrediction/BkgPrediction.root", "READ"); 
  TTree* bkgTree;
  bkgFile->GetObject( "ElAndMu", bkgTree);
  
  Float_t bkg = 0.;
  Float_t bkgTotUnc = 0.;

  Float_t obs = 0.;

  bkgTree->SetBranchAddress( "srData", &obs);
  bkgTree->SetBranchAddress( "srAllBkgCorr", &bkg);
  bkgTree->SetBranchAddress( "TotUnc", &bkgTotUnc);

  TFile* sigFile = new TFile( "../../SignalSystematics/SignalSys.root", "READ"); 

  std::vector<std::vector<TString> > sysColl;
  std::vector<TString> sys;
  sys.push_back(TString("JES_Up"));
  sys.push_back(TString("JES_Down"));
  sysColl.push_back(sys);

  sys.clear();
  sys.push_back(TString("BTagReweight_UpBC"));
  sys.push_back(TString("BTagReweight_DownBC"));
  sysColl.push_back(sys);

  sys.clear();
  sys.push_back(TString("BTagReweight_UpLight"));
  sys.push_back(TString("BTagReweight_DownLight"));
  sysColl.push_back(sys);

  std::vector<TString> decayMode;
  decayMode.push_back(TString("tt"));
  decayMode.push_back(TString("tb"));
  decayMode.push_back(TString("bb"));
		      
  Systematics* systematics[3];
  TString dirname;
  TString histoname;

  TDirectory* srDir;
  TDirectory* histoDir;

  TFile* outFile = new TFile( "InputForLimits.root", "RECREATE");
  TDirectory* outBRDir;
  TDirectory* outSRDir;

  TH1F* datah = new TH1F( "data", "data", 1, 0., 1.);
  TH1F* bkgh = new TH1F( "bkg", "bkg", 1, 0., 1.);

  TH2F* sigh;
  TH2F* sig_toth;
  TH2F* effh;
  
  TH2F* sigLh;
  TH2F* effLh;

  TH2F* sigRh;
  TH2F* effRh;
  
  TH2F* jesh;
  TH2F* btagBCh;
  TH2F* btagLighth;
  TH2F* btagh;
  
  TH2F* sysh;

  TH2F* unch;

  TH2F* jesPercenth;
  TH2F* btagBCPercenth;
  TH2F* btagLightPercenth;
  TH2F* btagPercenth;
  
  TH2F* sysPercenth;

  TH2F* uncPercenth;

  Float_t sig = 0.;
  Float_t stat = 0.;

  Float_t jes = 0.;
  Float_t bc = 0.;
  Float_t light = 0.;
  
  Float_t unc = 0.;

  int N = bkgTree->GetEntries();

  for ( int ibr = 0; ibr < NBR; ibr++){
    dirname = ""; dirname += BR[ibr];
    outFile->mkdir(dirname);
    outBRDir = outFile->GetDirectory(dirname);
//.........这里部分代码省略.........

int evtmatching(TString infdfinder=Form("/export/d00/scratch/jwang/Dmeson/evtmatchingInput/DfinderMC_Pyquen_D0tokaonpion_D0pt%.0fp0_Pthat%.0f_TuneZ2_Unquenched_5020GeV_GENSIM_75x_v2_20151027.root",pthat,pthat), TString infhlt=Form("/export/d00/scratch/jwang/Dmeson/evtmatchingInput/openHLT_HF_HLTHeavyFlavour_MVA_Pyquen_D0tokaonpion_D0pt%.0fp0_Pthat%.0f_1026.root",pthat,pthat),TString outfile=Form("/export/d00/scratch/jwang/Dmeson/DfinderMC_20151028_EvtMatching_Pyquen_D0tokaonpion_D0pt%.0fp0_Pthat%.0f_TuneZ2_Unquenched_5020GeV_GENSIM_75x_v2_20151027.root",pthat,pthat))
{
  TFile* fdfinder = new TFile(infdfinder);
  TFile* fhlt = new TFile(infhlt);
  TTree* root = (TTree*)fdfinder->Get("Dfinder/root");
  TTree* hiroot = (TTree*)fdfinder->Get("hiEvtAnalyzer/HiTree");
  TTree* hltroot = (TTree*)fhlt->Get("hltbitanalysis/HltTree");

  TFile* outf = new TFile(outfile,"recreate");
  TDirectory* droot = outf->mkdir("Dfinder","");
  TDirectory* dhltroot = outf->mkdir("hltbitanalysis","");
  TDirectory* dhiroot = outf->mkdir("hiEvtAnalyzer","");
  droot->cd();
  TTree* ntReco = root->CloneTree(0);
  dhltroot->cd();
  TTree* ntHlt = hltroot->CloneTree(0);
  dhiroot->cd();
  TTree* ntHi = hiroot->CloneTree(0);

  setEvtDBranch(root);
  setEvtHLTBranch(hltroot);
  setEvtHIBranch(hiroot);
  cout<<"---- Check evt no. for three trees"<<endl;
  cout<<root->GetEntries()<<", "<<hltroot->GetEntries()<<", "<<hiroot->GetEntries()<<endl;
  Long64_t nentriesReco = root->GetEntries();
  Long64_t nentriesHlt = hltroot->GetEntries();
  Long64_t nentriesHi = hiroot->GetEntries();
  Int_t aRecoRunNo[nentriesReco],aRecoEvtNo[nentriesReco],aRecoLumiNo[nentriesReco];
  Int_t aHltRun[nentriesHlt],aHltEvent[nentriesHlt],aHltLumiBlock[nentriesHlt];
  Int_t aHiRun[nentriesHi],aHiEvt[nentriesHi],aHiLumi[nentriesHi];
  cout<<"---- Event reading"<<endl;
  Long64_t nentriesMax;
  if(nentriesHi>((nentriesReco>nentriesHlt)?nentriesReco:nentriesHlt)) nentriesMax=nentriesHi;
  else nentriesMax=((nentriesReco>nentriesHlt)?nentriesReco:nentriesHlt);
  for(Int_t i=0;i<nentriesMax;i++)
    {
      if(i%10000==0)
	{
	  cout<<setw(7)<<i<<" / "<<nentriesMax<<endl;
	}
      if(i<nentriesReco)
	{
	  root->GetEntry(i);
	  aRecoRunNo[i] = RecoRunNo;
	  aRecoEvtNo[i] = RecoEvtNo;
	  aRecoLumiNo[i] = RecoLumiNo;
	}
      if(i<nentriesHlt)
	{
	  hltroot->GetEntry(i);
	  aHltRun[i] = HltRun;
	  aHltEvent[i] = (Int_t)HltEvent;
	  aHltLumiBlock[i] = HltLumiBlock;
	}
      if(i<nentriesHi)
	{
	  hiroot->GetEntry(i);
	  aHiRun[i] = HiRun;
	  aHiEvt[i] = HiEvt;
	  aHiLumi[i] = HiLumi;
	}
    }

  Int_t hltmatching[nentriesReco],himatching[nentriesReco];
  ofstream fout("evtmatchingResult/evtmatching.dat");
  cout<<"---- Event matching"<<endl;
  for(Int_t ievt=0;ievt<nentriesReco;ievt++)
    {
      if(ievt%10000==0) cout<<setw(7)<<ievt<<" / "<<nentriesReco<<endl;
      hltmatching[ievt] = -1;
      himatching[ievt] = -1;
      for(Int_t jevt=0;jevt<nentriesHlt;jevt++)
	{
	  if(aRecoRunNo[ievt]==aHltRun[jevt]&&aRecoEvtNo[ievt]==aHltEvent[jevt]&&aRecoLumiNo[ievt]==aHltLumiBlock[jevt])
	    {
	      hltmatching[ievt] = jevt;
	      break;
	    }
	}
      for(Int_t jevt=0;jevt<nentriesHi;jevt++)
	{
	  if(aRecoRunNo[ievt]==aHiRun[jevt]&&aRecoEvtNo[ievt]==aHiEvt[jevt]&&aRecoLumiNo[ievt]==aHiLumi[jevt])
	    {
	      himatching[ievt] = jevt;
	      break;
	    }
	}
      fout<<hltmatching[ievt]<<" "<<himatching[ievt]<<endl;
    }
  cout<<"---- Event matching done"<<endl;
  cout<<"---- Writing hlt tree"<<endl;
  for(Int_t i=0;i<nentriesReco;i++)
    {
      if (i%1000==0) cout<<setw(7)<<i<<" / "<<nentriesReco<<"  HI index:"<<setiosflags(ios::left)<<setw(7)<<himatching[i]<<" HLT index:"<<hltmatching[i]<<endl;
      if(hltmatching[i]<0)
	{
	  continue;
	}
      root->GetEntry(i);
      droot->cd();
//.........这里部分代码省略.........

void NormalizeRunTo1Paddle(const Int_t run, const TString& histname)
{

    cout << fnm.GetSclrFileName(true) << endl;

    ScalerTool st(fnm.GetSclrFileName(true).data());

    stringstream ss;
    ss << fnm.GetHistFileBase(true) << run << ".root";
    TFile f(ss.str().data(),"UPDATE");
    cout << "Reading " << ss.str() << endl;
    TDirectory *dir = f.GetDirectory(histname.Data());
    if (dir==NULL)
    {
        cout << histname.Data() <<" is not a directory" << endl;
        return;
    }
    else
    {
        dir->cd();
    }
    vector<TH1*> vec = ROOTUtils::GetAllTH1InFile(dir);
    // the order is not important
    //ROOTUtils::EnforceProperOrdering(vec);
    cout << "Found " << vec.size() << " histograms" << endl;

    stringstream ss1;
    ss1 << normed_file << run << ".root";
    TFile normf(ss1.str().data(),"UPDATE");
    TDirectory *normdir = normf.GetDirectory(histname.Data());
    if (normdir==NULL)
    {
        normdir = normf.mkdir(histname.Data());
    }
    normdir->cd();

    vector<TH1*>::iterator it;
    Int_t val = st.GetValue(run,14);
    cout << endl;
    cout << "Normalization by 1Pad/100 = " << val << endl;
    for (it=vec.begin(); it!=vec.end(); it++)
    {
        if (val!=0)
            (*it)->Scale(0.01/val);
        else
            continue;

        string name = (*it)->GetName();
        name.append("_norm_1pad");
        (*it)->SetName(name.data());
        name = (*it)->GetTitle();
        name.append(" Normed to 1Pad");
        (*it)->SetTitle(name.data());
    //    normf.cd();
        normdir->cd();
        (*it)->Write();
    }
    cout << "Normalized histograms saved in " << normed_file << endl;

    //  return 0;

}

//.........这里部分代码省略.........
              }
              else if( tcmet_ > 30 ){
                fillUnderOverFlow( hmueta_metgt30 , hyp_lt_p4()[hypIdx].eta() , weight );
                fillUnderOverFlow( hmueta_metgt30 , hyp_ll_p4()[hypIdx].eta() , weight );
              }
            }

            //look at track pT's near electron
            if( myType == 0 && tcmet_ > 30 ){
              for( unsigned int itrk = 0 ; itrk < trks_trk_p4().size() ; ++itrk ){
               
                if( isMuon( itrk ) )            continue;
                if( isElectron( itrk ) )        continue;
                if( !isGoodTrack( itrk ) )      continue;
             
                LorentzVector vtrk = trks_trk_p4().at(itrk);
                LorentzVector vlt  = hyp_lt_p4()[hypIdx];
                LorentzVector vll  = hyp_ll_p4()[hypIdx];

                float drll = dRbetweenVectors(vtrk, vll);
                float drlt = dRbetweenVectors(vtrk, vlt);

                if( drll < 0.3 || drlt < 0.3 ){
                  fillUnderOverFlow( htrkptnearel , trks_trk_p4().at(itrk).pt() , weight );
                  //printEvent();
                  //cout << "pt eta phi " << vtrk.pt() << " " << vtrk.eta() << " " << vtrk.phi() << endl;
                }
              }
            }

            string leptype[2]={"ee","mm"};

            if( isData && ( tcmet_ > 30. || pfmet_ > 30. ) ){

              if( calculateTCMET ){
                
                // calculate tcmet on-the-fly
                metStruct structMET = correctedTCMET( true, ofile );
                
              }

              cout << "|" << setw(8)  << evt_run()                   << setw(4) 
                   << "|" << setw(6)  << evt_lumiBlock()             << setw(4) 
                   << "|" << setw(12) << evt_event()                 << setw(4) 
                   << "|" << setw(6)  << leptype[myType]             << setw(4) 
                   << "|" << setw(6)  << nJets                       << setw(4) 
                   << "|" << setw(6)  << nBTags                      << setw(4) 
                   << "|" << setw(8)  << fround(tcmet_,1)            << setw(4) 
                   << "|" << setw(8)  << fround(tcmetNew_calo_,1)    << setw(4) 
                   << "|" << setw(8)  << fround(tcmetNew_pfc_,1)     << setw(4) 
                   << "|" << setw(8)  << fround(pfmet_,1)            << setw(4) 
                   << "|" << setw(8)  << fround(dphijetmet_,2)       << setw(4) << "|" << endl; 
            }

            if( imax > -1 ){

              if( tcmet_ < 20 )
                fillHistos( hdphijetmet_metlt20  , dphijetmet_    , weight , myType , nJets  );

              if( tcmet_ > 30 )
                fillHistos( hdphijetmet_metgt30  , dphijetmet_    , weight , myType , nJets  );

            }

            fillHistos( hnjets            , nJets            , weight , myType );
            fillHistos( htcmet            , tcmet_           , weight , myType , nJets );
            fillHistos( htcmetNew_calo    , tcmetNew_calo_   , weight , myType , nJets  );
            fillHistos( htcmetNew_pfc     , tcmetNew_pfc_    , weight , myType , nJets  );
            fillHistos( hpfmet            , pfmet_           , weight , myType , nJets  );
            
            njets_   = nJets;
            leptype_ = myType;

            eventTree_->Fill();
          }// end loop over hypotheses
        } // end loop over events
    } // end loop over files
  
  cout << "\n\n********************SUMMARY********************" << endl;
  cout << "Total number of events: " << nEventsTotal << endl;
  cout << "Total number of events that pass dup veto: " << nPassDuplicate 
       << " (" << 100*(double)nPassDuplicate/nEventsTotal << "% of total)" << endl;
  cout << "Total number of events that pass good run/lumi: " << nPassGoodRun 
       << " (" << 100*(double)nPassGoodRun/nEventsTotal << "% of total)" << endl;

  cout << nGoodEl << " ee events in Z mass window" << endl;
  cout << nGoodMu << " mm events in Z mass window" << endl;

  if (nEventsChain != nEventsTotal)
    std::cout << "ERROR: number of events from files is not equal to total number of events" << std::endl;

  CloseBabyNtuple();

  TDirectory *rootdir = gDirectory->GetDirectory("Rint:");
  rootdir->cd();
  saveHist( Form( "output/%s_%s_histos.root" , prefix , iter ) );
  deleteHistos();

  
} // end ScanChain

//.........这里部分代码省略.........
        if(vetoedInputFile)
        {
            vetoedInputFile->Close();
        }

        nonVetoInputFile->Close();
        macropulseFile->Close();
        return 1;
    }

    for(int i=0; i<numberOfEntries; i++)
    {
        macropulseTree->GetEntry(i);

        macropulseList.push_back(me);
    }

    // open gamma correction file
    TFile* gammaCorrectionFile = new TFile(gammaCorrectionFileName.c_str(),"READ");
    if(!gammaCorrectionFile->IsOpen())
    {
        cerr << "Error: failed to open " << gammaCorrectionFileName << "  to read gamma correction." << endl;

        if(vetoedInputFile)
        {
            vetoedInputFile->Close();
        }

        nonVetoInputFile->Close();
        macropulseFile->Close();
        return 1;
    }

    TDirectory* gammaDirectory = (TDirectory*)gammaCorrectionFile->Get(config.analysis.GAMMA_CORRECTION_TREE_NAME.c_str());
    if(!gammaDirectory)
    {
        cerr << "Error: failed to open summedDet directory in " << gammaCorrectionFileName << " for reading gamma corrections." << endl;

        if(vetoedInputFile)
        {
            vetoedInputFile->Close();
        }

        nonVetoInputFile->Close();
        macropulseFile->Close();
        gammaCorrectionFile->Close();
        return 1;
    }

    gammaDirectory->cd();

    TH1D* gammaCorrectionHisto = (TH1D*)gammaDirectory->Get("gammaCorrection");
    if(!gammaCorrectionHisto)
    {
        cerr << "Error: failed to open gammaCorrections histo in " << gammaCorrectionFileName << " for reading gamma corrections." << endl;

        if(vetoedInputFile)
        {
            vetoedInputFile->Close();
        }

        nonVetoInputFile->Close();
        macropulseFile->Close();
        gammaCorrectionFile->Close();

        return 1;

//.........这里部分代码省略.........
      lnsq_40[i] = constant_mass(i*250,squarkMasses);
      lngl_40[i] = constant_mass(i*250,gluinoMasses);
      sq_40_text[i] = constant_squark_text_tanBeta40(i*250,lnsq_40[i]);
      gl_40_text[i] = constant_gluino_text_tanBeta40(i*250,lngl_40[i]);;
    }
  }
  

  //Legends
  TLegend* legst  = makeStauLegend(0.05,tanBeta_);
//   TLegend* legNoEWSB  = makeNoEWSBLegend(0.05,tanBeta_);
  TLegend* legexp = makeExpLegend( *TEV_sg_cdf,*TEV_sg_d0,*LEP_ch,*LEP_sl,*TEV_sn_d0_1,0.035,tanBeta_);
  
 
  //make Canvas
  TCanvas* cvsSys = new TCanvas("cvsnm","cvsnm",0,0,800,600);
  gStyle->SetOptTitle(0);
  cvsSys->SetFillColor(0);
  cvsSys->GetPad(0)->SetRightMargin(0.07);
  cvsSys->Range(-120.5298,26.16437,736.0927,750);
  //  cvsSys->Range(-50.5298,26.16437,736.0927,500);
  cvsSys->SetFillColor(0);
  cvsSys->SetBorderMode(0);
  cvsSys->GetPad(0)->SetBorderMode(0);
  cvsSys->GetPad(0)->SetBorderSize(2);
  cvsSys->GetPad(0)->SetLeftMargin(0.1407035);
  cvsSys->GetPad(0)->SetTopMargin(0.08);
  cvsSys->GetPad(0)->SetBottomMargin(0.13);

  cvsSys->SetTitle("tan#beta="+tanb);
 
  output->cd();
  
  TDirectory* curDir = gDirectory;
//   TFile* f = new TFile("limits_binc_ht1000_met350-contours.root");
//   TFile* f = new TFile("limits_msugraNLO_multibtag_ht1000_met250_m12_0-550_HN_comb-contours.root");
  TFile* f = ContourFile;
  TGraph* ra4VieObs = (TGraph*)f->Get("gObs");
  TGraph* ra4VieExpM2 = (TGraph*)f->Get("gExpMinus2");
  TGraph* ra4VieExpM1 = (TGraph*)f->Get("gExpMinus1");
  TGraph* ra4VieExp = (TGraph*)f->Get("gExpMedian");
  TGraph* ra4VieExpP1 = (TGraph*)f->Get("gExpPlus1");
  TGraph* ra4VieExpP2 = (TGraph*)f->Get("gExpPlus2");
  TGraph* ra4VieExpArea(0);

  TGraph* ra4VieObsThM(0);
  TGraph* ra4VieExpThM(0);
  if ( ContourFileM ) {
    ra4VieObsThM = (TGraph*)ContourFileM->Get("gObs");
    ra4VieExpThM = (TGraph*)ContourFileM->Get("gExpMedian");
  }
  TGraph* ra4VieObsThP(0);
  TGraph* ra4VieExpThP(0);
  if ( ContourFileP ) {
    ra4VieObsThP = (TGraph*)ContourFileP->Get("gObs");
    ra4VieExpThP = (TGraph*)ContourFileP->Get("gExpMedian");
  }

  curDir->cd();


  double m0min = 0;
  double m0max=1600;
//   double m0max=1800;
  double xscale = m0max-m0min;
  if (tanBeta_ == 50) m0min=200;

// input: - Input file (result from TMVA)
//        - use of TMVA plotting TStyle
void deviations( TString fin = "TMVAReg.root", 
                 HistType htype = MVAType, Bool_t showTarget, Bool_t useTMVAStyle = kTRUE )
{
   // set style and remove existing canvas'
   TMVAGlob::Initialize( useTMVAStyle );
   gStyle->SetNumberContours(999);

   // switches
   const Bool_t Save_Images     = kTRUE;

   // checks if file with name "fin" is already open, and if not opens one
   TFile* file = TMVAGlob::OpenFile( fin );  

   // define Canvas layout here!
   Int_t xPad = 1; // no of plots in x
   Int_t yPad = 1; // no of plots in y
   Int_t noPad = xPad * yPad ; 
   const Int_t width = 650;   // size of canvas

   // this defines how many canvases we need
   TCanvas* c[100];

   // counter variables
   Int_t countCanvas = 0;

   // search for the right histograms in full list of keys
   //    TList* methods = new TMap();

   TIter next(file->GetListOfKeys());
   TKey *key(0);   
   while ((key = (TKey*)next())) {

      if (!TString(key->GetName()).BeginsWith("Method_")) continue;
      if (!gROOT->GetClass(key->GetClassName())->InheritsFrom("TDirectory")) continue;

      TString methodName;
      TMVAGlob::GetMethodName(methodName,key);
      cout << "--- Plotting deviation for method: " << methodName << endl;

      TObjString *mN  = new TObjString( methodName );
      TDirectory* mDir = (TDirectory*)key->ReadObj();

      TList* jobNames = new TList();

      TIter keyIt(mDir->GetListOfKeys());
      TKey *titkey;
      while ((titkey = (TKey*)keyIt())) {

         if (!gROOT->GetClass(titkey->GetClassName())->InheritsFrom("TDirectory")) continue;

         TDirectory *titDir = (TDirectory *)titkey->ReadObj();

         TObjString *jN = new TObjString( titDir->GetName() );
         if (!jobNames->Contains( jN )) jobNames->Add( jN );
         else delete jN;
	    
         TString methodTitle;
         TMVAGlob::GetMethodTitle(methodTitle,titDir);

         TString hname = "MVA_" + methodTitle;
         TIter   dirKeyIt( titDir->GetListOfKeys() );
         TKey*   dirKey;

         Int_t countPlots = 0;
         while ((dirKey = (TKey*)dirKeyIt())){
            if (dirKey->ReadObj()->InheritsFrom("TH2F")) {
               TString s(dirKey->ReadObj()->GetName());
               if (s.Contains("_reg_") && 
                   ( (showTarget && s.Contains("_tgt")) || (!showTarget && !s.Contains("_tgt")) ) && 
                   s.Contains( (htype == CompareType ? "train" : "test" ))) {
                  c[countCanvas] = new TCanvas( Form("canvas%d", countCanvas+1), 
                                                Form( "Regression output deviation versus %s for method: %s",
                                                      (showTarget ? "target" : "input variables"), methodName.Data() ),
                                                countCanvas*50+100, (countCanvas+1)*20, width, (Int_t)width*0.72 ); 
                  c[countCanvas]->SetRightMargin(0.10); // leave space for border
                  TH1* h = (TH1*)dirKey->ReadObj();
                  h->SetTitle( Form("Output deviation for method: %s (%s sample)", 
                                    hname.Data(), (htype == CompareType ? "training" : "test" )) );
                  //                                    methodName.Data(), (htype == CompareType ? "training" : "test" )) );
                  h->Draw("colz");
                  TLine* l = new TLine( h->GetXaxis()->GetXmin(), 0, h->GetXaxis()->GetXmax(), 0 );
                  l->SetLineStyle(2);
                  l->Draw();

                  // update and print
                  cout << "plotting logo" << endl;
                  TMVAGlob::plot_logo(1.058);
                  c[countCanvas]->Update();

                  TString fname = Form( "plots/deviation_%s_%s_%s_c%i", 
                                        methodName.Data(), 
                                        (showTarget ? "target" : "vars"),
                                        (htype == CompareType ? "training" : "test" ), countPlots );
                  TMVAGlob::imgconv( c[countCanvas], fname );

                  countPlots++;
                  countCanvas++;
               }
//.........这里部分代码省略.........

void LeptonPreselectionCMG( PreselType type, RooWorkspace * w ) {
	const Options & opt = Options::getInstance(); 
	if (type == ELE)
		cout << "Running Electron Preselection :" << endl;
	else if (type == MU)
		cout << "Running Muon Preselection :" << endl;
	else if (type == EMU)
		cout << "Running Electron-Muon Preselection() ..." << endl;
	else if (type == PHOT)
		cout << "Running Photon Preselection :" << endl;

	string systVar;
	try {
		systVar = opt.checkStringOption("SYSTEMATIC_VAR");
	} catch (const std::string & exc) {
		cout << exc << endl;
	}
	if (systVar == "NONE")
		systVar.clear();

#ifdef CMSSWENV
	JetCorrectionUncertainty jecUnc("Summer13_V4_MC_Uncertainty_AK5PFchs.txt");
#endif

	string inputDir = opt.checkStringOption("INPUT_DIR");
	string outputDir = opt.checkStringOption("OUTPUT_DIR");
	string sampleName = opt.checkStringOption("SAMPLE_NAME");
	string inputFile = inputDir + '/' + sampleName + ".root";
	cout << "\tInput file: " << inputFile << endl;

	bool isSignal = opt.checkBoolOption("SIGNAL");
	TGraph * higgsW = 0;
	TGraph * higgsI = 0;
	if (isSignal) {
		double higgsM = opt.checkDoubleOption("HIGGS_MASS");
		if (higgsM >= 400) {
			string dirName = "H" + double2string(higgsM);
			bool isVBF = opt.checkBoolOption("VBF");
			string lshapeHistName = "cps";
			string intHistName = "nominal";
			
			if (systVar == "LSHAPE_UP") {
				intHistName = "up";
			} else if (systVar == "LSHAPE_DOWN") {
				intHistName = "down";
			}

			if (isVBF) {
				TFile weightFile("VBF_LineShapes.root");
				higgsW = (TGraph *) ( (TDirectory *) weightFile.Get(dirName.c_str()))->Get( lshapeHistName.c_str() )->Clone();

			} else {
				TFile weightFile("GG_LineShapes.root");
				higgsW = (TGraph *) ( (TDirectory *) weightFile.Get(dirName.c_str()))->Get( lshapeHistName.c_str() )->Clone();
				TFile interfFile("newwgts_interf.root");
				higgsI = (TGraph *) ( (TDirectory *) interfFile.Get(dirName.c_str()))->Get( intHistName.c_str() )->Clone();
			}
		}
	}

	TFile * file = new TFile( inputFile.c_str() );
	if (!file->IsOpen())
		throw string("ERROR: Can't open the file: " + inputFile + "!");
	TDirectory * dir = (TDirectory *) file->Get("dataAnalyzer");
	TH1D * nEvHisto = (TH1D *) dir->Get("cutflow");
	TH1D * puHisto = (TH1D *) dir->Get("pileup");
	TTree * tree = ( TTree * ) dir->Get( "data" );
	Event ev( tree );
	const int * runP = ev.getSVA<int>("run"); 
	const int * lumiP = ev.getSVA<int>("lumi"); 
	const int * eventP = ev.getSVA<int>("event"); 
	const bool * trigBits = ev.getAVA<bool>("t_bits");
	const int * trigPres = ev.getAVA<int>("t_prescale");
	const float * metPtA = ev.getAVA<float>("met_pt");
	const float * metPhiA = ev.getAVA<float>("met_phi");
	const float * rhoP = ev.getSVA<float>("rho");
	const float * rho25P = ev.getSVA<float>("rho25");
	const int * nvtxP = ev.getSVA<int>("nvtx"); 
	const int * niP = ev.getSVA<int>("ngenITpu"); 
	
#ifdef PRINTEVENTS
	string eventFileName;
	if (type == ELE)
		eventFileName = "events_ele.txt";
	else if (type == MU)
		eventFileName = "events_mu.txt";
	else if (type == EMU)
		eventFileName = "events_emu.txt";

	EventPrinter evPrint(ev, type, eventFileName);
	evPrint.readInEvents("diff.txt");
	evPrint.printElectrons();
	evPrint.printMuons();
	evPrint.printZboson();
	evPrint.printJets();
	evPrint.printHeader();
#endif

	string outputFile = outputDir + '/' + sampleName;

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

//_____________________________________________________________________________
Int_t ProofAux::GenerateTree(const char *fnt, Long64_t ent, TString &fn)
{
   // Generate the main tree for the 'friends' tutorial; the tree is called
   // 'Tmain', has 'ent' entries and is saved to file 'fnt'.
   // The full file path is returned in 'fn'.
   // Return 0 on success, -1 on error.

   Int_t rc = -1;

   // Check the filename
   fn = fnt;
   if (fn.IsNull()) {
      Error("GenerateTree", "file name undefined!");
      return rc;
   }
   TUrl uu(fn, kTRUE);
   if (!strcmp(uu.GetProtocol(), "file") && !fn.BeginsWith("/")) {
      // Local file with relative path: create under the data directory
      if (!gProofServ ||
          !(gProofServ->GetDataDir()) || strlen(gProofServ->GetDataDir()) <= 0) {
         Error("GenerateTree", "data directory undefined!");
         return rc;
      }
      // Insert data directory
      fn.Insert(0, TString::Format("%s/", gProofServ->GetDataDir()));
      // Make sure the directory exists
      TString dir = gSystem->DirName(fn);
      if (gSystem->AccessPathName(dir, kWritePermission)) {
         if (gSystem->mkdir(dir, kTRUE) != 0) {
            Error("GenerateTree", "problems creating directory %s to store the file", dir.Data());
            return rc;
         }
      }
   }

   // Create the file
   TDirectory* savedir = gDirectory;
   TFile *f = new TFile(fn, "RECREATE");
   if (!f || f->IsZombie()) {
      Error("GenerateTree", "problems opening file %s", fn.Data());
      return rc;
   }
   savedir->cd();
   rc = 0;

   // Create the tree
   TTree *T = new TTree("Tmain","Main tree for tutorial friends");
   T->SetDirectory(f);
   Int_t Run = 1;
   T->Branch("Run",&Run,"Run/I");
   Long64_t Event = 0;
   T->Branch("Event",&Event,"Event/L");
   Float_t x = 0., y = 0., z = 0.;
   T->Branch("x",&x,"x/F");
   T->Branch("y",&y,"y/F");
   T->Branch("z",&z,"z/F");
   TRandom r;
   for (Long64_t i = 0; i < ent; i++) {
      if (i > 0 && i%1000 == 0) Run++;
      Event = i;
      x = r.Gaus(10,1);
      y = r.Gaus(20,2);
      z = r.Landau(2,1);
      T->Fill();
   }
   T->Print();
   f->cd();
   T->Write();
   T->SetDirectory(0);
   f->Close();
   delete f;
   delete T;

   // Notify success
   Info("GenerateTree", "file '%s' successfully created", fn.Data());

   // Add to the list
   TString fds(fn);
   if (!strcmp(uu.GetProtocol(), "file")) {
      if (gSystem->Getenv("LOCALDATASERVER")) {
         if (strcmp(TUrl(gSystem->Getenv("LOCALDATASERVER"), kTRUE).GetProtocol(), "file"))
            fds.Insert(0, TString::Format("%s/", gSystem->Getenv("LOCALDATASERVER")));
      } else {
         fds.Insert(0, TString::Format("root://%s/", gSystem->HostName()));
      }
   }
   fMainList->Add(new TObjString(fds));

   // Done
   return rc;
}

//
// Plot_ECvsP - plot the EC/P vs P distributions with cuts
//
//                  fAna = output from eg2a DMS
//                  target = target name
//
void Plot_ECvsP(char *fAna="Ana.root", char *fitParams="ECvsP_Fit.dat", char *target)
{
    Int_t i, j;
    Int_t iSector;
    
    const Int_t NPARAM = 3;
    TF1 *fitMean[NSECTORS];
    TF1 *fitCut[NSECTORS];
    TH2D *h2D[NSECTORS];
    
    char funcName[50];
    
    Double_t a, b, c, d, f;
    Double_t parMean[5];
    Double_t parSigma[2];
    
    // open text file for fit parameters
    ifstream fin(fitParams);

    // data files contain the trees
    printf("Analyzing file %s\n",fAna);
    TFile *fm = new TFile(fAna,"READ");
    TDirectory *tmp = fm->GetDirectory("ElectronID");
    
    // Canvas to plot histogram
    TCanvas *c1 = new TCanvas("c1","c1",0,0,1200,800);
    c1->SetBorderMode(1);  //Bordermode (-1=down, 0 = no border, 1=up)
    c1->SetBorderSize(5);
    gStyle->SetOptStat(1111);
    c1->SetFillStyle(4000);
    c1->Divide(3,2);
    
    for(j=0; j<NSECTORS; j++){
        sprintf(hname,"ECtotP_VS_P_Sector%i",j+1);
        h2D[j] = (TH2D*)tmp->Get(hname);
        
        c1->cd(j+1);
        gPad->SetLeftMargin(Lmar);
        gPad->SetRightMargin(Rmar);
        gPad->SetFillColor(0);
        
        h2D[j]->GetXaxis()->CenterTitle();
        h2D[j]->GetYaxis()->CenterTitle();
        h2D[j]->GetYaxis()->SetTitleOffset(yoff);
        
        h2D[j]->SetAxisRange(0.0,3.0,"X");
        h2D[j]->Draw("colz");
        
        fin >> iSector >> a >> b >> c >> d >> f;
        parMean[0] = a;
        parMean[1] = b;
        parMean[2] = c;
        parMean[3] = d;
        parMean[4] = f;
        
        sprintf(funcName,"fitMean%i",j+1);
        fitMean[j] = new TF1(funcName,polFit,0.0,3.0,3);
        fitMean[j]->SetParameters(&parMean[0]);
        fitMean[j]->SetLineWidth(2);
        fitMean[j]->Draw("same");
        
        sprintf(funcName,"fitCut%i",j+1);
        fitCut[j] = new TF1(funcName,CutBelow,0.0,3.0,5);
        fitCut[j]->SetParameters(parMean);
        fitCut[j]->SetLineColor(4);
        fitCut[j]->SetLineWidth(2);
        fitCut[j]->Draw("same");
        
        sprintf(funcName,"fitCut%i",j+1);
        fitCut[j] = new TF1(funcName,CutAbove,0.0,3.0,5);
        fitCut[j]->SetParameters(parMean);
        fitCut[j]->SetLineColor(4);
        fitCut[j]->SetLineWidth(2);
        fitCut[j]->Draw("same");
    }
    sprintf(OutCan,"Plot_ECvsP_%s.gif",target);
    c1->Print(OutCan);
    sprintf(OutCan,"Plot_ECvsP_%s.eps",target);
    c1->Print(OutCan);
    
    fin.close();
}

void makeVNDet(){

  bool testrun          = 0;
  const int norder_     = 4;
  const int QnBinOrder_ = 2;
  const double vtxCut   = 15.;

  static const int ptBinMin  = 0;
  static const int ptBinMax  = nptbinsDefault-1;
  static const int etaBinMin = 0; //0;
  static const int etaBinMax = netabinsDefault-1;

  TFile * fAna;
  TTree * tree;
  double centval;
  double vtx;
  TH2D * sumw;
  TH2D * sumwqx;
  TH2D * sumwqy;
  TH2I * hMult;
  double qnHFx_EP[NumEPNames];
  double qnHFy_EP[NumEPNames];
  double sumET_EP[NumEPNames];

  TFile * fQNDet;
  TH1D * hqnHFDet_x[NumEPNames];
  TH1D * hqnHFDet_y[NumEPNames];

  TFile * fQN;
  TH1D * hqbins[NCENT][NEPSymm];

  TFile * fOut;
  TDirectory * SubEvt_0;
  TDirectory * SubEvt_1;
  TDirectory * FullEvt;

  TH2D * hVNDetX_0[NQN];
  TH2D * hVNDetY_0[NQN];
  TH2D * hVNDetX_1[NQN];
  TH2D * hVNDetY_1[NQN];
  TH2D * hVNDetX_full[NQN];
  TH2D * hVNDetY_full[NQN];

  double VNRawX_0[NCENT][NEPSymm][NQN];
  double VNRawY_0[NCENT][NEPSymm][NQN];
  double VNRawX_1[NCENT][NEPSymm][NQN];
  double VNRawY_1[NCENT][NEPSymm][NQN];
  double VNRawX_full[NCENT][NEPSymm][NQN];
  double VNRawY_full[NCENT][NEPSymm][NQN];

  double sumw_0[NCENT][NEPSymm][NQN];
  double sumw_1[NCENT][NEPSymm][NQN];
  double sumw_full[NCENT][NEPSymm][NQN];

  int evtMult_0[NCENT][NEPSymm][NQN];
  int evtMult_1[NCENT][NEPSymm][NQN];
  int evtMult_full[NCENT][NEPSymm][NQN];
 
  double VNDetX_0[NCENT][NEPSymm][NQN];
  double VNDetY_0[NCENT][NEPSymm][NQN];
  double VNDetX_1[NCENT][NEPSymm][NQN];
  double VNDetY_1[NCENT][NEPSymm][NQN];
  double VNDetX_full[NCENT][NEPSymm][NQN];
  double VNDetY_full[NCENT][NEPSymm][NQN];

  int Nevents[NCENT][NEPSymm][NQN];
  int NFails[NCENT][NEPSymm][NQN];

  //
  // MAIN
  //


  //-- Set up the analyzer objects
  fAna = new TFile("/rfs/jcastle/PbPb2015/PixelTracking_MB2/EbyETree_pixel_tight.root");

  tree = (TTree*) fAna->Get("ebyeana/tree");
  sumwqx = new TH2D(Form("sumwqx%i", norder_), Form("sumwqx%i", norder_), nptbinsDefault, ptbinsDefault, netabinsDefault, etabinsDefault);
  sumwqy = new TH2D(Form("sumwqy%i", norder_), Form("sumwqy%i", norder_), nptbinsDefault, ptbinsDefault, netabinsDefault, etabinsDefault);
  sumw   = new TH2D("sumw",                    "sumw",                    nptbinsDefault, ptbinsDefault, netabinsDefault, etabinsDefault);
  hMult  = new TH2I("hMult",                   "hMult",                   nptbinsDefault, ptbinsDefault, netabinsDefault, etabinsDefault);
    
  tree->SetBranchAddress("Cent",                    &centval);
  tree->SetBranchAddress("Vtx",                     &vtx);
  tree->SetBranchAddress("mult",                    &hMult);
  tree->SetBranchAddress(Form("sumwqx%i", norder_), &sumwqx);
  tree->SetBranchAddress(Form("sumwqy%i", norder_), &sumwqy);
  tree->SetBranchAddress("sumw",                    &sumw);
  tree->SetBranchAddress("qnHFx_EP",                &qnHFx_EP);
  tree->SetBranchAddress("qnHFy_EP",                &qnHFy_EP);
  tree->SetBranchAddress("sumET_EP",                &sumET_EP);

  //-- Get the QN Detector histograms                                                                                                                                                                                             
  fQNDet = new TFile( Form("../../../../../../v%i/eta2.4/AnalyzerResults/Q%iDet.root", QnBinOrder_, QnBinOrder_) );
  for(int iEP = 0; iEP < NumEPNames; iEP++){
    int EPbin  = EPSymmPartnerBin[iEP];
    if( EPbin != EPSymmBin ) continue;
    hqnHFDet_x[iEP] = (TH1D*) fQNDet->Get( Form("hqnHFDet_x_%s", EPNames[iEP].data()) );
    hqnHFDet_y[iEP] = (TH1D*) fQNDet->Get( Form("hqnHFDet_y_%s", EPNames[iEP].data()) );
  }
//.........这里部分代码省略.........

//
// Analyze_ECvsP - fit the EC/P vs P distributions
//
//                  fAna = output from eg2a DMS
//                  target = target name
//
void Analyze_ECvsP(char *fAna="Ana.root", char *target)
{
    Int_t i, j;
    
    const Int_t NPARAM = 3;
    TH1D *h1D[NSECTORS][NPARAM];
    TH2D *h2D[NSECTORS];
    
    char HistName2D[50];
    char strname[50];
    char *yname[NPARAM] = {"Amplitude","Mean","Sigma"};
    
    // open text file for the yields
    char OutFile[100];
    sprintf(OutFile,"ECvsP_Fit_%s.dat",target);
    ofstream fout(OutFile);
    
    TCanvas *can[NSECTORS]; // Canvas to plot histogram
    
    Double_t par[3]={1.0,1.0,1.0};
    TF1 *pol;
    TF1 *sig;
    
    // data files contain the trees
    printf("Analyzing file %s\n",fAna);
    TFile *fm = new TFile(fAna,"READ");
    TDirectory *tmp = fm->GetDirectory("ElectronID");
    
    for(j=0; j<NSECTORS; j++){
        sprintf(cname,"can%i",j+1);
        sprintf(cname,"Canvas, Sector %i",j+1);
        can[j] = new TCanvas(cname,ctitle,50*j,0,600,600);
        can[j]->SetBorderMode(1);  //Bordermode (-1=down, 0 = no border, 1=up)
        can[j]->SetBorderSize(5);
        gStyle->SetOptStat(1111);
        can[j]->SetFillStyle(4000);
        can[j]->Divide(2,2);
        
        sprintf(HistName2D,"ECtotP_VS_P_Sector%i",j+1);
        h2D[j] = (TH2D*)tmp->Get(HistName2D);

        can[j]->cd(1);
        gPad->SetLeftMargin(Lmar);
        gPad->SetRightMargin(Rmar);
        gPad->SetFillColor(0);

        h2D[j]->GetXaxis()->CenterTitle();
        h2D[j]->GetYaxis()->CenterTitle();
        h2D[j]->GetYaxis()->SetTitleOffset(yoff);
    
        h2D[j]->SetAxisRange(0.0,3.0,"X");
        h2D[j]->Draw("colz");
    
        h2D[j]->FitSlicesY(0,0,-1,20);
    
        fout<<j+1<<"\t";
        
        for(i=0; i<NPARAM; i++){
            can[j]->cd(i+2);
            gPad->SetLeftMargin(Lmar);
            gPad->SetRightMargin(Rmar);
            gPad->SetFillColor(0);

            sprintf(strname,"%s_%i",HistName2D,i);
            h1D[j][i] = (TH1D*)gDirectory->Get(strname);
            h1D[j][i]->GetXaxis()->CenterTitle();
            h1D[j][i]->GetYaxis()->CenterTitle();
            h1D[j][i]->GetYaxis()->SetTitle(yname[i]);
            h1D[j][i]->GetYaxis()->SetTitleOffset(yoff);
            h1D[j][i]->SetAxisRange(0.0,2.5,"X");
            if (i==0){
                h1D[j][i]->SetAxisRange(0.0,40.0,"Y");
            }
            if (i==1) {
                h1D[j][i]->SetAxisRange(0.15,0.35,"Y");
                pol = new TF1("pol","pol2",0.5,3.0);
                h1D[j][i]->Fit("pol","R");
                fout<<pol->GetParameter(0)<<"\t"<<pol->GetParameter(1)<<"\t"<<pol->GetParameter(2)<<"\t";
            }
            if (i==2) {
                par[0] = 10.0;
                par[1] = 0.1;
                par[2] = 0.1;
                h1D[j][i]->SetAxisRange(0.0,0.3,"Y");
                sig = new TF1("sig",SigmaFit,0.5,3.0,2);
                sig->SetParameters(par);
                h1D[j][i]->Fit("sig","R");
                fout<<sig->GetParameter(0)<<"\t"<<sig->GetParameter(1)<<endl;
            }
            h1D[j][i]->Draw();
        }
       	sprintf(OutCan,"Analyze_ECvsP_S%i_%s.gif",j+1,target);
        can[j]->Print(OutCan);
        sprintf(OutCan,"Analyze_ECvsP_S%i_%s.eps",j+1,target);
//.........这里部分代码省略.........

/*************************************************************************************
 * createMFOutfile: moves the MLB distributions into an output file for use in 
 *                  R. Nally's MassFit.C code.
 *  input:  the main() arguments array
 *  output: writes to an output file the histograms, in a MassFit.C-readable format 
 *
 *  Structure-wise: can be implemented into class easily.
 ***********************************/
void createMFOutfile(const char* argv[]) {
  TFile *f = new TFile(argv[2]);

  //create the output file we'd like to write histograms to
  TFile *output = new TFile(outfileName, "RECREATE");
  
  //get the filesystem information from the file
  f->cd();
  TList *alokDirs = (TList*) f->GetListOfKeys();

  //create a list of "All" histograms and initialize (TODO)
  TList *allHists = new TList;

  //loop through the directories in the input file
  for(int idir=0; alokDirs->At(idir-1) != alokDirs->Last(); idir++) {
    //if it's not mlb, we don't care
    if(!TString(alokDirs->At(idir)->GetName()).Contains("_Mlb")) continue;
    
    //get the file directory information, its histograms
    TDirectory *cDir  = (TDirectory*) f->Get(alokDirs->At(idir)->GetName());
    TList *alokHistos = (TList*)      cDir->GetListOfKeys();

    //loop through the histograms in the current directory 
    for(int ihisto=0; alokHistos->At(ihisto-1) != alokHistos->Last(); ihisto++) {
      // don't consider the graph objects
      if(TString(alokHistos->At(ihisto)->GetName()).Contains("Graph_")) continue; 

      // clone the histogram, give it a new name
      TString cloneName = formatName(alokHistos->At(ihisto)->GetName(),nominalWidth);
      TH1F *thisto = (TH1F*) cDir->Get(alokHistos->At(ihisto)->GetName());
      TH1F *tclone = (TH1F*) thisto->Clone(cloneName);

      // open the outfile and write the histo in
      output->cd();
      TList *outkeys = (TList*) output->GetListOfKeys();

      // if the histogram already exists, add thisto to the existing one
      // messy, but removes the need for magic
      for(int iout=0; outkeys->At(iout-1) != outkeys->Last(); iout++) {
        if(outkeys->At(iout)->GetName() == cloneName) {
          cout<<" - found another histogram in output with name "<<cloneName<<endl;
          TH1F *theHisto = (TH1F*) output->Get(cloneName);
          cout<<" - got the same histogram from the output file"<<endl;
          TH1F *tHclone  = (TH1F*) theHisto->Clone(cloneName);
          cout<<" - cloned the histogram"<<endl;

          cout<<" - adding in clone"<<endl;
          tclone->Add(tHclone);

          cout<<" - deleting the original histogram from the output file"<<endl;
          output->Delete(cloneName + TString(";1"));
          cout<<" - deleted thing from output file"<<endl;
          cout<<" - tclone looks like "<<tclone<<endl;
        }
      }
      cout<<" - writing the tclone to file"<<endl;
      tclone->Write();

      // reopen the input root file and clean up
      f->cd();
      delete thisto;
      delete tclone;
    }
  }

  f->Close();
  output->cd();
  output->Close();

  // if we want to interpolate, start making more histograms
  if(interpolate) {
    // Have to reopen the outfile because ROOT is odd
    TFile *output2 = new TFile(outfileName, "UPDATE");
    output2->cd();

    std::pair<double,TString> maxPair = *moreFiles.begin();
    double maxWidth = maxPair.first;

    // check that it makes sense to interpolate with our settings
    if(maxWidth <= nominalWidth) {
      cout<<"\n\nERROR: Max width is less than or equal to the nominal width! Exiting."<<endl;
      exit(EXIT_FAILURE);
    }

    // open the input file and loop through the relevant directories
    TFile *maxFile = new TFile(maxPair.second);
    for(int i=0; i<lepsSize; i++) {
      // change directory and get the name of the folder we want to access
      maxFile->cd();
      char a[128];
      sprintf(a, "mlbwa_%s_Mlb", leps[i]);
      TDirectory