Dear experts,
I am trying to create the vector tree within Pythia 8.
My code is working fine with low statistics, however when I increase the number of events I faced following error message:
libc++abi.dylib: terminating with uncaught exception of type std::out_of_range: vector
Would you please help me to solve the problem?
Please see the below code, I am using.
// main72.cc is a part of the PYTHIA event generator.
// Copyright (C) 2016 Torbjorn Sjostrand.
// PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
// Please respect the MCnet Guidelines, see GUIDELINES for details.
// This is a simple test program.
// It compares SlowJet, FJcore and FastJet, showing that they
// find the same jets.
#include <iostream>
#include "TApplication.h"
#include "TTree.h"
#include <vector>
#include "TFile.h"
#include "Pythia8/Pythia.h"
#include "TROOT.h"
#include "TVirtualPad.h"
#include "TBranch.h"
#include <fstream>
// The FastJet3.h header enables automatic initialisation of
// fastjet::PseudoJet objects from Pythia8 Particle and Vec4 objects,
// as well as advanced features such as access to (a copy of)
// the original Pythia 8 Particle directly from the PseudoJet,
// and fastjet selectors that make use of the Particle properties.
// See the extensive comments in the header file for further details
// and examples.
#include "Pythia8Plugins/FastJet3.h"
using namespace Pythia8;
using namespace std;
vector<double> *id = new vector<double>;
vector<double> *events = new vector<double>;
vector<double> *nJetsS = new vector<double>;
vector<double> *jetpt = new vector<double>;
vector<double> *jetrapidity= new vector<double>;
vector<double> *jetphi= new vector<double>;
int main() {
// Number of events, generated and listed ones (for jets).
int nEvent = 1500000;
int nListJets = 3;
double eProton = 50000.;
double eElectron = 60;
double Q2min = 25.;
//double nFast=0.;
double PT, ETA, PHI=0.;
// Select common parameters for SlowJet and FastJet analyses.
int power = -1; // -1 = anti-kT; 0 = C/A; 1 = kT.
double R = 0.5; // Jet size.
double pTMin = 5.0; // Min jet pT.
double etaMax = 7.0; // Pseudorapidity range of detector.
int select = 2; // Which particles are included?
int massSet = 2; // Which mass are they assumed to have?
// Generator. Shorthand for event.
Pythia pythia;
Event& event = pythia.event;
// Beam parameters.
/*
pythia.readString("Beams:idA = 11");
pythia.readString("Beams:idB = 2212");
pythia.readString("Beams:eA = 60.");
pythia.readString("Beams:eB = 50000.");
*/
pythia.readString("Beams:frameType = 2");
// BeamA = proton.
pythia.readString("Beams:idA = 2212");
pythia.settings.parm("Beams:eA", eProton);
// BeamB = electron.
pythia.readString("Beams:idB = 11");
pythia.settings.parm("Beams:eB", eElectron);
//pythia.readString("Beams:eCM = 100000.");
pythia.readString("Random:setSeed = on");
pythia.readString("Random:seed = 10000002");
// Process selection.
/// pythia.readString("HardQCD:all = on");
/// pythia.readString("PhaseSpace:pTHatMin = 10.");
pythia.readString("WeakBosonExchange:ff2ff(t:gmZ) = on");
pythia.readString("WeakBosonExchange:ff2ff(t:W) = on");
// No event record printout.
pythia.readString("Next:numberShowInfo = 0");
pythia.readString("Next:numberShowProcess = 0");
pythia.readString("Next:numberShowEvent = 0");
// Phase-space cut: minimal Q2 of process.
pythia.settings.parm("PhaseSpace:Q2Min", Q2min);
// Set dipole recoil on. Necessary for DIS + shower.
pythia.readString("SpaceShower:dipoleRecoil = on");
// Allow emissions up to the kinematical limit,
// since rate known to match well to matrix elements everywhere.
pythia.readString("SpaceShower:pTmaxMatch = 2");
// QED radiation off lepton not handled yet by the new procedure.
pythia.readString("PDF:lepton = off");
pythia.readString("TimeShower:QEDshowerByL = off");
pythia.init();
// Set up FastJet jet finder.
// one can use either explicitly use antikt, cambridge, etc., or
// just use genkt_algorithm with specification of power
//fastjet::JetAlgorithm algorithm;
//if (power == -1) algorithm = fastjet::antikt_algorithm;
//if (power == 0) algorithm = fastjet::cambridge_algorithm;
//if (power == 1) algorithm = fastjet::kt_algorithm;
//fastjet::JetDefinition jetDef(algorithm, R);
// there's no need for a pointer to the jetDef (it's a fairly small object)
fastjet::JetDefinition jetDef(fastjet::genkt_algorithm, R, power);
std::vector <fastjet::PseudoJet> fjInputs;
gROOT->Reset();
TFile *file = new TFile("mn_fastjets_tree_FCCenergy_1p5M_seed02.root","recreate"); //root
TTree *t1 = new TTree ("Tree","Tree"); // Tree
//Long64_t fgMaxTreeSize = 100*1024*1024;//set file size
//t1->SetMaxTreeSize(fgMaxTreeSize);// setting maximum tree size
//Branches
t1->Branch("id",&id);
//event number
t1->Branch("events",&events);
t1->Branch("nJetsS",&nJetsS);
t1->Branch("jetpt",&jetpt);
t1->Branch("jetrapidity",&jetrapidity);
t1->Branch("jetphi",&jetphi);
// Begin event loop. Generate event. Skip if error.
for(int iEvent=0;iEvent<=nEvent;iEvent++){
//clock_t befGen = clock();
if (!pythia.next()) continue;
//clock_t aftGen = clock();
events->push_back(iEvent);
//cout<<"Event tree is filled"<<endl;
// Begin FastJet analysis: extract particles from event record.
//clock_t befFast = clock();
fjInputs.resize(0);
Vec4 pTemp;
double mTemp;
int nAnalyze = 0;
for (int i = 0; i < event.size(); ++i) if (event[i].isFinal()) {
// Require visible/charged particles inside detector.
if (select > 2 && event[i].isNeutral() ) continue;
else if (select == 2 && !event[i].isVisible() ) continue;
if (etaMax < 20. && abs(event[i].eta()) > etaMax) continue;
// Create a PseudoJet from the complete Pythia particle.
fastjet::PseudoJet particleTemp = event[i];
// Optionally modify mass and energy.
pTemp = event[i].p();
mTemp = event[i].m();
if (massSet < 2) {
mTemp = (massSet == 0 || event[i].id() == 22) ? 0. : 0.13957;
pTemp.e( sqrt(pTemp.pAbs2() + mTemp*mTemp) );
particleTemp.reset_momentum( pTemp.px(), pTemp.py(),
pTemp.pz(), pTemp.e() );
}
// Store acceptable particles as input to Fastjet.
// Conversion to PseudoJet is performed automatically
// with the help of the code in FastJet3.h.
fjInputs.push_back( particleTemp);
++nAnalyze;
}
// Run Fastjet algorithm and sort jets in pT order.
vector <fastjet::PseudoJet> inclusiveJets, sortedJets;
fastjet::ClusterSequence clustSeq(fjInputs, jetDef);
inclusiveJets = clustSeq.inclusive_jets(pTMin);
sortedJets = sorted_by_pt(inclusiveJets);
//clock_t aftFast = clock();
// List first few FastJet jets and some info about them.
// Note: the final few columns are illustrative of what information
// can be extracted, but does not exhaust the possibilities.
if (iEvent < nListJets) {
cout << "\n -------- FastJet jets, p = " << setw(2) << power
<< " --------------------------------------------------\n\n "
<< " i pT y phi mult chgmult photons"
<< " hardest pT in neutral " << endl
<< " "
<< " constituent hadrons " << endl;
for (int i = 0; i < int(sortedJets.size()); ++i) {
vector<fastjet::PseudoJet> constituents
= sortedJets[i].constituents();
fastjet::PseudoJet hardest
= fastjet::SelectorNHardest(1)(constituents)[0];
vector<fastjet::PseudoJet> neutral_hadrons
= ( fastjet::SelectorIsHadron()
&& fastjet::SelectorIsNeutral())(constituents);
double neutral_hadrons_pt = join(neutral_hadrons).perp();
cout << setw(4) << i << fixed << setprecision(3) << setw(11)
<< sortedJets[i].perp() << setw(9) << sortedJets[i].rap()
<< setw(9) << sortedJets[i].phi_std()
<< setw(6) << constituents.size()
<< setw(8) << fastjet::SelectorIsCharged().count(constituents)
<< setw(8) << fastjet::SelectorId(22).count(constituents)
<< setw(13) << hardest.user_info<Particle>().name()
<< " " << setw(10) << neutral_hadrons_pt << endl;
}
cout << "\n -------- End FastJet Listing ------------------"
<< "---------------------------------" << endl;
}
int eventNo=-1;
int jetNo=-1;
for(int i = 0; i < pythia.event.size(); i++){
// cout<<"\r"<<double(iEvent)/nEvent*100<<"%"<<flush;//progress
if (sortedJets.empty()) continue;
id->push_back(pythia.event[i].id());
eventNo=i;
// Fill distribution of FastJet jets .
int nFast = sortedJets.size();
nJetsS->push_back(sortedJets.size());
for (int j = 0; j < nFast; ++j) {
jetNo=j;
jetpt->push_back(sortedJets[j].perp());
PT=sortedJets[j].perp();
jetrapidity->push_back(sortedJets[j].rap());
ETA=sortedJets[j].rap();
jetphi->push_back(sortedJets[j].phi_std());
PHI=sortedJets[j].phi_std();
}
}//eventsize
//cout<<"event no ="<<eventNo<<" jet no ="<<jetNo<<" jetpt ="<<PT<<" jetrapidity ="<<ETA<<" jetphi ="<<PHI<<endl;
t1->Fill();
id->clear();
events->clear();
nJetsS->clear();
jetpt->clear();
jetrapidity->clear();
jetphi->clear();
//cout<<" jet size ="<<sortedJets.size()<<endl;
} //event loop
// Statistics on event generation.
pythia.stat();
file = t1->GetCurrentFile();// splitting root file
t1->Write();
file->Close();
return 0;
}
_ROOT Version:ROOT 6.06/06