void Rivet_Interface::SetEventWeight(const Rivet_Scale_Variation* rsv,
                                     HepMC::GenEvent& evt, const int& idx)
{
  double wgt(idx<0?rsv->Weight():rsv->Weight(idx));
  DEBUG_FUNC(rsv->Name()<<": "<<wgt);
  evt.weights()[0]=wgt;
#ifdef HEPMC_HAS_CROSS_SECTION
  evt.cross_section()->set_cross_section(rsv->TotalXS(),rsv->TotalErr());
#endif
}

int
DirectPhotonPythia::process_event(PHCompositeNode *topNode)
{
  _ievent ++;

  PHHepMCGenEventMap * geneventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
  if (!geneventmap)
  {
    std::cout <<PHWHERE<<" - Fatal error - missing node PHHepMCGenEventMap"<<std::endl;
    return Fun4AllReturnCodes::ABORTRUN;
  }

  PHHepMCGenEvent *genevt = geneventmap->get(_embedding_id);
  if (!genevt)
  {
    std::cout <<PHWHERE<<" - Fatal error - node PHHepMCGenEventMap missing subevent with embedding ID "<<_embedding_id;
    std::cout <<". Print PHHepMCGenEventMap:";
    geneventmap->identify();
    return Fun4AllReturnCodes::ABORTRUN;
  }

  HepMC::GenEvent* theEvent = genevt->getEvent();
  

for (HepMC::GenEvent::particle_const_iterator p = theEvent->particles_begin();
      p != theEvent->particles_end(); ++p)
    {
      TParticlePDG * pdg_p = TDatabasePDG::Instance()->GetParticle(
          (*p)->pdg_id());
      if (pdg_p)
        {
          if (TString(pdg_p->GetName()) == "gamma")
            {
              // cout <<"Find gamma at eta = "<< (*p)->momentum().pseudoRapidity() <<endl;

	      if((*p)->momentum().perp()>5.) _h2->Fill((*p)->momentum().perp(), (*p)->momentum().pseudoRapidity());
              if((*p)->momentum().perp()>5.) _h1->Fill((*p)->momentum().perp());
              
              if((*p)->momentum().perp()>5.) {
		const float px = (*p)->momentum().x(), 
                            py = (*p)->momentum().y(), 
                            pz = (*p)->momentum().z(),
		            pt = (*p)->momentum().perp(), 
                           eta = (*p)->momentum().pseudoRapidity(),
                           phi = (*p)->momentum().phi();
                float shower_data[70] = {(float)_ievent,px,py,pz,pt,eta,phi};
               
                _ntp_gamma -> Fill(shower_data);
             }
            }
        }
      }

      return 0;;
}

void Output_Delphes::AnalyseEvent
(ExRootTreeBranch *branch,const HepMC::GenEvent& evt,
 const Long64_t eventNumber,const double weight)
{
  TRootLHEFEvent *element;
  element = static_cast<TRootLHEFEvent*>(branch->NewEntry());
  element->Number = eventNumber;
  element->Nparticles = evt.particles_size();
  element->ProcessID = evt.signal_process_id();
  element->Weight = weight;
  element->ScalePDF = evt.event_scale();
  element->CouplingQED = evt.alphaQED();
  element->CouplingQCD = evt.alphaQCD();
}

/// Smearing function
StatusCode GaussSmearVertex::smearVertex(HepMC::GenEvent& theEvent) {

  double dx = m_gaussDist() * sqrt(m_xsig) + m_xmean;
  double dy = m_gaussDist() * sqrt(m_ysig) + m_ymean;
  double dz = m_gaussDist() * sqrt(m_zsig) + m_zmean;
  double dt = m_gaussDist() * sqrt(m_tsig) + m_tmean;

  Gaudi::LorentzVector dpos(dx, dy, dz, dt);

  debug() << "Smearing vertices by " << dpos << endmsg;

  for (auto vit = theEvent.vertices_begin(); vit != theEvent.vertices_end(); ++vit) {
    Gaudi::LorentzVector pos((*vit)->position());
    pos += dpos;
    (*vit)->set_position(HepMC::FourVector(pos.x(), pos.y(), pos.z(), pos.t()));
  }

  return StatusCode::SUCCESS;
}

int main(int argc, char *argv[]){

  TFile* fout = new TFile("cascade.result.root","RECREATE");

  HepMC::IO_GenEvent fin("../CASCADE_HEPMC/bin/example_test_out.dat",std::ios::in);
  
  HepMC::GenEvent* evt = new HepMC::GenEvent;

  TTree *tree = new TTree("tree","cascade result");
  tree->Branch("event",&evt);

  Int_t nevt = 0;

  while (evt = fin.read_next_event() ) {

    nevt++;

    if(debug){ 
      cout<<"nevt: "<<nevt<<endl;
      getchar();
      evt->print();
    }

    tree->Fill();
  }

  tree->Print();

  fout->cd();
  tree->Write();
  fout->Close();
  
  delete fout;

  return(0);  
}

StatusCode ConstPtParticleGun::getNextEvent(HepMC::GenEvent& theEvent) {
  Gaudi::LorentzVector theFourMomentum;
  Gaudi::LorentzVector origin;
  // note: pgdid is set in function generateParticle
  int thePdgId;
  generateParticle(theFourMomentum, origin, thePdgId);

  // create HepMC Vertex --
  // by calling add_vertex(), the hepmc event is given ownership of the vertex
  HepMC::GenVertex* v = new HepMC::GenVertex(HepMC::FourVector(origin.X(), origin.Y(), origin.Z(), origin.T()));
  // create HepMC particle --
  // by calling add_particle_out(), the hepmc vertex is given ownership of the particle
  HepMC::GenParticle* p = new HepMC::GenParticle(
      HepMC::FourVector(theFourMomentum.Px(), theFourMomentum.Py(), theFourMomentum.Pz(), theFourMomentum.E()),
      thePdgId,
      1);  // hepmc status code for final state particle

  v->add_particle_out(p);

  theEvent.add_vertex(v);
  theEvent.set_signal_process_vertex(v);

  return StatusCode::SUCCESS;
}

StatusCode EDMToHepMCConverter::execute() {

  const fcc::MCParticleCollection* particles = m_genphandle.get();
  // ownership of event given to data service at the end of execute
  HepMC::GenEvent* event = new HepMC::GenEvent;

  // conversion of units to EDM standard units:
  // First cover the case that hepMC file is not in expected units and then convert to EDM default
  double hepmc2EdmLength = conversion_factor(event->length_unit(), gen::hepmcdefault::length) * gen::hepmc2edm::length;
  double hepmc2EdmEnergy =
      conversion_factor(event->momentum_unit(), gen::hepmcdefault::energy) * gen::hepmc2edm::energy;

  for (auto p : *(particles)) {
    if (p.status() == 1) {  // only final state particles
      GenParticle* pHepMC =
          new GenParticle(HepMC::FourVector(p.p4().px, p.p4().py, p.p4().pz, p.p4().mass / hepmc2EdmEnergy),
                          p.pdgId(),
                          p.status());  // hepmc status code for final state particle

      fcc::ConstGenVertex vStart = p.startVertex();
      if (p.startVertex().isAvailable()) {
        HepMC::GenVertex* v =
            new HepMC::GenVertex(HepMC::FourVector(vStart.position().x / hepmc2EdmLength,
                                                   vStart.position().y / hepmc2EdmLength,
                                                   vStart.position().z / hepmc2EdmLength,
                                                   vStart.ctau() / Gaudi::Units::c_light / hepmc2EdmLength));

        v->add_particle_out(pHepMC);
        event->add_vertex(v);
      }
    }
  }

  m_hepmchandle.put(event);
  return StatusCode::SUCCESS;
}

void Output_Delphes::AnalyseParticles(ExRootTreeBranch *branch, const HepMC::GenEvent& evt)
{
  TRootC::GenParticle *element;
  TLorentzVector momentum;
  Double_t signPz;
  ReadStats();
  for(int n=1; n<=evt.particles_size(); n++) {
    int mo1, mo2, da1, da2, status, pid;
    getStatsFromTuple(mo1,mo2,da1,da2,status,pid,n);
    element = static_cast<TRootC::GenParticle*>(branch->NewEntry());
    element->PID = pid;
    element->Status = status;
    element->M1 = mo1 - 1; // added -1 as the numbering in the tree starts from 0
    element->M2 = mo2 - 1;
    element->D1 = da1 - 1;
    element->D2 = da2 - 1;
    element->E = index_to_particle[n]->momentum().e();
    element->Px = index_to_particle[n]->momentum().px();
    element->Py = index_to_particle[n]->momentum().py();
    element->Pz = index_to_particle[n]->momentum().pz();
    element->PT = sqrt(pow(element->Px,2)+pow(element->Py,2));
    momentum.SetPxPyPzE(element->Px, element->Py, element->Pz, element->E);
    signPz = (element->Pz >= 0.0) ? 1.0 : -1.0;
    element->Eta = element->PT < 1e-6 ? signPz*999.9 : momentum.Eta();
    element->Phi = index_to_particle[n]->momentum().phi();
    HepMC::GenVertex* vrtI = (index_to_particle[n])->production_vertex();
    HepMC::GenVertex::particles_in_const_iterator partI;
    if(vrtI) {
      element->T = vrtI->position().t();
      element->X = vrtI->position().x();
      element->Y = vrtI->position().y();
      element->Z = vrtI->position().z();
    }
    else {
      element->T = 0.;
      element->X = 0.;
      element->Y = 0.;
      element->Z = 0.;
    }  
  }
}

//.........这里部分代码省略.........
    if (verbosity > 2) {
      cout << "PHSartre::process_event - triggersize: " << _registeredTriggers.size() << endl;
    }

    for (unsigned int tr = 0; tr < _registeredTriggers.size(); tr++) { 
      bool trigResult = _registeredTriggers[tr]->Apply(event);

      if (verbosity > 2) {
	cout << "PHSartre::process_event trigger: "
	     << _registeredTriggers[tr]->GetName() << "  " << trigResult << endl;
      }

      if (_triggersOR && trigResult) {
	passedTrigger = true;
	break;
      } else if (_triggersAND) {
	andScoreKeeper &= trigResult;
      }
      
      if (verbosity > 2 && !passedTrigger) {
	cout << "PHSartre::process_event - failed trigger: "
	     << _registeredTriggers[tr]->GetName() <<  endl;
      }
    }

    if ((andScoreKeeper && _triggersAND) || (_registeredTriggers.size() == 0)) {
      passedTrigger = true;
    }

  }

  // fill HepMC object with event
  
  HepMC::GenEvent *genevent = new HepMC::GenEvent(HepMC::Units::GEV, HepMC::Units::MM);

  // add some information to the event
  genevent->set_event_number(_eventcount);

  // Set the PDF information
  HepMC::PdfInfo pdfinfo;
  pdfinfo.set_scalePDF(event->Q2);
  genevent->set_pdf_info(pdfinfo); 

  // We would also like to save:
  //
  // event->t;
  // event->x;
  // event->y;
  // event->s;
  // event->W;
  // event->xpom;
  // (event->polarization == transverse ? 0 : 1);
  // (event->diffractiveMode == coherent ? 0 : 1);
  // 
  // but there doesn't seem to be a good place to do so 
  // within the HepMC event information?
  //
  // t, W and Q^2 form a minial set of good variables for diffractive events
  // Maybe what I do is record the input particles to the event at the HepMC
  // vertices and reconstruct the kinematics from there? 
  
  // Create HepMC vertices and add final state particles to them

  // First, the emitter(electron)-virtual photon vertex:

  HepMC::GenVertex* egammavtx = new HepMC::GenVertex(CLHEP::HepLorentzVector(0.0,0.0,0.0,0.0));

int main(int argc, char* argv[])
{
  TApplication theApp(srcName.Data(), &argc, argv);
//=============================================================================

  if (argc<5) return -1;
  TString sPath = argv[1]; if (sPath.IsNull()) return -1;
  TString sFile = argv[2]; if (sFile.IsNull()) return -1;
  TString sJetR = argv[3]; if (sJetR.IsNull()) return -1;
  TString sSjeR = argv[4]; if (sSjeR.IsNull()) return -1;
//=============================================================================

  sPath.ReplaceAll("#", "/");
//=============================================================================

  double dJetR = -1.;
  if (sJetR=="JetR02") dJetR = 0.2;
  if (sJetR=="JetR03") dJetR = 0.3;
  if (sJetR=="JetR04") dJetR = 0.4;
  if (sJetR=="JetR05") dJetR = 0.5;

  if (dJetR<0.) return -1;
  cout << "Jet R = " << dJetR << endl;
//=============================================================================

  double dSjeR = -1.;
  if (sSjeR=="SjeR01") dSjeR = 0.1;
  if (sSjeR=="SjeR02") dSjeR = 0.2;
  if (sSjeR=="SjeR03") dSjeR = 0.3;
  if (sSjeR=="SjeR04") dSjeR = 0.4;

  if (dSjeR<0.) return -1;
  cout << "Sub-jet R = " << dSjeR << endl;
//=============================================================================

  const double dJetsPtMin  = 0.001;
  const double dCutEtaMax  = 1.6;
  const double dJetEtaMax  = 1.;
  const double dJetAreaRef = TMath::Pi() * dJetR * dJetR;

  fastjet::GhostedAreaSpec areaSpc(dCutEtaMax);
  fastjet::JetDefinition   jetsDef(fastjet::antikt_algorithm, dJetR, fastjet::BIpt_scheme, fastjet::Best);
  fastjet::AreaDefinition  areaDef(fastjet::active_area_explicit_ghosts,areaSpc);

  fastjet::Selector selectJet = fastjet::SelectorAbsEtaMax(dJetEtaMax);
  fastjet::JetDefinition subjDef(fastjet::kt_algorithm, dSjeR, fastjet::BIpt_scheme, fastjet::Best);
//=============================================================================

  std::vector<fastjet::PseudoJet> fjInput;
  const Double_t dCut = TMath::TwoPi() / 3.;
//=============================================================================

  enum { kWgt, kXsc, kLje, kLjr, kLtk, kLtr, kJet, kAje, kMje, k1sz, k1sA, k1sm, k1sr, k2sz, k2sA, k2sm, k2sr, kDsm, kDsr, kVar };

  TFile *file = TFile::Open(Form("%s.root",sFile.Data()), "NEW");
  TNtuple *nt = new TNtuple("nt", "", "fWgt:fXsc:fLje:fLjr:fLtk:fLtr:fJet:fAje:fMje:f1sj:f1sA:f1sm:f1sr:f2sj:f2sA:f2sm:f2sr:fDsm:fDsr");
//=============================================================================

  HepMC::IO_GenEvent ascii_in(Form("%s/%s.hepmc",sPath.Data(),sFile.Data()), std::ios::in);
  HepMC::GenEvent *evt = ascii_in.read_next_event();

  while (evt) {
    fjInput.resize(0);

    double dLtk = -1.;
    TVector3 vPar, vLtk;
    for (HepMC::GenEvent::particle_const_iterator p=evt->particles_begin(); p!=evt->particles_end(); ++p) if ((*p)->status()==1) {
      double dEta = (*p)->momentum().eta(); if (TMath::Abs(dEta)>dCutEtaMax) continue;

      double dTrk = (*p)->momentum().perp();
      double dPhi = (*p)->momentum().phi();
      vPar.SetPtEtaPhi(dTrk, dEta, dPhi);
      fjInput.push_back(fastjet::PseudoJet(vPar.Px(), vPar.Py(), vPar.Pz(), vPar.Mag()));
      if (dTrk>dLtk) { dLtk = dTrk; vLtk.SetPtEtaPhi(dTrk, dEta, dPhi); }
    }
//=============================================================================

    fastjet::ClusterSequenceArea clustSeq(fjInput, jetsDef, areaDef);
    std::vector<fastjet::PseudoJet> includJets = clustSeq.inclusive_jets(dJetsPtMin);
    std::vector<fastjet::PseudoJet> selectJets = selectJet(includJets);
//=============================================================================

    if (selectJets.size()>0) {
      std::vector<fastjet::PseudoJet> sortedJets = fastjet::sorted_by_pt(selectJets);
      TVector3 vLje; vLje.SetPtEtaPhi(sortedJets[0].pt(), sortedJets[0].eta(), sortedJets[0].phi());

      TVector3 vJet;
      int kJrl = -1; double dJrl = -1.;
      int kTrl = -1; double dTrl = -1.;
      for (int j=0; j<sortedJets.size(); j++) {
        double dJet = sortedJets[j].pt();
        sortedJets[j].set_user_index(-1);
        vJet.SetPtEtaPhi(dJet, sortedJets[j].eta(), sortedJets[j].phi());
        if (TMath::Abs(vJet.DeltaPhi(vLje))>dCut) { sortedJets[j].set_user_index(1); if (dJet>dJrl) { dJrl = dJet; kJrl = j; } }
        if (TMath::Abs(vJet.DeltaPhi(vLtk))>dCut) { sortedJets[j].set_user_index(2); if (dJet>dTrl) { dTrl = dJet; kTrl = j; } }
      }
//=============================================================================

      TVector3 v1sj, v2sj;
      for (int j=0; j<sortedJets.size(); j++) {
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
  TApplication theApp(srcName.Data(), &argc, argv);
//=============================================================================

  if (argc<5) return -1;
  TString sPath = argv[1]; if (sPath.IsNull()) return -1;
  TString sFile = argv[2]; if (sFile.IsNull()) return -1;
  TString sJetR = argv[3]; if (sJetR.IsNull()) return -1;
  TString sSjeR = argv[4]; if (sSjeR.IsNull()) return -1;
//=============================================================================

  sPath.ReplaceAll("#", "/");
//=============================================================================

  double dJetR = -1.;
  if (sJetR=="JetR02") dJetR = 0.2;
  if (sJetR=="JetR03") dJetR = 0.3;
  if (sJetR=="JetR04") dJetR = 0.4;
  if (sJetR=="JetR05") dJetR = 0.5;

  if (dJetR<0.) return -1;
  cout << "Jet R = " << dJetR << endl;
//=============================================================================

  double dSjeR = -1.;
  if (sSjeR=="SjeR01") dSjeR = 0.1;
  if (sSjeR=="SjeR02") dSjeR = 0.2;
  if (sSjeR=="SjeR03") dSjeR = 0.3;
  if (sSjeR=="SjeR04") dSjeR = 0.4;

  if (dSjeR<0.) return -1;
  cout << "Sub-jet R = " << dSjeR << endl;
//=============================================================================

  const double dJetsPtMin  = 0.001;
  const double dCutEtaMax  = 1.6;
  const double dJetEtaMax  = 1.;
  const double dJetAreaRef = TMath::Pi() * dJetR * dJetR;

  fastjet::GhostedAreaSpec areaSpc(dCutEtaMax);
  fastjet::JetDefinition   jetsDef(fastjet::antikt_algorithm, dJetR, fastjet::E_scheme, fastjet::Best);

//fastjet::AreaDefinition  areaDef(fastjet::active_area,areaSpc);
  fastjet::AreaDefinition  areaDef(fastjet::active_area_explicit_ghosts,areaSpc);

//fastjet::JetDefinition   bkgsDef(fastjet::kt_algorithm, 0.2, fastjet::BIpt_scheme, fastjet::Best);
//fastjet::AreaDefinition  aBkgDef(fastjet::active_area_explicit_ghosts, areaSpc);

  fastjet::Selector selectJet = fastjet::SelectorAbsEtaMax(dJetEtaMax);
//fastjet::Selector selectRho = fastjet::SelectorAbsEtaMax(dCutEtaMax-0.2);
//fastjet::Selector selecHard = fastjet::SelectorNHardest(2);
//fastjet::Selector selectBkg = selectRho * (!(selecHard));
//fastjet::JetMedianBackgroundEstimator bkgsEstimator(selectBkg, bkgsDef, aBkgDef);
//fastjet::Subtractor                   bkgSubtractor(&bkgsEstimator);

  fastjet::JetDefinition subjDef(fastjet::kt_algorithm, dSjeR, fastjet::E_scheme, fastjet::Best);
//=============================================================================

  std::vector<fastjet::PseudoJet> fjInput;
  const double dMass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
//=============================================================================

  TList *list = new TList();
  TH1D *hJet = new TH1D("hJet", "", 1000, 0., 1000.); hJet->Sumw2(); list->Add(hJet);

  enum                    {  kJet,   kMje,   kDsz,   kIsm,   kZsm,   kDsr, kVar };
  const TString sHist[] = { "aJet", "aMje", "aDsz", "aIsm", "aZsm", "aDsr" };
  const Int_t      nv[] = {   1000,    150,  120,      150,    150,    500 };
  const Double_t dMin[] = {     0.,     0.,   0.,       0.,     0.,     0. };
  const Double_t dMax[] = {  1000.,   150.,  1.2,     150.,    1.5,     5. };

  THnSparseD *hs = new THnSparseD("hs", "", kVar, nv, dMin, dMax); hs->Sumw2();
  for (Int_t i=0; i<kVar; i++) hs->GetAxis(i)->SetName(sHist[i].Data()); list->Add(hs);
//=============================================================================

  HepMC::IO_GenEvent ascii_in(Form("%s/%s.hepmc",sPath.Data(),sFile.Data()), std::ios::in);
  HepMC::GenEvent *evt = ascii_in.read_next_event();

  while (evt) {
    fjInput.resize(0);

    TLorentzVector vPar;
    for (HepMC::GenEvent::particle_const_iterator p=evt->particles_begin(); p!=evt->particles_end(); ++p) if ((*p)->status()==1) {
      vPar.SetPtEtaPhiM((*p)->momentum().perp(), (*p)->momentum().eta(), (*p)->momentum().phi(), dMass);

      if ((TMath::Abs(vPar.Eta())<dCutEtaMax)) {
        fjInput.push_back(fastjet::PseudoJet(vPar.Px(), vPar.Py(), vPar.Pz(), vPar.E()));
      }
    }
//=============================================================================

    fastjet::ClusterSequenceArea clustSeq(fjInput, jetsDef, areaDef);
    std::vector<fastjet::PseudoJet> includJets = clustSeq.inclusive_jets(dJetsPtMin);
//  std::vector<fastjet::PseudoJet> subtedJets = bkgSubtractor(includJets);
    std::vector<fastjet::PseudoJet> selectJets = selectJet(includJets);
//  std::vector<fastjet::PseudoJet> sortedJets = fastjet::sorted_by_pt(selectJets);

    for (int j=0; j<selectJets.size(); j++) {
      double dJet = selectJets[j].pt();
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
  TApplication theApp(srcName.Data(), &argc, argv);
//=============================================================================

  if (argc<5) return -1;
  TString sPath = argv[1]; if (sPath.IsNull()) return -1;
  TString sFile = argv[2]; if (sFile.IsNull()) return -1;
  TString sJetR = argv[3]; if (sJetR.IsNull()) return -1;
  TString sSjeR = argv[4]; if (sSjeR.IsNull()) return -1;
//=============================================================================

  sPath.ReplaceAll("#", "/");
//=============================================================================

  double dJetR = -1.;
  if (sJetR=="JetR02") dJetR = 0.2;
  if (sJetR=="JetR03") dJetR = 0.3;
  if (sJetR=="JetR04") dJetR = 0.4;
  if (sJetR=="JetR05") dJetR = 0.5;

  if (dJetR<0.) return -1;
  cout << "Jet R = " << dJetR << endl;
//=============================================================================

  double dSjeR = -1.;
  if (sSjeR=="SjeR01") dSjeR = 0.1;
  if (sSjeR=="SjeR02") dSjeR = 0.2;
  if (sSjeR=="SjeR03") dSjeR = 0.3;
  if (sSjeR=="SjeR04") dSjeR = 0.4;

  if (dSjeR<0.) return -1;
  cout << "Sub-jet R = " << dSjeR << endl;
//=============================================================================

  const double dJetsPtMin  = 0.1;
  const double dJetEtaMax  = 2.;
  const double dCutEtaMax  = 2.6;

  fastjet::GhostedAreaSpec areaSpc(dCutEtaMax);
  fastjet::JetDefinition   jetsDef(fastjet::antikt_algorithm, dJetR, fastjet::BIpt_scheme, fastjet::Best);
  fastjet::AreaDefinition  areaDef(fastjet::active_area_explicit_ghosts,areaSpc);

  fastjet::Selector selectJet = fastjet::SelectorAbsEtaMax(dJetEtaMax);
  fastjet::JetDefinition subjDef(fastjet::kt_algorithm, dSjeR, fastjet::BIpt_scheme, fastjet::Best);
//=============================================================================

  std::vector<fastjet::PseudoJet> fjInput;
//=============================================================================

  TFile *file = TFile::Open(Form("%s.root",sFile.Data()), "NEW");
  TList *list = new TList();

  TH1D *hWeightSum = new TH1D("hWeightSum", "", 1, 0., 1.); list->Add(hWeightSum);

  TH2D *hTrkPtEta = new TH2D("hTrkPtEta", "", 1000, 0., 500., 60, -3., 3.); hTrkPtEta->Sumw2(); list->Add(hTrkPtEta);
  TH2D *hTrkPtPhi = new TH2D("hTrkPtPhi", "", 1000, 0., 500., 20, -1., 1.); hTrkPtPhi->Sumw2(); list->Add(hTrkPtPhi);

  TH2D *hJetPtNsj = new TH2D("hJetPtNsj", "", 1000, 0., 1000., 100, -0.5, 99.5); hJetPtNsj->Sumw2(); list->Add(hJetPtNsj);


  TH2D *hJetPtEta = new TH2D("hJetPtEta", "", 1000, 0., 1000., 60, -3., 3.); hJetPtEta->Sumw2(); list->Add(hJetPtEta);
  TH2D *hJetPtPhi = new TH2D("hJetPtPhi", "", 1000, 0., 1000., 20, -1., 1.); hJetPtPhi->Sumw2(); list->Add(hJetPtPhi);

  TH2D *hLjePtEta = new TH2D("hLjePtEta", "", 1000, 0., 1000., 60, -3., 3.); hLjePtEta->Sumw2(); list->Add(hLjePtEta);
  TH2D *hLjePtPhi = new TH2D("hLjePtPhi", "", 1000, 0., 1000., 20, -1., 1.); hLjePtPhi->Sumw2(); list->Add(hLjePtPhi);

  TH2D *hNjePtEta = new TH2D("hNjePtEta", "", 1000, 0., 1000., 60, -3., 3.); hNjePtEta->Sumw2(); list->Add(hNjePtEta);
  TH2D *hNjePtPhi = new TH2D("hNjePtPhi", "", 1000, 0., 1000., 20, -1., 1.); hNjePtPhi->Sumw2(); list->Add(hNjePtPhi);


  TH2D *hJe2PtEta = new TH2D("hJe2PtEta", "", 1000, 0., 1000., 60, -3., 3.); hJe2PtEta->Sumw2(); list->Add(hJe2PtEta);
  TH2D *hJe2PtPhi = new TH2D("hJe2PtPhi", "", 1000, 0., 1000., 20, -1., 1.); hJe2PtPhi->Sumw2(); list->Add(hJe2PtPhi);

  TH2D *hLj2PtEta = new TH2D("hLj2PtEta", "", 1000, 0., 1000., 60, -3., 3.); hLj2PtEta->Sumw2(); list->Add(hLj2PtEta);
  TH2D *hLj2PtPhi = new TH2D("hLj2PtPhi", "", 1000, 0., 1000., 20, -1., 1.); hLj2PtPhi->Sumw2(); list->Add(hLj2PtPhi);

  TH2D *hNj2PtEta = new TH2D("hNj2PtEta", "", 1000, 0., 1000., 60, -3., 3.); hNj2PtEta->Sumw2(); list->Add(hNj2PtEta);
  TH2D *hNj2PtPhi = new TH2D("hNj2PtPhi", "", 1000, 0., 1000., 20, -1., 1.); hNj2PtPhi->Sumw2(); list->Add(hNj2PtPhi);
//=============================================================================

  HepMC::IO_GenEvent ascii_in(Form("%s/%s.hepmc",sPath.Data(),sFile.Data()), std::ios::in);
  HepMC::GenEvent *evt = ascii_in.read_next_event();

  while (evt) {
    fjInput.resize(0);
    double dWeight = evt->weights().back();
    double dXsect  = evt->cross_section()->cross_section() / 1e9;
    double dNorm   = dWeight * dXsect;
    hWeightSum->Fill(0.5, dWeight);

    TVector3 vTrk;
    for (HepMC::GenEvent::particle_const_iterator p=evt->particles_begin(); p!=evt->particles_end(); ++p) if ((*p)->status()==1) {
      double dTrkPt  = (*p)->momentum().perp(); if (dTrkPt<0.5) continue;
      double dTrkEta = (*p)->momentum().eta();  if (TMath::Abs(dTrkEta)>dCutEtaMax) continue;
      double dTrkPhi = (*p)->momentum().phi();  vTrk.SetPtEtaPhi(dTrkPt, dTrkEta, dTrkPhi);
      double dTrkCos = TMath::Cos(dTrkPhi); if (dTrkCos==1.) dTrkCos = 1. - 1e-6;

      fjInput.push_back(fastjet::PseudoJet(vTrk.Px(), vTrk.Py(), vTrk.Pz(), vTrk.Mag()));

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

bool HepMC2_Interface::Sherpa2HepMC(ATOOLS::Blob_List *const blobs,
                                    HepMC::GenEvent& event, double weight)
{
#ifdef USING__HEPMC2__UNITS
  event.use_units(HepMC::Units::GEV,
                  HepMC::Units::MM);
#endif
  event.set_event_number(ATOOLS::rpa->gen.NumberOfGeneratedEvents());
  size_t decid(11);
  std::map<size_t,size_t> decids;
  Blob *sp(blobs->FindFirst(btp::Signal_Process));
  if (sp) {
    Blob_Data_Base *info((*sp)["Decay_Info"]);
    if (info) {
      DecayInfo_Vector decs(info->Get<DecayInfo_Vector>());
      for (size_t i(0);i<decs.size();++i) decids[decs[i]->m_id]=++decid;
    }
  }
  m_blob2genvertex.clear();
  m_particle2genparticle.clear();
  HepMC::GenVertex * vertex;
  std::vector<HepMC::GenParticle*> beamparticles;
  for (ATOOLS::Blob_List::iterator blit=blobs->begin();
       blit!=blobs->end();++blit) {
    if (Sherpa2HepMC(*(blit),vertex,decids)) {
      event.add_vertex(vertex);
      if ((*blit)->Type()==ATOOLS::btp::Signal_Process) {
        if ((**blit)["NLO_subeventlist"]) {
          THROW(fatal_error,"Events containing correlated subtraction events"
                +std::string(" cannot be translated into the full HepMC event")
                +std::string(" format.\n")
                +std::string("   Try 'EVENT_OUTPUT=HepMC_Short' instead."));
        }
        event.set_signal_process_vertex(vertex);
        if((*blit)->NInP()==2) {
          kf_code fl1=(*blit)->InParticle(0)->Flav().HepEvt();
          kf_code fl2=(*blit)->InParticle(1)->Flav().HepEvt();
          double x1=(*blit)->InParticle(0)->Momentum()[0]/rpa->gen.PBeam(0)[0];
          double x2=(*blit)->InParticle(1)->Momentum()[0]/rpa->gen.PBeam(1)[0];
          double q(0.0), p1(0.0), p2(0.0);
          Blob_Data_Base *facscale((**blit)["Factorisation_Scale"]);
	  if (facscale) q=sqrt(facscale->Get<double>());
	  Blob_Data_Base *xf1((**blit)["XF1"]);
          Blob_Data_Base *xf2((**blit)["XF2"]);
          if (xf1) p1=xf1->Get<double>();
          if (xf2) p2=xf2->Get<double>();
	  HepMC::PdfInfo pdfinfo(fl1, fl2, x1, x2, q, p1, p2);
          event.set_pdf_info(pdfinfo);
        }
      }
      else if ((*blit)->Type()==ATOOLS::btp::Beam || 
	       (*blit)->Type()==ATOOLS::btp::Bunch) {
        for (HepMC::GenVertex::particles_in_const_iterator 
	       pit=vertex->particles_in_const_begin();
             pit!=vertex->particles_in_const_end(); ++pit) {
          if ((*pit)->production_vertex()==NULL) {
            beamparticles.push_back(*pit);
          }
        }
      }
    }
  }
  if (beamparticles.size()==2) {
    event.set_beam_particles(beamparticles[0],beamparticles[1]);
  }
  std::vector<double> weights;
  weights.push_back(weight);
  if (sp) {
    Blob_Data_Base *info((*sp)["MEWeight"]);
    if (!info) THROW(fatal_error,"Missing weight info.");
    double meweight(info->Get<double>());
    weights.push_back(meweight);
    Blob_Data_Base *ofni((*sp)["Weight_Norm"]);
    if (!ofni) THROW(fatal_error,"Missing weight normalisation.");
    double weightnorm(ofni->Get<double>());
    weights.push_back(weightnorm);
    ofni=(*sp)["Trials"];
    if (!ofni) THROW(fatal_error,"Missing nof trials.");
    double trials(ofni->Get<double>());
    weights.push_back(trials);
  }
  event.weights()=weights;
  return true;
}

bool HepMC2_Interface::Sherpa2ShortHepMC(ATOOLS::Blob_List *const blobs,
                                         HepMC::GenEvent& event, double weight)
{
#ifdef USING__HEPMC2__UNITS
  event.use_units(HepMC::Units::GEV,
                  HepMC::Units::MM);
#endif
  event.set_event_number(ATOOLS::rpa->gen.NumberOfGeneratedEvents());
  HepMC::GenVertex * vertex=new HepMC::GenVertex();
  std::vector<HepMC::GenParticle*> beamparticles;
  std::vector<std::pair<HepMC::FourVector,int> > beamparts,
                                                 remnantparts1, remnantparts2;
  Blob *sp(blobs->FindFirst(btp::Signal_Process));
  NLO_subevtlist* subevtlist(NULL);
  ME_wgtinfo* wgtinfo(0);
  
  if (sp) {
    Blob_Data_Base* seinfo=(*sp)["ME_wgtinfo"];
    if (seinfo) wgtinfo=seinfo->Get<ME_wgtinfo*>();
   
    Blob_Data_Base * bdb((*sp)["NLO_subeventlist"]);
    if (bdb) subevtlist=bdb->Get<NLO_subevtlist*>();
  }
  for (ATOOLS::Blob_List::iterator blit=blobs->begin();
       blit!=blobs->end();++blit) {
    Blob* blob=*blit;
    for (int i=0;i<blob->NInP();i++) {
      if (blob->InParticle(i)->ProductionBlob()==NULL) {
        Particle* parton=blob->InParticle(i);
        ATOOLS::Vec4D mom  = parton->Momentum();
        HepMC::FourVector momentum(mom[1],mom[2],mom[3],mom[0]);
        HepMC::GenParticle* inpart = 
	  new HepMC::GenParticle(momentum,parton->Flav().HepEvt(),2);
        vertex->add_particle_in(inpart);
        // distinct because SHRIMPS has no bunches for some reason
        if (blob->Type()==btp::Beam || blob->Type()==btp::Bunch) {
          beamparticles.push_back(inpart);
          beamparts.push_back(std::make_pair(momentum,parton->Flav().HepEvt()));
        }
      }
    }
    for (int i=0;i<blob->NOutP();i++) {
      if (blob->OutParticle(i)->DecayBlob()==NULL) {
        Particle* parton=blob->OutParticle(i);
        ATOOLS::Vec4D mom  = parton->Momentum();
        HepMC::FourVector momentum(mom[1],mom[2],mom[3],mom[0]);
        HepMC::GenParticle* outpart = 
	  new HepMC::GenParticle(momentum,parton->Flav().HepEvt(),1);
        vertex->add_particle_out(outpart);
        if (blob->Type()==btp::Beam) {
          if      (mom[3]>0) 
	    remnantparts1.push_back(std::make_pair(momentum,
						   parton->Flav().HepEvt()));
          else if (mom[3]<0) 
	    remnantparts2.push_back(std::make_pair(momentum,
						   parton->Flav().HepEvt()));
          else THROW(fatal_error,"Ill defined beam remnants.");
        }
      }
    }
    
    if ((*blit)->Type()==ATOOLS::btp::Signal_Process) {
      if((*blit)->NInP()==2) {
        kf_code fl1=(*blit)->InParticle(0)->Flav().HepEvt();
        kf_code fl2=(*blit)->InParticle(1)->Flav().HepEvt();
        double x1=(*blit)->InParticle(0)->Momentum()[0]/rpa->gen.PBeam(0)[0];
        double x2=(*blit)->InParticle(1)->Momentum()[0]/rpa->gen.PBeam(1)[0];
        double q(0.0), p1(0.0), p2(0.0);
        Blob_Data_Base *facscale((**blit)["Factorisation_Scale"]);
        if (facscale) q=sqrt(facscale->Get<double>());
	Blob_Data_Base *xf1((**blit)["XF1"]);
        Blob_Data_Base *xf2((**blit)["XF2"]);
        if (xf1) p1=xf1->Get<double>();
        if (xf2) p2=xf2->Get<double>();
	HepMC::PdfInfo pdfinfo(fl1, fl2, x1, x2, q, p1, p2);
	event.set_pdf_info(pdfinfo);
      }
    }
  }
  event.add_vertex(vertex);
  if (beamparticles.size()==2) {
    event.set_beam_particles(beamparticles[0],beamparticles[1]);
  }
  std::vector<double> weights; weights.push_back(weight);
  if (sp && !subevtlist) {
    Blob_Data_Base *info;
    info=((*sp)["MEWeight"]);
    if (!info) THROW(fatal_error,"Missing weight info.");
    double meweight(info->Get<double>());
    weights.push_back(meweight);
    info=((*sp)["Weight_Norm"]);
    if (!info) THROW(fatal_error,"Missing weight normalisation.");
    double weightnorm(info->Get<double>());
    weights.push_back(weightnorm);
    info=(*sp)["Trials"];
    if (!info) THROW(fatal_error,"Missing nof trials.");
    double trials(info->Get<double>());
    weights.push_back(trials);
    //alphaS value && power
    double rscale2 = (*sp)["Renormalization_Scale"]->Get<double>();
//.........这里部分代码省略.........

bool Rivet_Interface::Run(ATOOLS::Blob_List *const bl)
{
  DEBUG_FUNC("");
  Particle_List pl=bl->ExtractParticles(1);
  for (Particle_List::iterator it=pl.begin(); it!=pl.end(); ++it) {
    if ((*it)->Momentum().Nan()) {
      msg_Error()<<METHOD<<" encountered NaN in momentum. Ignoring event:"
                 <<endl<<*bl<<endl;
      return false;
    }
  }

  double weight(bl->Weight());
  HepMC::GenEvent event;
  if (m_usehepmcshort)  m_hepmc2.Sherpa2ShortHepMC(bl, event, weight);
  else                  m_hepmc2.Sherpa2HepMC(bl, event, weight);
  std::vector<HepMC::GenEvent*> subevents(m_hepmc2.GenSubEventList());
#ifdef HEPMC_HAS_CROSS_SECTION
  // leave this, although will be overwritten later
  HepMC::GenCrossSection xs;
  xs.set_cross_section(p_eventhandler->TotalXS(), p_eventhandler->TotalErr());
  event.set_cross_section(xs);
  for (size_t i(0);i<subevents.size();++i) {
    subevents[i]->set_cross_section(xs);
  }
#endif

  // 1st event build index map, thereafter only lookup
  ExtractVariations(event,subevents);
  for (RivetScaleVariationMap::iterator it(m_rivet.begin());
       it!=m_rivet.end();++it) {
    msg_Debugging()<<"Running rivet for "<<it->first<<" with "
                   <<it->second->RivetMap().size()<<" histograms."<<std::endl;
    Rivet_Map& rivetmap(it->second->RivetMap());
    if (subevents.size()) {
      it->second->SynchroniseCrossSection();
      for (size_t i(0);i<subevents.size();++i) {
        SetEventWeight(it->second,*subevents[i],i);
        GetRivet(rivetmap,"", 0)->analyze(*subevents[i]);
      }
    }
    else {
      SetEventWeight(it->second,event);
      GetRivet(rivetmap,"",0)->analyze(event);
      Blob *sp(bl->FindFirst(btp::Signal_Process));
      size_t parts=0;
      if (sp) {
        std::string multi(sp?sp->TypeSpec():"");
        if (multi[3]=='_') multi=multi.substr(2,1);
        else multi=multi.substr(2,2);
        parts=ToType<size_t>(multi);
      }
      if (m_splitjetconts && sp) {
        GetRivet(rivetmap,"",parts)->analyze(event);
      }
      if (m_splitcoreprocs && sp) {
        GetRivet(rivetmap,GetCoreProc(sp->TypeSpec()),0)
            ->analyze(event);
        if (m_splitjetconts) {
          GetRivet(rivetmap,GetCoreProc(sp->TypeSpec()),parts)
              ->analyze(event);
        }
      }
      if (m_splitSH && sp) {
        std::string typespec=sp->TypeSpec();
        typespec=typespec.substr(typespec.length()-2, 2);
        std::string type="";
        if (typespec=="+S") type="S";
        else if (typespec=="+H") type="H";

        if (type!="") {
          GetRivet(rivetmap,type,0)->analyze(event);
          if (m_splitjetconts) {
            GetRivet(rivetmap,type,parts)->analyze(event);
          }
        }
      }
    }
  }
  if (subevents.size()) {
    ResetRivetScaleVariationMapRSWeights();
    m_hepmc2.DeleteGenSubEventList();
  }

  ++m_nevt;
  for (RivetScaleVariationMap::iterator it(m_rivet.begin());
       it!=m_rivet.end();++it) {
    msg_Debugging()<<"Checking rivet for "<<it->first<<" with "
                   <<it->second->RivetMap().size()<<" histograms."<<std::endl;
  }
  if (m_evtbyevtxs) {
    for (RivetScaleVariationMap::iterator mit(m_rivet.begin());
         mit!=m_rivet.end();++mit) {
      std::string out=m_outpath;
      if (mit->first!="" && mit->first!="nominal") out += "."+mit->first;
      std::string namestr(m_rivet.size()>1?" for "+mit->first:"");
      std::string output(std::string("**  Total XS")+namestr
                         +std::string(" = ( ")
                         +ToString(mit->second->TotalXS(),m_xsoutputprecision)
                         +std::string(" +- ")
//.........这里部分代码省略.........

void Rivet_Interface::ExtractVariations(const HepMC::GenEvent& evt)
{
  DEBUG_FUNC("");
  const HepMC::WeightContainer& wc(evt.weights());
  std::map<std::string,double> wgtmap;
  double ntrials(1.);
  size_t xstype(0);
#ifdef HEPMC_HAS_NAMED_WEIGHTS
#ifdef HEPMC_HAS_WORKING_NAMED_WEIGHTS // replace by final HepMC-2.07 variable
  std::vector<std::string> keys(wc.keys());
  msg_Debugging()<<keys<<std::endl;
  for (size_t i(0);i<keys.size();++i) {
    std::string cur(keys[i]);
    if (m_splitvariations && cur.find("MUR")!=std::string::npos &&
                             cur.find("MUF")!=std::string::npos &&
                             cur.find("PDF")!=std::string::npos) {
      wgtmap[cur]=wc[cur];
    }
    else if (cur=="Weight")  wgtmap["nominal"]=wc[cur];
    else if (cur=="NTrials") ntrials=wc[cur];
    else if (cur=="Reweight_Type" && wc[cur]&64) xstype=1;
  }
#else
  // lookup all evt-wgts with name "MUR<fac>_MUF<fac>_PDF<id>"
  // at the moment the only way to do that is to filter the printout
  // accuracy limited to print out accu of 6 digits, must suffice
  MyStrStream str;
  str.precision(m_hepmcoutputprecision);
  wc.print(str);

  // need a temp object first, as we need to get ntrials first
  while (str) {
    double wgt(0.);
    std::string cur("");
    str>>cur;
    if (cur.length()==0) continue;
    // weight is between "," and trailing bracket
    // name is between leading bracket and ","
    wgt=ToType<double>(cur.substr(cur.find(",")+1,cur.find(")")-1));
    cur=cur.substr(1,cur.find(",")-1);
    if (m_splitvariations && cur.find("MUR")!=std::string::npos &&
                             cur.find("MUF")!=std::string::npos &&
                             cur.find("PDF")!=std::string::npos) {
      wgtmap[cur]=wgt;
    }
    else if (cur=="Weight")  wgtmap["nominal"]=wgt;
    else if (cur=="NTrials") ntrials=wgt;
    else if (cur=="Reweight_Type" && ((int)wgt)&64) xstype=1;
  }
#endif /* HEPMC_HAS_WORKING_NAMED_WEIGHTS */
#else
  wgtmap["nominal"]=wc[0];
  ntrials=wc[3];
  xstype=(((wc.size()==5&&((int)wc[4]&64))||(wc.size()==11&&((int)wc[10]&64)))?1:0);
#endif /* HEPMC_HAS_NAMED_WEIGHTS */
  if (msg_LevelIsDebugging()) {
    for (std::map<std::string,double>::iterator wit(wgtmap.begin());
         wit!=wgtmap.end();++wit)
      msg_Out()<<wit->first<<" : "<<wit->second<<std::endl;
  }
  // now construct or fill into the scale variation map
  for (std::map<std::string,double>::iterator wit(wgtmap.begin());
       wit!=wgtmap.end();++wit) {
    RivetScaleVariationMap::iterator rit=m_rivet.find(wit->first);
    if (rit==m_rivet.end()) {
      msg_Debugging()<<"creating new entry in m_rivet"<<std::endl;
      m_rivet[wit->first]