C++ ExRootTreeReader类代码示例

void analyzer()
{

  TString processName = "ZJets";

  TFile* f = TFile::Open(Form("hist_%s.root", processName.Data()), "recreate");

  // Create chain of root trees
  TChain chain("DelphesMA5tune");
  //kisti
  //chain.Add("/cms/home/tjkim/fcnc/sample/ZToLL50-0Jet_sm-no_masses/events_*.root");
  //hep
  chain.Add("/Users/tjkim/Work/fcnc/samples/ZToLL50-0Jet_sm-no_masses/events_*.root");  

  // Create object of class ExRootTreeReader
  ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
  Long64_t numberOfEntries = treeReader->GetEntries();
  
  // Get pointers to branches used in this analysis
  TClonesArray *branchMuon = treeReader->UseBranch("DelphesMA5tuneMuon");
  
  // Book histograms
  TH1 *histDiMuonMass = new TH1F("dimuon_mass","Di-Muon Invariant Mass (GeV)",100, 50, 150);
 
  // Loop over all events

  for(Int_t entry = 0; entry < numberOfEntries; ++entry)
  {
    // Load selected branches with data from specified event
    treeReader->ReadEntry(entry);
 
    int nmuon = 0; 
    for( int i = 0; i < branchMuon->GetEntries(); i++)
    {
      Muon *muon = (Muon*) branchMuon->At(i);
      if( muon->PT <= 20 || abs(muon->Eta) >= 2.4 ) continue;
      nmuon = nmuon + 1 ;
    }

    if( nmuon >= 2){

      Muon *mu1 = (Muon*) branchMuon->At(0);
      Muon *mu2 = (Muon*) branchMuon->At(1);

      // Plot di-muon invariant mass 
      histDiMuonMass->Fill(((mu1->P4()) + (mu2->P4())).M());
    }
 
  }

  // Show resulting histograms
  histDiMuonMass->Write();

  f->Close();

}

plotDelphes3(TString inputName){

  gSystem->Load("libDelphes");

  // Create chain of root trees
  TChain chain("Delphes");
  chain.Add(inputName.Data());
  
  // Create object of class ExRootTreeReader
  ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
  Long64_t numberOfEntries = treeReader->GetEntries();
  
  // Get pointers to branches used in this analysis
  TClonesArray *b_Jet = treeReader->UseBranch("Jet");
  TClonesArray *b_Electron = treeReader->UseBranch("Electron");
  TClonesArray *b_Muon = treeReader->UseBranch("Muon");
  TClonesArray *b_MET = treeReader->UseBranch("MissingET");
  TClonesArray *b_Truth = treeReader->UseBranch("Particle");

  TString outputName(inputName);
  outputName="test.root";
  TFile* output=new TFile(outputName,"RECREATE");

  // Loop over all events
  for(Int_t entry = 0; entry < numberOfEntries; ++entry)
    {
      if(entry==maxEvents) break;
      if(entry%10000==0)cout<<entry<<endl;

      // Load selected branches with data from specified event
      treeReader->ReadEntry(entry);

      if(b_MET->GetEntries() > 0){
	MissingET* met = (MissingET*) b_MET->At(0);
	h_MET->Fill(met->MET);
      }

      int nJets=0, nBJets=0;
      for(int jetNo=0; jetNo<b_Jet->GetEntries(); jetNo++){
	Jet *jet=(Jet*) b_Jet->At(jetNo);
	if(jet->PT>jet_pT_min && fabs(jet->Eta)<jet_eta_max){
	  nJets++;
	  h_jet_pT->Fill(jet->PT);
	  h_jet_eta->Fill(jet->Eta);
	    
	  if(jet->BTag){
	    nBJets++;
	    h_bJet_pT->Fill(jet->PT);
	    h_bJet_eta->Fill(jet->Eta);
	  }
	}
      }
      h_jet_mult->Fill(nJets);
      h_bJet_mult->Fill(nBJets);

      int nElectrons=0;
      for(int electronNo=0; electronNo<b_Electron->GetEntries(); electronNo++){
	Electron *electron=(Electron*) b_Electron->At(electronNo);
	if(electron->PT>electron_pT_min && fabs(electron->Eta)<electron_eta_max){
	  nElectrons++;
	  h_electron_pT->Fill(electron->PT);
	  h_electron_eta->Fill(electron->Eta);
	}
      }
      h_electron_mult->Fill(nElectrons);

      int nMuons=0;
      for(int muonNo=0; muonNo<b_Muon->GetEntries(); muonNo++){
	Muon *muon=(Muon*) b_Muon->At(muonNo);
	if(muon->PT>muon_pT_min && fabs(muon->Eta)<muon_eta_max){
	  nMuons++;
	  h_muon_pT->Fill(muon->PT);
	  h_muon_eta->Fill(muon->Eta);
	}
      }
      h_muon_mult->Fill(nMuons);

      int nGenMuons=0, nTMRMuonsIDIso=0;
      for(int genParticleNo=0; genParticleNo<b_Truth->GetEntries(); genParticleNo++){
	GenParticle *particle=(GenParticle*) b_Truth->At(genParticleNo);
	if (particle->Status==3){
	  if (abs(particle->PID)==13){
	    nGenMuons++;
	    h_genMuon_pT->Fill(particle->PT);
	    h_genMuon_eta->Fill(fabs(particle->Eta));
	    h_genMuon_eta_pT->Fill(fabs(particle->Eta),particle->PT);
	      
	    for(int muonNo=0; muonNo<b_Muon->GetEntries(); muonNo++){
	      Muon *muon=(Muon*) b_Muon->At(muonNo);
	      if(truthMatch(1,muon->Eta,muon->Phi,1,particle->Eta,particle->Phi)){
		nTMRMuonsIDIso++;
		h_TMRMuon_ID_Iso_pT->Fill(muon->PT);
		h_TMRMuon_ID_Iso_eta->Fill(fabs(muon->Eta));
		h_TMRMuon_ID_Iso_eta_pT->Fill(fabs(muon->Eta),muon->PT);
		break;  //To prevent 2 reco mu per truth mu                                                                                                                    
	      }
	    }
	  }
	}
      }
//.........这里部分代码省略.........

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





	//----------------------------------------------------------------------------------------
	////////////////////////////////////////  My Data STructure  /////////////////////////////
	//----------------------------------------------------------------------------------------

	TH1F::SetDefaultSumw2(true);

	gSystem->Load("libExRootAnalysis");
	gSystem->Load("libDelphes");


	// Create chain of root trees
	TChain chain("Delphes");
	//chain.Add(InputFileName.c_str());
	char filename[1000];
	FILE *input;
	input = fopen(InputFileName.c_str(),"r");
	if (input != NULL)
	{
		// lets read each line and get the filename from it
		while (fscanf(input,"%s\n",filename) != EOF)
		{
			printf("%s\n",filename);
			chain.Add(filename);
		}
	}


	// Create object of class ExRootTreeReader
	ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
	Long64_t numberOfEntries = treeReader->GetEntries();

	// Get pointers to branches used in this analysis
	TClonesArray *branchJet    = treeReader->UseBranch("Jet");
 	TClonesArray *branchGenJet = treeReader->UseBranch("GenJet");
 	//TClonesArray *branchSHT    = treeReader->UseBranch("ScalarHT");


	//----------------------------------------------------------------------------------------
	/////////////////////////////////////  LOOP  Over the EVENTS  ////////////////////////////
	//----------------------------------------------------------------------------------------

	cout << "Set Branch Addresses" << endl;

	int decade = 0;
	unsigned entries = 0;

	if (NENTRIES==-1) entries=numberOfEntries;
	else entries=NENTRIES;

	cout << "Reading TREE: " << numberOfEntries << " events available, \n";
	cout << "\t\t" << NENTRIES << " of them will be analyzed." << endl;



//Jet *genjet;
//Jet *calojet;
vector<Jet *> mycalojet;
vector<Jet *> mygenjet;


//TLorentzVector jet_2[2];

//.........这里部分代码省略.........
  myTree_ -> Branch("JetMass", JetMass_, "Jet[nJet]/F");	
  myTree_ -> Branch("JetEta", JetEta_, "JetEta[nJet]/F");	
  myTree_ -> Branch("JetPhi", JetPhi_, "JetPhi[nJet]/F");	
  myTree_ -> Branch("JetNtrk", JetNtrk_, "JetNtrk[nJet]/F");


  myTree_ -> Branch("ngghQuark",   &ngghQuark_,     "ngghQuark/I")   ;
  myTree_ -> Branch("gghQuarkPt",   gghQuarkPt_,    "gghQuarkPt[ngghQuark]/F")   ;
  myTree_ -> Branch("gghQuarkEta" , gghQuarkEta_,    "gghQuarkEta[ngghQuark]/F")  ;
  myTree_ -> Branch("gghQuarkPhi",  gghQuarkPhi_,    "gghQuarkPhi[ngghQuark]/F")  ;
  myTree_ -> Branch("gghQuarkMass", gghQuarkMass_,   "gghQuarkMass[ngghQuark]/F") ;
  myTree_ -> Branch("gghQuarkStatus", gghQuarkStatus_,  "gghQuarkStatus[ngghQuark]/F");
  myTree_ -> Branch("gghQuarkPdgId", gghQuarkPdgId_,  "gghQuarkPdgId[ngghQuark]/I");
  myTree_ -> Branch("gghQuarkMothPdgId", gghQuarkMothPdgId_,  "gghQuarkMothPdgId[ngghQuark]/I");
  
  myTree_ -> Branch("nZ",   &nZ_,     "nZ/I")   ;
  myTree_ -> Branch("ZPt",   ZPt_,    "ZPt[nZ]/F")   ;
  myTree_ -> Branch("ZEta" , ZEta_,    "ZEta[nZ]/F")  ;
  myTree_ -> Branch("ZPhi",  ZPhi_,    "ZPhi[nZ]/F")  ;
  myTree_ -> Branch("ZMass", ZMass_,   "ZMass[nZ]/F") ;
  myTree_ -> Branch("ZStatus", ZStatus_,  "ZStatus[nZ]/F");
  myTree_ -> Branch("ZPdgId", ZPdgId_,  "ZPdgId[nZ]/I");
  //myTree_ -> Branch("ZMothPdgId", ZMothPdgId_,  "ZMothPdgId[nZ]/I");
 
  myTree_ -> Branch("nMu",   &nMu_,     "nMu/I")   ;
  myTree_ -> Branch("MuPt",   MuPt_,    "MuPt[nMu]/F")   ;
  myTree_ -> Branch("MuEta" , MuEta_,    "MuEta[nMu]/F")  ;
  myTree_ -> Branch("MuPhi",  MuPhi_,    "MuPhi[nMu]/F")  ;
  myTree_ -> Branch("MuMass", MuMass_,   "MuMass[nMu]/F") ;
  myTree_ -> Branch("MuStatus", MuStatus_,  "MuStatus[nMu]/F");
  myTree_ -> Branch("MuPdgId", MuPdgId_,  "MuPdgId[nMu]/I");
  //myTree_ -> Branch("MuMothPdgId", MuMothPdgId_,  "MuMothPdgId[nMu]/I");
  
   
  myTree_->Branch("npu", &npu_, "npu/I");
  myTree_->Branch("pu_zpos", &pu_zpos_, "pu_zpos[npu]/F");
  myTree_->Branch("pu_sumpt_lowpt", &pu_sumpt_lowpt_, "pu_sumpt_lowpt[npu]/F");
  myTree_->Branch("pu_sumpt_highpt", &pu_sumpt_highpt_, "pu_sumpt_highpt[npu]/F");
  myTree_->Branch("pu_ntrks_lowpt", &pu_ntrks_lowpt_, "pu_ntrks_lowpt[npu]/F");
  myTree_->Branch("pu_ntrks_highpt", &pu_ntrks_highpt_, "pu_ntrks_highpt[npu]/F");


  myTree_->Branch("vxMC",&vxMC_,"vxMC/F");
  myTree_->Branch("vyMC",&vyMC_,"vyMC/F");
  myTree_->Branch("vzMC",&vzMC_,"vzMC/F");
  myTree_->Branch("nvertex",&nvertex_,"nvertex/I");
  myTree_->Branch("vx",&vx_,"vx[nvertex]/F");
  myTree_->Branch("vy",&vy_,"vy[nvertex]/F");
  myTree_->Branch("vz",&vz_,"vz[nvertex]/F");
  myTree_->Branch("vntracks",&vntracks_,"vntracks[nvertex]/F");
  myTree_->Branch("vchi2",&vchi2_,"vchi2[nvertex]/F");
  myTree_->Branch("vndof",&vndof_,"vndof[nvertex]/F");

  myTree_->Branch("mH", &mH_, "mH/I");

  signal (SIGINT, signalHandler);
 
  //Create chain of root trees
  TChain chain("Delphes");
  for (const auto &file : fileList)
    chain.Add (file.c_str ());

  // Create object of class ExRootTreeReader
  ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
  Long64_t numberOfEntries = treeReader->GetEntries();
  cout << "================================================================================" << endl;
  cout << "processing " << numberOfEntries << " events..." << endl;
  cout << "================================================================================" << endl << endl;

  // Get pointers to branches used in this analysis
  unordered_map<string, TClonesArray *> branches;
  branches["Event"] = treeReader->UseBranch("Event");
  branches["Particle"] = treeReader->UseBranch("Particle");
  branches["VBFquarks"] = treeReader->UseBranch("VBFquarks");
  branches["UnsortedCluster"] = treeReader->UseBranch("UnsortedCluster");
  branches["NPU"] = treeReader->UseBranch("NPU");
  branches["Jet"] = treeReader->UseBranch("Jet");
  branches["JetTRK"] = treeReader->UseBranch("JetTRK");
  //branches["JetGenClosest"] = treeReader->UseBranch("JetGenClosest");
  //branches["JetGenBest"] = treeReader->UseBranch("JetGenBest");
  branches["JetNoCHS"] = treeReader->UseBranch("JetNoCHS");
  branches["BeamSpotParticle"] = treeReader->UseBranch("BeamSpotParticle");
  branches["Track"] = treeReader->UseBranch("Track");
  branches["Rho"] = treeReader->UseBranch("Rho");
  branches["GenJet"] = treeReader->UseBranch("GenJet");

  // Loop over all events
  for(Int_t entry = 0; entry < numberOfEntries && !interrupted; ++entry)
    {
      if (!(entry % REPORT_EVERY))
        cout << "processing event " << (entry + 1) << "..." << endl;
      treeReader->ReadEntry(entry);
      analyze (branches, myTree_);
    }

  myTree_->Write();
  fout->Close ();

  delete treeReader;
}

void vbf_bbbb_select(const TString confname="test_vbf.txt",
	       const Float_t xsec=1.0,
	       const Float_t totalEvents=100,
	       const TString outputfile="test.root") {

  // declare constants
  const Int_t B_ID_CODE = 5;
  const Int_t H_ID_CODE = 25;

  const Float_t MAX_MATCH_DIST = 0.5;

 // read input input file
  TChain chain("Delphes");
  TString inputfile;
  
  ifstream ifs;
  ifs.open(confname);
  assert(ifs.is_open());
  string line;
  while (getline(ifs, line)) {
    stringstream ss(line);;
    ss >> inputfile;
    chain.Add(inputfile);    
  }
  ifs.close();
  
  ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
  Long64_t numberOfEntries = treeReader->GetEntries();

  TClonesArray *branchJet = treeReader->UseBranch("Jet");
  TClonesArray *branchGenJet = treeReader->UseBranch("GenJet");
  TClonesArray *branchParticle = treeReader->UseBranch("Particle");

  //set up loop variables
  GenParticle *genParticle=0;
  Jet *genJet=0;
  Jet *jet=0;

  // set up storage variables
  Jet         *jetB1=0, *jetB2=0, *jetB3=0, *jetB4=0;
  Jet         *jet1=0, *jet2=0;
  GenParticle *genB1=0, *genB2=0, *genB3=0, *genB4=0;
  GenParticle *genH1=0, *genH2=0;
  Jet *genJetB1=0, *genJetB2=0, *genJetB3=0, *genJetB4=0;
  Jet *genJetVBF1=0, *genJetVBF2=0;

  Int_t iB1=-1,       iB2=-1,       iB3=-1,       iB4=-1;
  Int_t iH1=-1,       iH2=-1;
  Int_t iJet1=-1,     iJet2=-1;
  Int_t iGenB1=-1,    iGenB2=-1,    iGenB3=-1,    iGenB4=-1;
  Int_t iGenJetB1=-1, iGenJetB2=-1, iGenJetB3=-1, iGenJetB4=-1;
  Int_t iGenJetVBF1=-1, iGenJetVBF2=-1;
  Int_t iHmatch1=-1,   iHmatch2=-1,  iHmatch3=-1,  iHmatch4=-1;

  LorentzVector *sRecoB1=0, *sRecoB2=0, *sRecoB3=0, *sRecoB4=0;
  LorentzVector *sGenJetB1=0, *sGenJetB2=0, *sGenJetB3=0, *sGenJetB4=0;
  LorentzVector *sGenJetVBF1=0, *sGenJetVBF2=0;
  LorentzVector *sGenB1=0, *sGenB2=0, *sGenB3=0, *sGenB4=0;
  LorentzVector *sGenH1=0, *sGenH2=0;
  LorentzVector *sRecoJet1=0, *sRecoJet2=0;

  TFile *outFile = new TFile(outputfile, "RECREATE");

  // tree to hold information about selected events
  TTree *outTree = new TTree("Events", "Events");

  outTree->Branch("iHmatch1",       &iHmatch1,                                                     "iHmatch1/i");   // which Higgs does b-jet 1 come from
  outTree->Branch("iHmatch2",       &iHmatch2,                                                     "iHmatch2/i");   // which Higgs does b-jet 2 come from
  outTree->Branch("iHmatch3",       &iHmatch3,                                                     "iHmatch3/i");   // which Higgs does b-jet 3 come from
  outTree->Branch("iHmatch4",       &iHmatch4,                                                     "iHmatch4/i");   // which Higgs does b-jet 4 come from

  outTree->Branch("sGenB1",         "ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double> >", &sGenB1);        // 4-vector for generator leading b-jet
  outTree->Branch("sGenB2",         "ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double> >", &sGenB2);        // 4-vector for generator b-jet
  outTree->Branch("sGenB3",         "ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double> >", &sGenB3);        // 4-vector for generator b-jet
  outTree->Branch("sGenB4",         "ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double> >", &sGenB4);        // 4-vector for generator b-jet

  outTree->Branch("sGenH1",         "ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double> >", &sGenH1);        // 4-vector for generator higgs
  outTree->Branch("sGenH2",         "ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double> >", &sGenH2);        // 4-vector for generator higgs

  outTree->Branch("sGenJetB1",      "ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double> >", &sGenJetB1);     // 4-vector for generator leading b-jet
  outTree->Branch("sGenJetB2",      "ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double> >", &sGenJetB2);     // 4-vector for generator b-jet
  outTree->Branch("sGenJetB3",      "ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double> >", &sGenJetB3);     // 4-vector for generator b-jet
  outTree->Branch("sGenJetB4",      "ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double> >", &sGenJetB4);     // 4-vector for generator b-jet

  outTree->Branch("sGenJetVBF1",    "ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double> >", &sGenJetVBF1);   // 4-vector for generator leading b-jet
  outTree->Branch("sGenJetVBF2",    "ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double> >", &sGenJetVBF2);   // 4-vector for generator b-jet

  outTree->Branch("sRecoB1",        "ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double> >", &sRecoB1);       // 4-vector for reconstructed leading b-jet
  outTree->Branch("sRecoB2",        "ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double> >", &sRecoB2);       // 4-vector for reconstructed b-jet
  outTree->Branch("sRecoB3",        "ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double> >", &sRecoB3);       // 4-vector for reconstructed b-jet
  outTree->Branch("sRecoB4",        "ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double> >", &sRecoB4);       // 4-vector for reconstructed b-jet

  outTree->Branch("sRecoJet1",      "ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double> >", &sRecoJet1);     // 4-vector for reconstructed leading forward-jet
  outTree->Branch("sRecoJet2",      "ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<double> >", &sRecoJet2);     // 4-vector for reconstructed second forward-jet

  // define placeholder vector for things that don't exist
  LorentzVector nothing(-999,-999,0,-999);

  Int_t iMatched=0;
  Int_t iNot=0;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    chain.Add("/afs/cern.ch/work/g/gkarapin/sample/400.root");
    chain.Add("/afs/cern.ch/work/g/gkarapin/sample/401.root");
    chain.Add("/afs/cern.ch/work/g/gkarapin/sample/402.root");
    chain.Add("/afs/cern.ch/work/g/gkarapin/sample/403.root");
    chain.Add("/afs/cern.ch/work/g/gkarapin/sample/405.root");
    chain.Add("/afs/cern.ch/work/g/gkarapin/sample/407.root");
    chain.Add("/afs/cern.ch/work/g/gkarapin/sample/409.root");
    chain.Add("/afs/cern.ch/work/g/gkarapin/sample/410.root");
    chain.Add("/afs/cern.ch/work/g/gkarapin/sample/411.root");
    chain.Add("/afs/cern.ch/work/g/gkarapin/sample/412.root");
    chain.Add("/afs/cern.ch/work/g/gkarapin/sample/413.root");
    chain.Add("/afs/cern.ch/work/g/gkarapin/sample/414.root");
    chain.Add("/afs/cern.ch/work/g/gkarapin/sample/415.root");
    chain.Add("/afs/cern.ch/work/g/gkarapin/sample/416.root");
    chain.Add("/afs/cern.ch/work/g/gkarapin/sample/417.root");
    chain.Add("/afs/cern.ch/work/g/gkarapin/sample/419.root");
    chain.Add("/afs/cern.ch/work/g/gkarapin/sample/420.root");
    chain.Add("/afs/cern.ch/work/g/gkarapin/sample/421.root");
    chain.Add("/afs/cern.ch/work/g/gkarapin/sample/424.root");
    chain.Add("/afs/cern.ch/work/g/gkarapin/sample/425.root");
    chain.Add("/afs/cern.ch/work/g/gkarapin/sample/426.root");
    chain.Add("/afs/cern.ch/work/g/gkarapin/sample/427.root");
    chain.Add("/afs/cern.ch/work/g/gkarapin/sample/428.root");
    chain.Add("/afs/cern.ch/work/g/gkarapin/sample/429.root");
    chain.Add("/afs/cern.ch/work/g/gkarapin/sample/430.root");
    chain.Add("/afs/cern.ch/work/g/gkarapin/sample/431.root");
    chain.Add("/afs/cern.ch/work/g/gkarapin/sample/432.root");
    //chain.Add("/afs/cern.ch/work/g/gkarapin/sample/433.root");
    chain.Add("/afs/cern.ch/work/g/gkarapin/sample/434.root");
    chain.Add("/afs/cern.ch/work/g/gkarapin/sample/435.root");
     //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    
    // Create object of class ExRootTreeReader
    ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
    Long64_t nentries = treeReader->GetEntries();

    std::cout << "can you read my tree ?" << std::endl;

    // Get pointers to branches used in this analysis
    TClonesArray *branchElectron = treeReader->UseBranch("Electron");
    TClonesArray *branchMissingET = treeReader->UseBranch("MissingET");
    TClonesArray *branchNPU = treeReader->UseBranch("NPU");
    
    Long64_t nentries = treeReader->GetEntries();
    
    const int qTnbins = 35;
    double qTxbins[qTnbins+1] = {0., 2, 4., 6, 8., 10, 12.5, 15, 17.5, 20., 22.5, 25.,28, 31,34,37, 40,43.5, 47, 51,55, 60, 65,70,76,82,89,97,106,116,127,139,152, 166,181,200}
    
    //Define histograms
    
    TH1D * histElectronEta = new TH1D("histElectronEta",";Electron #eta; number of events", 60, -3, 3);
    TH1D * histElectron1Pt = new TH1D("histElectron1Pt",";Electron1 P_{T} [GeV]; number of events", 200, 0.0, 200.0); 
    TH1D * histElectron2Pt = new TH1D("histElectron2Pt",";Electron2 P_{T} [GeV]; number of events", 200, 0.0, 200.0); 
    TH1D * histZmass   = new TH1D("histZmass","; Invariant mass of Z [GeV]; number of events", 120, 60.0, 120.0);
    TH1D * histZpT  = new TH1D("histZpT",";transverse momentum of Z, q_{T} [GeV]; number of events", 200, 0.0, 400.0);
    TH1D * histnvtx  = new TH1D("histnvtx",";number of vertices; number of events", 200, 0.0, 200);
    
    TH1D * histMet  = new TH1D("histMet",";#slash{E}_{T} [GeV]; number of events", 200, 0.0, 800.0);
    TH1D * histMetx = new TH1D("histMetx",";#slash{E}_{T}x [GeV]; number of events", 200, -400.0, 400.0);
    TH1D * histMety = new TH1D("histMety",";#slash{E}_{T}y [GeV]; number of events", 200, -400.0, 400.0); 
    
    TH1D * histuperp = new TH1D("histuperp",";u_{#perp}  [GeV]; number of events",100,-400.,400.);
    TH1D * histux = new TH1D("histux",";u_{x} [GeV]; number of events",100,-400.,400.);
    TH1D * histuy = new TH1D("histuy",";u_{y} [GeV]; number of events",100,-400.,400.); 
    TH1D * histupara = new TH1D("histupara",";u_{#parallel}  [GeV]; number of events",100,-400,+400); 
    TH1D * histuparaqT = new TH1D("histuparaqT","; -(u_{#parallel}+q_{T}) [GeV]; number of events",100,-400,+400);

int main(int argc, char*argv[])
{
 gSystem->Load("libDelphes");
 gSystem->Load("libExRootAnalysis");
 //arg
 string inputFile=argv[1];
 string outputFile=argv[2];
 TFile *out = TFile::Open(outputFile.c_str(),"RECREATE");
 // Create chain of root trees
 TChain chain("Delphes");
 chain.Add(inputFile.c_str());

 // Create object of class ExRootTreeReader
 ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
 Long64_t numberOfEntries = treeReader->GetEntries();

 // Get pointers to branches used in this analysis
 TClonesArray *branchJet  = treeReader->UseBranch("Jet");
 TClonesArray *branchGenJet = treeReader->UseBranch("GenJet");
 TClonesArray *scalarht = treeReader->UseBranch("ScalarHT");

 /////////////////////////////
 //Histograms/////////////////
 /////////////////////////////

 TH2 *hist_deltaY_60_100 = new TH2F("hist_deltaY_60_100","hist_deltaY_60_100",500,-5.0,5.0,500,-0.2,0.2);
 TH2 *hist_deltaY_100_150 = new TH2F("hist_deltaY_100_150","hist_deltaY_100_150",500,-5.0,5.0,500,-0.2,0.2);
 TH2 *hist_deltaY_150_200 = new TH2F("hist_deltaY_150_200","hist_deltaY_150_200",500,-5.0,5.0,500,-0.2,0.2);

 TH1 *histevent60100 = new TH1F("histevent60100","histevent60100",500,-0.5,0.5);
 TH1 *histevent100150 = new TH1F("histevent100150","histevent60100",500,-0.5,0.5);
 TH1 *histevent150200 = new TH1F("histevent150200","histevent60100",500,-0.5,0.5);

 TH1 *hist_PT_resolution = new TH1F("hist_PT_resolution","hist_PT_resolution",500,-1.0,1);

 TH2 *hist_PT_resolution_for_each[3];
 hist_PT_resolution_for_each[0] = new TH2F("hist_leading_pt_res","hist_PT_resolution_for_each",500,0.0,4000,500,-1.0,1.0);
 hist_PT_resolution_for_each[1] = new TH2F("hist_nleading_pt_res","hist_PT_resolution_for_each",500,0.0,4000,500,-1.0,1.0);
 hist_PT_resolution_for_each[2] = new TH2F("hist_nnleading_pt_res","hist_PT_resolution_for_each",500,0.0,4000,500,-1.0,1.0);

 TH1 *histleadingcolojetpt = new TH1F("leadingcolojetpt", "leadingcolojetpt", 50, 0.0, 400.0);
 TH1 *histnleadingcolojetpt = new TH1F("nleadingcolojetpt", "nleadingcolojetpt", 50, 0.0, 400.0);
 TH1 *histnnleadingcolojetpt = new TH1F("nnleadingcolojetpt", "nnleadingcolojetpt", 50, 0.0, 400.0);

 TH1 *histcolojetht = new TH1F("histcolojetht","histcolojetht",1000, 0.0, 10000.0);




 TH2 *histresleadingjet = new TH2F("histresleadingjet","histresleadingjet",500,0.0,4000,500,-1.0,1.0);
 TH2 *histres2leadingjet = new TH2F("histres2leadingjet","histresleadingjet",500,0.0,4000,500,-1.0,1.0);
 TH2 *histres3leadingjet = new TH2F("histres3leadingjet","histresleadingjet",500,0.0,2500,500,-1.0,1.0);

 TH1 *histcolo2jetht = new TH1F("histcolo2jetht",">2jet ht",500,0.0,10000.0);
 TH1 *histcolo3jetht = new TH1F("histcolo3jetht",">3jet ht",500,0.0,10000.0);

 TH1 *histgen2jetht = new TH1F("histgen2jetht",">2jet ht",500,0.0,10000.0);
 TH1 *histgen3jetht = new TH1F("histgen3jetht",">3jet ht",500,0.0,10000.0);

  //Definition
 vector<Jet *> colojet;
 vector<Jet *> genjet;

 ScalarHT* scht=0;
 Jet *cjet;
 Jet *gjet;

 int nocolojet = 0;
 //double ptrescut = 0;


  for(int entry = 0; entry < numberOfEntries; ++entry)
  {
     //her 10000 olayda bıze bılgı verıyor
     if(entry%10000 == 0) cout << "event number: " << entry << endl;
     // Load selected branches with data from specified event
     treeReader->ReadEntry(entry);

     //////////////////////////////////////
     ///Scalerht yı scht ye dolduruyor.////
     //////////////////////////////////////

     for(int i =0 ;i<scalarht->GetEntriesFast();++i )
     {
       scht =(ScalarHT*) scalarht->At(i);
     }

     ///////////////////////////////////////////////
     ////genjet bilgilerini alıyor//////////////////
     ///////////////////////////////////////////////

     genjet.clear();
     for(int i=0;i < branchGenJet->GetEntriesFast();++i)
     {
       gjet = (Jet*) branchGenJet->At(i);
       genjet.push_back(gjet);
     }

     //////////////////////////////////////////////
     ///colojet bilgilerini alıyor ////////////////
//.........这里部分代码省略.........

//void LeptonEnergy(const char *inputFile = "sourceFiles/LO/ttbar_LO_total.root")
void LeptonEnergy(const TString & file)
{

  //const char *inputFile = Form("/home/tjkim/work/pheno/topmass/sourceFiles/LO/fromSayaka/ttbar_%s.root",file.Data());
  //gSystem->Load("/export/apps/delphes//libDelphes");
  const char *inputFile = Form("/data/users/seohyun/analysis/ttbar_%s.root",file.Data());
  gSystem->Load("/home/seohyun/delphes/libDelphes.so");
/*
  TFile *f2 = TFile::Open("weightfunc2.root");
  TFile *f3 = TFile::Open("weightfunc3.root");
  TFile *f5 = TFile::Open("weightfunc5.root");
  TFile *f15 = TFile::Open("weightfunc15.root");

  const int nmass = 151;
  float mymass[ nmass ];
  float integral2[ nmass ];
  float integral3[ nmass ];
  float integral5[ nmass ];
  float integral15[ nmass ];
  for(int i=0; i < nmass ; i++){
    integral2[i] = 0.0;
    integral3[i] = 0.0;
    integral5[i] = 0.0;
    integral15[i] = 0.0;
  }

  TGraph * g2[nmass];
  TGraph * g3[nmass];
  TGraph * g5[nmass];
  TGraph * g15[nmass];

  TIter next(f2->GetListOfKeys());
  TKey *key;
  int i = 0;
  while( (key = (TKey*) next())) {
    TClass *cl = gROOT->GetClass( key->GetClassName());
    if( !cl->InheritsFrom("TGraph")) continue;
    g2[i] = (TGraph*) key->ReadObj();
    string mass = g2[i]->GetName(); 
    float temp = atof(mass.c_str()); 
    mymass[i] = temp;
    i++;
  }

  TIter next(f3->GetListOfKeys());
  i = 0;
  while( (key = (TKey*) next())) {
    TClass *cl = gROOT->GetClass( key->GetClassName());
    if( !cl->InheritsFrom("TGraph")) continue;
    g3[i] = (TGraph*) key->ReadObj();
    i++;
  }

  TIter next(f5->GetListOfKeys());
  i = 0;
  while( (key = (TKey*) next())) {
    TClass *cl = gROOT->GetClass( key->GetClassName());
    if( !cl->InheritsFrom("TGraph")) continue;
    g5[i] = (TGraph*) key->ReadObj();
    i++;
  }

  TIter next(f15->GetListOfKeys());
  i = 0;
  while( (key = (TKey*) next())) {
    TClass *cl = gROOT->GetClass( key->GetClassName());
    if( !cl->InheritsFrom("TGraph")) continue;
    g15[i] = (TGraph*) key->ReadObj();
    i++;
  } 

  TFile * res = TFile::Open("hist_LO_res_v3.root");
  TH1F * h_acc = (TH1F*) res->Get("h_totalacc_lepton"); 
*/

  //TFile* f = TFile::Open("hist_LO_res_60.root", "recreate");
  TFile* f = TFile::Open(Form("170717/hist_%s.root",file.Data()), "recreate");

  // Create chain of root trees
  TChain chain("Delphes");
  chain.Add(inputFile);
 
  // Create object of class ExRootTreeReader
  ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
  Long64_t numberOfEntries = treeReader->GetEntries();

  // Get pointers to branches used in this analysis
  TClonesArray *branchParticle = treeReader->UseBranch("Particle");
  TClonesArray *branchMuon = treeReader->UseBranch("Muon");
  TClonesArray *branchElectron = treeReader->UseBranch("Electron");
  TClonesArray *branchJet = treeReader->UseBranch("Jet");

  TClonesArray *branchEvent = treeReader->UseBranch("Event");
 
  GenParticle *particle;
  GenParticle *daughter1;
  GenParticle *daughter2;
  GenParticle *granddaughter1_1;
  GenParticle *granddaughter1_2;
//.........这里部分代码省略.........

void whbb_trigger(const TString input="whbb.txt",
		  const TString outputfile="/afs/cern.ch/work/j/jlawhorn/whbb.root") {

  // declare constants
  //const Double_t MUON_MASS = 0.105658369;
  //const Double_t ELE_MASS  = 0.000511;
  //const Double_t TAU_MASS  = 1.77682;
  
  const Int_t GAM_ID_CODE = 22;
  const Int_t TAU_ID_CODE = 15;
  const Int_t MU_ID_CODE = 13;
  const Int_t ELE_ID_CODE = 11;
  
  const Float_t MAX_MATCH_DIST = 0.4;
  
  TChain chain("Delphes");
  
  ifstream ifs;
  ifs.open(input.Data()); assert(ifs.is_open());
  string line;
  while(getline(ifs,line)) {
    TString inputfile;
    stringstream ss(line);
    ss >> inputfile;
    chain.Add(inputfile);
  }
  ifs.close();
  
  ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
  Long64_t n = treeReader->GetEntries();
  
  //TClonesArray *branchJet = treeReader->UseBranch("Jet");
  
  //if (!(branchJet)) {
  //cout << "File broken" << endl;
  //return;
  //}
  
  //TClonesArray *branchElectron = treeReader->UseBranch("Electron");
  //TClonesArray *branchMuon = treeReader->UseBranch("Muon");
  TClonesArray *branchGenJet = treeReader->UseBranch("GenJet");
  TClonesArray *branchGenMET = treeReader->UseBranch("GenMissingET");
  TClonesArray *branchHT = treeReader->UseBranch("ScalarHT");
  TClonesArray *branchParticle = treeReader->UseBranch("Particle");

  GenParticle *part=0, *mama=0; Jet *jet=0;
  MissingET *metObj=0;
  ScalarHT *htObj=0;
  //Electron *ele=0; Muon *mu=0;

  UInt_t isReco=0;
  Float_t e1pt, e1eta, e1phi;
  Float_t e2pt, e2eta, e2phi;
  Float_t m1pt, m1eta, m1phi;
  Float_t m2pt, m2eta, m2phi;
  Float_t g1pt, g1eta, g1phi;
  Float_t g2pt, g2eta, g2phi;
  Float_t t1pt, t1eta, t1phi;
  Float_t t2pt, t2eta, t2phi;
  Float_t j1pt, j1eta, j1phi;
  Float_t j2pt, j2eta, j2phi;
  Float_t j3pt, j3eta, j3phi;
  Float_t j4pt, j4eta, j4phi;
  Float_t mht, ht;
  UInt_t triggerBits126;
  UInt_t triggerBits180;
  UInt_t triggerBits250;
  UInt_t triggerBits350;
  
  TFile *outfile = new TFile(outputfile, "RECREATE");
  TTree *outtree = new TTree("Events", "Events");
  Trigger::Event data;
  outtree->Branch("Events", &data.isReco, "isReco/i:e1pt/F:e1eta:e1phi:e2pt:e2eta:e2phi:m1pt:m1eta:m1phi:m2pt:m2eta:m2phi:g1pt:g1eta:g1phi:g2pt:g2eta:g2phi:t1pt:t1eta:t1phi:t2pt:t2eta:t2phi:j1pt:j1eta:j1phi:j2pt:j2eta:j2phi:j3pt:j3eta:j3phi:j4pt:j4eta:j4phi:mht:ht:triggerBits126/i:triggerBits180:triggerBits250:triggerBits350");

  //n=50;
  
  for (Int_t iEntry=0; iEntry<n; iEntry++) { // entry loop                                                                          
    treeReader->ReadEntry(iEntry);
    
    e1pt=-99; e1eta=-99; e1phi=-99; 
    e2pt=-99; e2eta=-99; e2phi=-99;
    m1pt=-99; m1eta=-99; m1phi=-99; 
    m2pt=-99; m2eta=-99; m2phi=-99; 
    g1pt=-99; g1eta=-99; g1phi=-99; 
    g2pt=-99; g2eta=-99; g2phi=-99; 
    t1pt=-99; t1eta=-99; t1phi=-99; 
    t2pt=-99; t2eta=-99; t2phi=-99; 
    j1pt=-99; j1eta=-99; j1phi=-99; 
    j2pt=-99; j2eta=-99; j2phi=-99; 
    j3pt=-99; j3eta=-99; j3phi=-99; 
    j4pt=-99; j4eta=-99; j4phi=-99; 
    mht=-99; ht=-99; 

    if (branchGenMET) {
      metObj=(MissingET*) branchGenMET->At(0);
      mht=metObj->MET;
    }
    else mht=3;

    if (branchHT) {
//.........这里部分代码省略.........

int main(int argc, char *argv[])
{
  char appName[] = "root2lhco";
  stringstream message;
  FILE *outputFile = 0;
  TChain *inputChain = 0;
  LHCOWriter *writer = 0;
  ExRootTreeReader *treeReader = 0;
  Long64_t entry, allEntries;

  if(argc < 2 || argc > 3)
  {
    cerr << " Usage: " << appName << " input_file" << " [output_file]" << endl;
    cerr << " input_file - input file in ROOT format," << endl;
    cerr << " output_file - output file in LHCO format," << endl;
    cerr << " with no output_file, or when output_file is -, write to standard output." << endl;
    return 1;
  }

  signal(SIGINT, SignalHandler);

  gROOT->SetBatch();

  int appargc = 1;
  char *appargv[] = {appName};
  TApplication app(appName, &appargc, appargv);

  try
  {
    cerr << "** Reading " << argv[1] << endl;
    inputChain = new TChain("Delphes");
    inputChain->Add(argv[1]);

    ExRootTreeReader *treeReader = new ExRootTreeReader(inputChain);

    if(argc == 2 || strcmp(argv[2], "-") == 0)
    {
      outputFile = stdout;
    }
    else
    {
      outputFile = fopen(argv[2], "w");

      if(outputFile == NULL)
      {
        message << "can't open " << argv[2];
        throw runtime_error(message.str());
      }
    }

    fprintf(outputFile, "   #  typ      eta      phi      pt    jmas   ntrk   btag  had/em   dum1   dum2\n");

    allEntries = treeReader->GetEntries();
    cerr << "** Input file contains " << allEntries << " events" << endl;

    if(allEntries > 0)
    {
      // Create LHC Olympics converter:
      writer = new LHCOWriter(treeReader, outputFile);

      ExRootProgressBar progressBar(allEntries - 1);
      // Loop over all events
      for(entry = 0; entry < allEntries && !interrupted; ++entry)
      {
        if(!treeReader->ReadEntry(entry))
        {
          cerr << "** ERROR: cannot read event " << entry << endl;
          break;
        }

        writer->ProcessEvent();

        progressBar.Update(entry);
      }
      progressBar.Finish();

      delete writer;
    }

    cerr << "** Exiting..." << endl;

    if(outputFile != stdout) fclose(outputFile);
    delete treeReader;
    delete inputChain;
    return 0;
  }
  catch(runtime_error &e)
  {
    if(writer) delete writer;
    if(treeReader) delete treeReader;
    if(inputChain) delete inputChain;
    cerr << "** ERROR: " << e.what() << endl;
    return 1;
  }
}

// Takes as input a given process with corresponding luminosity, as well as sample CME, cuts, directory location, and output file name - saves corresponding histogram in root file.
void processsimulation(TString sample, Double_t crosssection, TLorentzVector CME, Double_t luminosity, TString filename, TString directory, std::map<TString, Int_t> cut, std::map<TString, Int_t> particletype, std::clock_t start){
    
    cout << "Processing " << sample << endl;
    // Defines which sample to read from
    TString samplename = directory + sample + root;
    TChain chain("Delphes");
    chain.Add(samplename); 
    ExRootTreeReader *treeReader = new ExRootTreeReader(&chain); // Reads sample from delphes root fiile

    //Branches to be read from the collision file
    std::array<TString, 5> branchnames = {(TString)"Muon",(TString)"Electron",(TString)"Photon",(TString)"Jet", (TString)"Particle" };
                                       // (TString)"Tower",(TString)"Track",(TString)"EFlowTrack",(TString)"EFlowPhoton",(TString)"EFlowNeutralHadron"};
    std::map<TString, TClonesArray*> branches;
    for (int i=0; i<(int)branchnames.size(); i++){
        branches.insert ( std::pair<TString, TClonesArray*> (branchnames[i], treeReader->UseBranch(branchnames[i])));
    }

    // Variables to store particles
    Electron *electron = 0;
    Muon *muon = 0;
    Jet *jet = 0;

    Double_t numberOfEntries = treeReader->GetEntries();
    Double_t norm_const = (crosssection*luminosity)/numberOfEntries; // Normalizes samples based on cross section and luminosity

    TTree *ttree = new TTree(sample, sample);
    

    //Creates a branch for each variable listed in "parameters" in the main function
    std::map<TString, Float_t> *parameters = new std::map<TString, Float_t>;
    std::map<TString, Int_t>::iterator setparameters;
    for ( setparameters = cut.begin(); setparameters != cut.end(); setparameters++){
        parameters->insert( std::pair<TString, Float_t> (setparameters->first, 0));
    }

    std::map<TString, Float_t>::iterator setbranch;
    for (setbranch = parameters->begin(); setbranch != parameters->end(); setbranch++){
        TString type = "/F";
        TString vartype = setbranch->first + type;
        ttree->Branch(setbranch->first, &(*parameters)[setbranch->first],vartype);
    }

    //Timer
    double duration = 0;
    double durationnew = 0;
    for (Int_t iEntry = 0; iEntry < numberOfEntries; iEntry++){
        treeReader->ReadEntry(iEntry);
        TString skip = "";

        durationnew = (( std::clock() - start ) / (double) CLOCKS_PER_SEC)/60;
        if (durationnew - duration > 1){
            cout << "Still running - it has been "<< durationnew << " minutes" <<'\n';
            duration = durationnew;
        }

        if (particletype["muon"]==1){
            // Event selection process on Jets
            skip = parameter_functions(branches["Muon"], muon, "muon", branches, cut, CME, ttree, parameters);
            if (skip == "fill" || skip == "skip") continue;
        }
        if (particletype["electron"]==1){
            // Event selection process on Jets
            skip = parameter_functions(branches["Electron"], electron, "electron", branches, cut, CME, ttree, parameters);
            if (skip == "fill" || skip == "skip") continue;
        }
        // if (particletype["jet"]==1){
        //     // Event selection process on Jets
        //     skip = parameter_functions(branches["Jet"], jet, "jet", branches, cut, CME, ttree, parameters);
        //     if (skip == "fill" || skip == "skip") continue;
        // }
    }

    TFile f(filename, "update");
    ttree->Write();
    f.Close();

    ofstream myfile;
    myfile.open ("normalization.txt", ios::app);
        cout << "Normalization constant for " << sample << " is " << norm_const << endl;
        cout << "The sample " << sample << " has unweighted " << numberOfEntries << " entries" << endl;
        myfile << sample << " " << norm_const << "\n"; //write to file
        myfile << sample << " " << numberOfEntries << "\n";
    myfile.close();
}

void e6_rootMacro1vTree(const char *inputFile) {
    gSystem->Load("libDelphes");

    // Create chain of root trees
    //    TChain chain("Delphes");
    //    chain.Add(inputFile);
    TChain *chain = new TChain("Delphes"); // adopted from Example3.C of Delphes
    chain->Add(inputFile);

    // Create object of class ExRootTreeReader
    //ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
    ExRootTreeReader *treeReader = new ExRootTreeReader(chain); // adopted from Example3.C of Delphes
    Long64_t numberOfEntries = treeReader->GetEntries();


    string histoFile_str = "e6_rootMacro1vTree.root";
    // TFile constructor accepts type "const char*"
    const char* histoFile_char = histoFile_str.c_str();
    // overwrite existing ".root" file
    TFile f(histoFile_char, "recreate");

    // number of particles in an event
    Long64_t numParticles = -1;
    Long64_t numElectrons = -1;
    Long64_t numMuons = -1;
    Long64_t numJets = -1;

    Double_t maxPT = -1;
    Double_t max2ndPT = 0; // 2nd highest PT
    Double_t currentPT = 0; // initially currentPT>maxPT so that first particle in the loop will have max. PT

    Double_t currentMass = 0;

    // Get pointers to branches used in this analysis
    TClonesArray *branchParticle = treeReader->UseBranch("Particle");
    TClonesArray *branchElectron = treeReader->UseBranch("Electron");
    TClonesArray *branchMuon = treeReader->UseBranch("Muon");
    TClonesArray *branchJet = treeReader->UseBranch("Jet");

    // Book histograms
    //////////// I will use a tree to study electron with maximum PT.
    TTree t_electron_with_max_PT("electron_with_max_PT", "ein Baum über Elektronen mit höchsten PT");

    int electron_id = 11;
    Double_t ch_e, mass_e, E_e, px_e, py_e, pz_e, pt_e, eta_e, phi_e, rapidity_e, stat_e, pid_e;

    t_electron_with_max_PT.Branch("electron.charge", &ch_e, "ch_e/D");
    //    t_electron_with_max_PT.Branch("electron.mass", &mass_e, "mass_e/D");
    //    t_electron_with_max_PT.Branch("electron.energy", &E_e, "E_e/D");
    //    t_electron_with_max_PT.Branch("electron.Px", &px_e, "px_e/D");
    //    t_electron_with_max_PT.Branch("electron.Py", &py_e, "py_e/D");
    //    t_electron_with_max_PT.Branch("electron.Pz", &pz_e, "pz_e/D");
    t_electron_with_max_PT.Branch("electron.PT", &pt_e, "pt_e/D");
    t_electron_with_max_PT.Branch("electron.Eta", &eta_e, "eta_e/D");
    t_electron_with_max_PT.Branch("electron.Phi", &phi_e, "phi_e/D");
    //    t_electron_with_max_PT.Branch("electron.Rapidity", &rapidity_e, "rapidity_e/D");
    //    t_electron_with_max_PT.Branch("electron.status", &stat_e, "stat_e/D");
    //    t_electron_with_max_PT.Branch("electron.PID", &pid_e, "pid_e/D");

    //////////// I will use a tree to study muon with maximum PT.
    TTree t_muon_with_max_PT("muon_with_max_PT", "ein Baum über Muonen mit höchsten PT");

    int muon_id = 13;
    Double_t ch_muon, mass_muon, E_muon, px_muon, py_muon, pz_muon, pt_muon, eta_muon, phi_muon, rapidity_muon, stat_muon, pid_muon;

    t_muon_with_max_PT.Branch("muon.charge", &ch_muon, "ch_muon/D");
    //    t_muon_with_max_PT.Branch("muon.mass", &mass_muon, "mass_muon/D");
    //    t_muon_with_max_PT.Branch("muon.energy", &E_muon, "E_muon/D");
    //    t_muon_with_max_PT.Branch("muon.Px", &px_muon, "px_muon/D");
    //    t_muon_with_max_PT.Branch("muon.Py", &py_muon, "py_muon/D");
    //    t_muon_with_max_PT.Branch("muon.Pz", &pz_muon, "pz_muon/D");
    t_muon_with_max_PT.Branch("muon.PT", &pt_muon, "pt_muon/D");
    t_muon_with_max_PT.Branch("muon.Eta", &eta_muon, "eta_muon/D");
    t_muon_with_max_PT.Branch("muon.Phi", &phi_muon, "phi_muon/D");
    //    t_muon_with_max_PT.Branch("muon.Rapidity", &rapidity_muon, "rapidity_muon/D");
    //    t_muon_with_max_PT.Branch("muon.status", &stat_muon, "stat_muon/D");
    //    t_muon_with_max_PT.Branch("muon.PID", &pid_muon, "pid_muon/D");

    //////////// I will use a tree to study jet with maximum PT.
    TTree t_jet_with_max_PT("jet_with_max_PT", "ein Baum über jets mit höchsten PT");

    Double_t ch_jet, mass_jet, E_jet, px_jet, py_jet, pz_jet, pt_jet, eta_jet, phi_jet, rapidity_jet, stat_jet, pid_jet;

    t_jet_with_max_PT.Branch("jet.charge", &ch_jet, "ch_jet/D");
    t_jet_with_max_PT.Branch("jet.mass", &mass_jet, "mass_jet/D");
    //    t_jet_with_max_PT.Branch("jet.energy", &E_jet, "E_jet/D");
    //    t_jet_with_max_PT.Branch("jet.Px", &px_jet, "px_jet/D");
    //    t_jet_with_max_PT.Branch("jet.Py", &py_jet, "py_jet/D");
    //    t_jet_with_max_PT.Branch("jet.Pz", &pz_jet, "pz_jet/D");
    t_jet_with_max_PT.Branch("jet.PT", &pt_jet, "pt_jet/D");
    t_jet_with_max_PT.Branch("jet.Eta", &eta_jet, "eta_jet/D");
    t_jet_with_max_PT.Branch("jet.Phi", &phi_jet, "phi_jet/D");
    //    t_jet_with_max_PT.Branch("jet.Rapidity", &rapidity_jet, "rapidity_jet/D");
    //    t_jet_with_max_PT.Branch("jet.status", &stat_jet, "stat_jet/D");
    //    t_jet_with_max_PT.Branch("jet.PID", &pid_jet, "pid_jet/D");

    //////////// I will use a tree to study jet with 2nd maximum PT.
    TTree t_jet_with_2ndmax_PT("jet_with_2ndmax_PT", "ein Baum über jets mit zweiten höchsten PT");

    Double_t ch_jet2nd, mass_jet2nd, E_jet2nd, px_jet2nd, py_jet2nd, pz_jet2nd, pt_jet2nd, eta_jet2nd, phi_jet2nd, rapidity_jet2nd, stat_jet2nd, pid_jet2nd;
//.........这里部分代码省略.........

void fillHistogram(TH1F * hist, std::string filename, float xs){


  hist->Sumw2();

  //see https://answers.launchpad.net/mg5amcnlo/+question/200336
  int n_events_nonzero_met = 0;

    // Load shared library
  gSystem->Load("lib/libExRootAnalysis.so");
  gSystem->Load("libPhysics");

  // Create chain of root trees
  TChain chain("LHCO");
  //chain.Add("/scratch3/anlevin/MG5_aMC_v2_0_0/ww_qed_4_qcd_99_pythia_matching/Events/run_01/tag_1_pgs_events.root");
  chain.Add(filename.c_str());
  
  // Create object of class ExRootTreeReader
  ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
  Long64_t numberOfEntries = treeReader->GetEntries();
  
  // Get pointers to branches used in this analysis
  TClonesArray *branchJet = treeReader->UseBranch("Jet");
  TClonesArray *branchElectron = treeReader->UseBranch("Electron");
  TClonesArray *branchMuon = treeReader->UseBranch("Muon");
  TClonesArray *branchMET = treeReader->UseBranch("MissingET");

  
  // Loop over all events
  for(Int_t entry = 0; entry < numberOfEntries; ++entry) {

    // Load selected branches with data from specified event
    treeReader->ReadEntry(entry);
  
    if (branchMET->GetEntries() == 0)
      continue;
    else 
      assert (branchMET->GetEntries() == 1);

    n_events_nonzero_met++;

    // If event contains at least 1 jet
    if(branchJet->GetEntries() <= 2)
      continue;

    TRootJet *jet1 = 0; 
    TRootJet *jet2 = 0; 

    for (int i = 0; i < branchJet->GetEntries(); i++){
      TRootJet *jet = (TRootJet*) branchJet->At(i);
      if (jet1 == 0)
	jet1 = jet; 
      else if (jet2 == 0 )
	jet2 = jet;
      else if (jet->PT > jet1->PT)
	jet1 = jet;
      else if (jet->PT > jet2->PT)
	jet2 = jet;
    }
    
    if (!jet1 || !jet2 )
      continue;

    if (jet1->PT < 30)
      continue;

    if (jet2->PT < 30)
      continue;

    TLorentzVector vec1, vec2;
    bool found1 = false;
    bool found2 = false;
    int charge1 = 0;
    int charge2 = 0;

    for (int i = 0; i < branchMuon->GetEntries(); i++){
      TRootMuon * muon = (TRootMuon*) branchJet->At(i);
      TLorentzVector vec;
      vec.SetPtEtaPhiM(muon->PT, muon->Eta, muon->Phi, muon_mass);

      if (!found1){
	found1 = true;
	vec1 = vec;
	charge1=muon->Charge;
      }
      else if (!found2){
	found2 = true;
	vec2 = vec;
	charge2=muon->Charge;
      }
      else if (vec.Pt() > vec1.Pt()) {
	vec1 = vec;
	charge1=muon->Charge;
      }
      else if (vec.Pt() > vec2.Pt()) {
	vec2 = vec;
	charge2=muon->Charge;
      }

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

int main(int argc, char *argv[])
{
  char appName[] = "root2pileup";
  stringstream message;
  TChain *inputChain = 0;
  ExRootTreeReader *treeReader = 0;
  TClonesArray *branchParticle = 0;
  TIterator *itParticle = 0;
  GenParticle *particle = 0;
  DelphesPileUpWriter *writer = 0;
  Long64_t entry, allEntries;
  Int_t i;

  if(argc < 3)
  {
    cout << " Usage: " << appName << " output_file" << " input_file(s)" << endl;
    cout << " output_file - output binary pile-up file," << endl;
    cout << " input_file(s) - input file(s) in ROOT format." << endl;
    return 1;
  }

  signal(SIGINT, SignalHandler);

  gROOT->SetBatch();

  int appargc = 1;
  char *appargv[] = {appName};
  TApplication app(appName, &appargc, appargv);

  try
  {
    inputChain = new TChain("Delphes");
    for(i = 2; i < argc && !interrupted; ++i)
    {
      inputChain->Add(argv[i]);
    }

    treeReader = new ExRootTreeReader(inputChain);
    branchParticle = treeReader->UseBranch("Particle");
    itParticle = branchParticle->MakeIterator();

    writer = new DelphesPileUpWriter(argv[1]);

    allEntries = treeReader->GetEntries();
    cout << "** Input file(s) contain(s) " << allEntries << " events" << endl;

    if(allEntries > 0)
    {
      ExRootProgressBar progressBar(allEntries - 1);
      // Loop over all events in the input file
      for(entry = 0; entry < allEntries && !interrupted; ++entry)
      {
        if(!treeReader->ReadEntry(entry))
        {
          cerr << "** ERROR: cannot read event " << entry << endl;
          break;
        }

        itParticle->Reset();
        while((particle = static_cast<GenParticle*>(itParticle->Next())))
        {
          writer->WriteParticle(particle->PID,
            particle->X, particle->Y, particle->Z, particle->T,
            particle->Px, particle->Py, particle->Pz, particle->E);
        }
        
        writer->WriteEntry();

        progressBar.Update(entry);
      }
      progressBar.Finish();

      writer->WriteIndex();
    }

    cout << "** Exiting..." << endl;

    delete writer;
    delete itParticle;
    delete treeReader;
    delete inputChain;
    return 0;
  }
  catch(runtime_error &e)
  {
    if(writer) delete writer;
    if(itParticle) delete itParticle;
    if(treeReader) delete treeReader;
    if(inputChain) delete inputChain;
    cerr << "** ERROR: " << e.what() << endl;
    return 1;
  }
}

int
main (int argc, char *argv[])
{
  if (argc != 3)
    {
      cout << "Usage: " << argv[0] << " INPUT_FILE OUTPUT_FILE" << endl;
      cout << "  or:  " << argv[0] << " INPUT_LIST OUTPUT_FILE" << endl;
      cout << "Run vertexing performance analysis on INPUT_FILE and store output to OUTPUT_FILE. If the" << endl;
      cout << "first argument is a text file, analysis is run on each file listed inside." << endl;
      cout << endl;

      return 0;
    }
  string inputFile = argv[1], outputFile = argv[2], upperInputFile;
  upperInputFile.resize (inputFile.length ());
  transform (inputFile.begin (), inputFile.end (), upperInputFile.begin (), ::toupper);

  // so canvases don't appear on the screen when being created
  // very useful when running on the OSU T3 from CERN
  gROOT->SetBatch();

  //Create chain of root trees
  TChain chain("Delphes");

  if (upperInputFile.length () < 5 || upperInputFile.substr (upperInputFile.length () - 5, 5) != ".ROOT")
  {
    ifstream fin (inputFile);
    string line;
    while(getline(fin, line))
      {
	chain.Add(line.c_str());
      }
    fin.close();
  }
  else
    chain.Add(inputFile.c_str());

  
  // Create object of class ExRootTreeReader
  ExRootTreeReader *treeReader = new ExRootTreeReader(&chain);
  Long64_t numberOfEntries = treeReader->GetEntries();
  cout << "processing " << numberOfEntries << " events..." << endl << endl;

  // Get pointers to branches used in this analysis
  TClonesArray *branchTrack = treeReader->UseBranch("Track");
  TClonesArray *branchNPU = treeReader->UseBranch("NPU");
  TClonesArray *branchCluster = treeReader->UseBranch("Cluster");
  TClonesArray *branchVBFquarks = treeReader->UseBranch("VBFquarks");


  ///////////////////////////////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////////////////////////////
  // Book histograms
  ///////////////////////////////////////////////////////////////////////////////////////////////
  ///////////////////////////////////////////////////////////////////////////////////////////////

  TH1::SetDefaultSumw2();

  // VBF Quark Histograms
  TH1D *hVBFQuarkEta   = new TH1D("VBFQuarkEta", ";VBF quark #eta", 1000, -5.0, 5.0);
  TH1D *hVBFQuarkPt   = new TH1D("VBFQuarkPt", ";VBF quark p_{T} [GeV]", 1000, 0.0, 1000.0);
  TH2D *hVBFQuarkEtaVsPhi   = new TH2D("VBFQuarkPt", ";VBF quark #phi;VBF quark #eta", 500, -3.5, 3.5, 500, -5.0, 5.0);

  // Track Histograms
  TH1D *hTrackPt   = new TH1D("TrackPt", "; Track p_{T} [GeV]", 400, 0.0, 200.0);
  TH1D *hVBFTrackPt   = new TH1D("VBFTrackPt", ";VBF track p_{T} [GeV]", 400, 0.0, 200.0);
  TH1D *hPUTrackPt   = new TH1D("PUTrackPt", ";PU track p_{T} [GeV]", 400, 0.0, 200.0);
  TH1D *hTrackX = new TH1D("TrackX", "; Track x [mm]", 300, -3, 3);
  TH1D *hTrackY = new TH1D("TrackY", "; Track y [mm]", 300, -3, 3);
  TH1D *hTrackZ = new TH1D("TrackZ", "; Track z [mm]", 300, -300, 300);

  TH1D *hTrackEX = new TH1D("TrackEX", "; Track #sigma_{x} [mm]", 100, 0, 1);
  TH1D *hTrackEY = new TH1D("TrackEY", "; Track #sigma_{y} [mm]", 100, 0, 1);
  TH1D *hTrackEZ = new TH1D("TrackEZ", "; Track #sigma_{z} [mm]", 100, 0, 20);

  TH1D *hTrackXSmearingShift = new TH1D("TrackXSmearingShift", "; Track x_{smeared} - x_{true} [mm]", 100, -2, 2);
  TH1D *hTrackYSmearingShift = new TH1D("TrackYSmearingShift", "; Track y_{smeared} - y_{true} [mm]", 100, -2, 2);
  TH1D *hTrackZSmearingShift = new TH1D("TrackZSmearingShift", "; Track z_{smeared} - z_{true} [mm]", 100, -50, 50);

  TH1D *hTrackDisplacementClusterX = new TH1D("TrackDisplacementClusterX", "; |Track x - Cluster x| [mm]", 100, 0, 10);
  TH1D *hTrackDisplacementClusterY = new TH1D("TrackDisplacementClusterY", "; |Track y - Cluster y| [mm]", 100, 0, 10);
  TH1D *hTrackDisplacementClusterZ = new TH1D("TrackDisplacementClusterZ", "; |Track z - Cluster z| [mm]", 100, 0, 300);

  TH1D *hTrackDisplacementInteractionX = new TH1D("TrackDisplacementInteractionX", "; |Track x - Interaction x| [mm]", 100, 0, 10);
  TH1D *hTrackDisplacementInteractionY = new TH1D("TrackDisplacementInteractionY", "; |Track y - Interaction y| [mm]", 100, 0, 10);
  TH1D *hTrackDisplacementInteractionZ = new TH1D("TrackDisplacementInteractionZ", "; |Track z - Interaction z| [mm]", 100, 0, 50);

  // Cluster Histograms
  TH1D *hClusterX = new TH1D("ClusterX", "; Cluster x [mm]", 300, -3, 3);
  TH1D *hClusterEX = new TH1D("ClusterEX", "; Cluster #sigmax [mm]", 100, 0, 0.1);
  TH1D *hClusterY = new TH1D("ClusterY", "; Cluster y [mm]", 300, -3, 3);
  TH1D *hClusterEY = new TH1D("ClusterEY", "; Cluster #sigmay [mm]", 100, 0, 0.1);
  TH1D *hClusterZ = new TH1D("ClusterZ", "; Cluster z [mm]", 300, -300, 300);
  TH1D *hClusterEZ = new TH1D("ClusterEZ", "; Cluster #sigmaz [mm]", 100, 0, 0.1);
  TH2D *hClusterTransversePosition = new TH2D("ClusterTransversePosition", "; Cluster x [mm]; Cluster y [mm]", 300, -3,3, 300, -3, 3);
  TH2D *hClusterRZPosition = new TH2D("ClusterRZPosition", "; Cluster z [mm]; Cluster r [mm]", 100, -300, 300, 100, 0, 1.5);

  // Reco PV Histograms
  TH1D *hRecoPVX = new TH1D("RecoPVX", "; Reco PV x [mm]", 300, -3, 3);
  TH1D *hRecoPVY = new TH1D("RecoPVY", "; Reco PV y [mm]", 300, -3, 3);
//.........这里部分代码省略.........