C++ GenEvent::weights方法代码示例

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 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 int multiLHC = 3000;
  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,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::BIpt_scheme, fastjet::Best);
//=============================================================================

  TRandom3 *r3 = new TRandom3(0);
  TF1 *fBkg = BackgroundSpec();
//=============================================================================

  std::vector<fastjet::PseudoJet> fjInputVac;
  std::vector<fastjet::PseudoJet> fjInputHyd;
//=============================================================================

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


//=============================================================================

  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) {
    fjInputVac.resize(0);
    fjInputHyd.resize(0);

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

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

      if ((TMath::Abs(vPar.Eta())<dCutEtaMax)) {
        fjInputVac.push_back(fastjet::PseudoJet(vPseudo.Px(), vPseudo.Py(), vPseudo.Pz(), vPseudo.E()));
        fjInputVac.back().set_user_info(new UserInfoTrk(false));
        fjInputVac.back().set_user_index(iCount); iCount+=1;
      }
    }
//=============================================================================

    for (int i=0; i<=multiLHC; i++) {
//.........这里部分代码省略.........

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]=new Rivet_Scale_Variation(wit->first);
      m_rivet[wit->first]->AddPoint(wit->second,ntrials,xstype);
    }
    else rit->second->AddPoint(wit->second,ntrials,xstype);
  }
  if (msg_LevelIsDebugging()) {
    for (RivetScaleVariationMap::iterator rit(m_rivet.begin());
         rit!=m_rivet.end();++rit)
      msg_Out()<<rit->first<<" : "<<rit->second->Weight()<<std::endl;
  }
}

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>();
//.........这里部分代码省略.........

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,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::antikt_algorithm, dSjeR, fastjet::BIpt_scheme, fastjet::Best);
//=============================================================================

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

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

  TH1D *hJet = new TH1D("hJet", "", 1000, 0., 1000.); hJet->Sumw2(); list->Add(hJet);
  TH2D *hJetNsj = new TH2D("hJetNsj", "", 1000, 0., 1000., 101, -0.5, 100.5); hJetNsj->Sumw2(); list->Add(hJetNsj);

  TH2D *hJetIsj = new TH2D("hJetIsj", "", 1000, 0., 1000., 1000, 0., 1000.); hJetIsj->Sumw2(); list->Add(hJetIsj);
  TH2D *hJet1sj = new TH2D("hJet1sj", "", 1000, 0., 1000., 1000, 0., 1000.); hJet1sj->Sumw2(); list->Add(hJet1sj);
  TH2D *hJet2sj = new TH2D("hJet2sj", "", 1000, 0., 1000., 1000, 0., 1000.); hJet2sj->Sumw2(); list->Add(hJet2sj);
  TH2D *hJetDsj = new TH2D("hJetDsj", "", 1000, 0., 1000., 1000, 0., 1000.); hJetDsj->Sumw2(); list->Add(hJetDsj);

  TH2D *hJetIsz = new TH2D("hJetIsz", "", 1000, 0., 1000., 120, 0., 1.2); hJetIsz->Sumw2(); list->Add(hJetIsz);
  TH2D *hJet1sz = new TH2D("hJet1sz", "", 1000, 0., 1000., 120, 0., 1.2); hJet1sz->Sumw2(); list->Add(hJet1sz);
  TH2D *hJet2sz = new TH2D("hJet2sz", "", 1000, 0., 1000., 120, 0., 1.2); hJet2sz->Sumw2(); list->Add(hJet2sz);
  TH2D *hJetDsz = new TH2D("hJetDsz", "", 1000, 0., 1000., 120, 0., 1.2); hJetDsz->Sumw2(); list->Add(hJetDsz);
//=============================================================================

  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 dXsect  = evt->cross_section()->cross_section() / 1e9;
    double dWeight = evt->weights().back();
    double dNorm = dWeight * dXsect;
    hWeightSum->Fill(0.5, dWeight);

    TVector3 vPar;
    for (HepMC::GenEvent::particle_const_iterator p=evt->particles_begin(); p!=evt->particles_end(); ++p) if ((*p)->status()==1) {
      vPar.SetXYZ((*p)->momentum().px(), (*p)->momentum().py(), (*p)->momentum().pz());

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

    fastjet::ClusterSequenceArea clustSeq(fjInput, jetsDef, areaDef);
//.........这里部分代码省略.........

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()));

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