From af7c9208184b7474dd4d6b860e568d02690ebd3f Mon Sep 17 00:00:00 2001 From: Marta Verweij Date: Fri, 26 Feb 2021 09:46:36 +0100 Subject: [PATCH] add gluon splitting to events + analyze --- Makefile | 20 +++- README_ForBScStudents.md | 2 +- include/pythiaEvent.hh | 23 ++++- include/vinciaEvent.hh | 93 ++++++++++++++++++ runGluonSplittingAnalysis.cc | 185 +++++++++++++++++++++++++++++++++++ 5 files changed, 316 insertions(+), 7 deletions(-) create mode 100644 include/vinciaEvent.hh create mode 100644 runGluonSplittingAnalysis.cc diff --git a/Makefile b/Makefile index 9b1c583..3bc9938 100644 --- a/Makefile +++ b/Makefile @@ -33,14 +33,14 @@ COMMONSRC = F77SRC = COMMONOBJ = -PROGSRC = runConstituentSubtraction.cc runConversionQPYTHIA.cc runCreatePythiaEvents.cc runCreatePythiaEventsPartonLevel.cc runCreatePythiaEventsWJet.cc runCreateThermalEvents.cc runCreateVinciaEvents.cc runCSVariations.cc runFromFile.cc runHybridJetAnalysis.cc runJetPerformance.cc runJetProfile.cc runJewelSub.cc runSDGenVarious.cc runSDGenVariousJewelSub.cc runSharedLayerSubtraction.cc runSimpleJetAnalysis.cc runSoftDrop.cc runtest.cc -PROGOBJ = runConstituentSubtraction.o runConversionQPYTHIA.o runCreatePythiaEvents.o runCreatePythiaEventsPartonLevel.o runCreatePythiaEventsWJet.o runCreateThermalEvents.o runCreateVinciaEvents.o runCSVariations.o runFromFile.o runHybridJetAnalysis.o runJetPerformance.o runJetProfile.o runJewelSub.o runSDGenVarious.o runSDGenVariousJewelSub.o runSharedLayerSubtraction.o runSimpleJetAnalysis.o runSoftDrop.o runtest.o +PROGSRC = runConstituentSubtraction.cc runConversionQPYTHIA.cc runCreatePythiaEvents.cc runCreatePythiaEventsPartonLevel.cc runCreatePythiaEventsWJet.cc runCreateThermalEvents.cc runCreateVinciaEvents.cc runCSVariations.cc runFromFile.cc runGluonSplittingAnalysis.cc runHybridJetAnalysis.cc runJetPerformance.cc runJetProfile.cc runJewelSub.cc runSDGenVarious.cc runSDGenVariousJewelSub.cc runSharedLayerSubtraction.cc runSimpleJetAnalysis.cc runSoftDrop.cc runtest.cc +PROGOBJ = runConstituentSubtraction.o runConversionQPYTHIA.o runCreatePythiaEvents.o runCreatePythiaEventsPartonLevel.o runCreatePythiaEventsWJet.o runCreateThermalEvents.o runCreateVinciaEvents.o runCSVariations.o runFromFile.o runGluonSplittingAnalysis.o runHybridJetAnalysis.o runJetPerformance.o runJetProfile.o runJewelSub.o runSDGenVarious.o runSDGenVariousJewelSub.o runSharedLayerSubtraction.o runSimpleJetAnalysis.o runSoftDrop.o runtest.o INCLUDE += LIBRARIES += -LPU14 -lPU14 -lz -all: runConstituentSubtraction runConversionQPYTHIA runCreatePythiaEvents runCreatePythiaEventsPartonLevel runCreatePythiaEventsWJet runCreateThermalEvents runCreateVinciaEvents runCSVariations runFromFile runHybridJetAnalysis runJetPerformance runJetProfile runJewelSub runSDGenVarious runSDGenVariousJewelSub runSharedLayerSubtraction runSimpleJetAnalysis runSoftDrop runtest +all: runConstituentSubtraction runConversionQPYTHIA runCreatePythiaEvents runCreatePythiaEventsPartonLevel runCreatePythiaEventsWJet runCreateThermalEvents runCreateVinciaEvents runCSVariations runFromFile runGluonSplittingAnalysis runHybridJetAnalysis runJetPerformance runJetProfile runJewelSub runSDGenVarious runSDGenVariousJewelSub runSharedLayerSubtraction runSimpleJetAnalysis runSoftDrop runtest runConstituentSubtraction: runConstituentSubtraction.o $(COMMONOBJ) @@ -70,6 +70,9 @@ runCSVariations: runCSVariations.o $(COMMONOBJ) runFromFile: runFromFile.o $(COMMONOBJ) $(CXX) $(LDFLAGS) -o $@ $@.o $(COMMONOBJ) $(LIBRARIES) +runGluonSplittingAnalysis: runGluonSplittingAnalysis.o $(COMMONOBJ) + $(CXX) $(LDFLAGS) -o $@ $@.o $(COMMONOBJ) $(LIBRARIES) + runHybridJetAnalysis: runHybridJetAnalysis.o $(COMMONOBJ) $(CXX) $(LDFLAGS) -o $@ $@.o $(COMMONOBJ) $(LIBRARIES) @@ -108,7 +111,7 @@ clean: rm -vf $(COMMONOBJ) $(PROGOBJ) realclean: clean - rm -vf runConstituentSubtraction runConversionQPYTHIA runCreatePythiaEvents runCreatePythiaEventsPartonLevel runCreatePythiaEventsWJet runCreateThermalEvents runCreateVinciaEvents runCSVariations runFromFile runHybridJetAnalysis runJetPerformance runJetProfile runJewelSub runSDGenVarious runSDGenVariousJewelSub runSharedLayerSubtraction runSimpleJetAnalysis runSoftDrop runtest + rm -vf runConstituentSubtraction runConversionQPYTHIA runCreatePythiaEvents runCreatePythiaEventsPartonLevel runCreatePythiaEventsWJet runCreateThermalEvents runCreateVinciaEvents runCSVariations runFromFile runGluonSplittingAnalysis runHybridJetAnalysis runJetPerformance runJetProfile runJewelSub runSDGenVarious runSDGenVariousJewelSub runSharedLayerSubtraction runSimpleJetAnalysis runSoftDrop runtest .cc.o: $< $(CXX) $(CXXFLAGS) $(INCLUDE) -c $< -o $@ @@ -171,6 +174,15 @@ runFromFile.o: include/softDropGroomer.hh include/jetCollection.hh runFromFile.o: include/jewelMatcher.hh include/treeWriter.hh runFromFile.o: include/jetMatcher.hh include/randomCones.hh runFromFile.o: include/Angularity.hh include/jewelMatcher.hh +runGluonSplittingAnalysis.o: include/ProgressBar.h PU14/EventMixer.hh +runGluonSplittingAnalysis.o: PU14/CmdLine.hh PU14/EventSource.hh +runGluonSplittingAnalysis.o: PU14/CmdLine.hh PU14/PU14.hh +runGluonSplittingAnalysis.o: PU14/HepPID/ParticleIDMethods.hh +runGluonSplittingAnalysis.o: include/extraInfo.hh include/jetCollection.hh +runGluonSplittingAnalysis.o: include/softDropGroomer.hh PU14/PU14.hh +runGluonSplittingAnalysis.o: include/jetCollection.hh include/jewelMatcher.hh +runGluonSplittingAnalysis.o: include/treeWriter.hh include/jetMatcher.hh +runGluonSplittingAnalysis.o: include/Angularity.hh runHybridJetAnalysis.o: include/ProgressBar.h PU14/EventMixer.hh runHybridJetAnalysis.o: PU14/CmdLine.hh PU14/EventSource.hh PU14/CmdLine.hh runHybridJetAnalysis.o: PU14/PU14.hh PU14/HepPID/ParticleIDMethods.hh diff --git a/README_ForBScStudents.md b/README_ForBScStudents.md index 863daca..9a80151 100644 --- a/README_ForBScStudents.md +++ b/README_ForBScStudents.md @@ -184,7 +184,7 @@ More info on the quark cluster can be found here: https://uugrasp.github.io/UUCo Event samples can be found in the jet quenching CERNBOX: * From lxplus (CERN account required): /eos/project/j/jetquenching/www * Webbrowser CERNBOX (CERN account required): https://cernbox.cern.ch/index.php/s/kRy9M7NC9iilE9Z -* Webbrowser (publicly accessible): http://jetquenchingtools.web.cern.ch/jetquenchingtools/ (You can use wget and curl on this) +* Webbrowser (publicly accessible): https://jetquenchingtools.web.cern.ch/JetQuenchingTools/ (You can use wget and curl on this) * Mount eos on a laptop or local desktop (CERN account required): https://cern.service-now.com/service-portal/article.do?n=KB0003493 You will find samples from various event generators. For underlying event we have: 'thermal' which is independent particle production using a Boltzmann distribution with a fixed multiplicity and mean pT (indicated in the file names). For the hard signal we have PYTHIA8 and JEWEL events with various pT,hat settings. diff --git a/include/pythiaEvent.hh b/include/pythiaEvent.hh index 3097be6..d0d04b1 100644 --- a/include/pythiaEvent.hh +++ b/include/pythiaEvent.hh @@ -76,18 +76,37 @@ std::vector pythiaEvent::createPythiaEvent() { partons.clear(); //empty list before storing partons of new event for (int i = 0; i < pythia.event.size(); ++i) { - if (pythia.event[i].isFinal()) { + if (pythia.event[i].isFinal()) { //all final state particles fastjet::PseudoJet p(pythia.event[i].px(),pythia.event[i].py(),pythia.event[i].pz(),pythia.event[i].e()); p.set_user_info(new extraInfo(pythia.event[i].id(), 0)); if(p.rap()>rapMin_ && p.rap() +#include +#include +#include +#include + +#include "Pythia8/Pythia.h" + +#include "extraInfo.hh" + +//using namespace std; + +//--------------------------------------------------------------- +// Description +// This class generates a pythia8 event +// Author: M. Verweij +//--------------------------------------------------------------- + +class vinciaEvent { + +private : + Pythia8::Pythia pythia; + double pthat_; + unsigned int tune_; + double rapMin_; + double rapMax_; + bool partonLevel_; + + std::vector partons; + +public : + vinciaEvent(double pthat = 120., unsigned int tune = 14, double rapMin = -3., double rapMax = 3., bool partonLevel = false); + std::vector createPythiaEvent(); + + std::vector getPartonList() const { return partons; } + + void getStat() {pythia.stat();} + +}; + +vinciaEvent::vinciaEvent(double pthat, unsigned int tune, double rapMin, double rapMax, bool partonLevel) : + pthat_(pthat), tune_(tune), rapMin_(rapMin), rapMax_(rapMax), partonLevel_(partonLevel) +{ + + // Generator. LHC process and output selection. Initialization. + // tunes: http://home.thep.lu.se/~torbjorn/pythia82html/Tunes.html + pythia.readString("Beams:eCM = 5002."); + pythia.readString("HardQCD:all = on"); + pythia.readString(Form("PhaseSpace:pTHatMin = %.1f",pthat_)); + pythia.readString("Next:numberShowInfo = 0"); + pythia.readString("Next:numberShowProcess = 0"); + pythia.readString("Next:numberShowEvent = 0"); + pythia.readString("PartonLevel:MPI = on"); + //pythia.readString(Form("Tune:pp = %d",tune_)); //MV: I think this is just relevant for pythia and not vincia + pythia.readString("Random:setSeed = on"); + pythia.readString("Random:seed = 0"); + if(partonLevel_) { + pythia.readString("HadronLevel:all = off"); + } + pythia.readString("PartonShowers:Model = 2"); //activate the VINCIA parton shower + + pythia.init(); + +} + +std::vector vinciaEvent::createPythiaEvent() { + + pythia.next(); //generate next event + + std::vector particles; + partons.clear(); //empty list before storing partons of new event + + for (int i = 0; i < pythia.event.size(); ++i) { + if (pythia.event[i].isFinal()) { + fastjet::PseudoJet p(pythia.event[i].px(),pythia.event[i].py(),pythia.event[i].pz(),pythia.event[i].e()); + p.set_user_info(new extraInfo(pythia.event[i].id(), 0)); + if(p.rap()>rapMin_ && p.rap() +#include + +#include "TFile.h" +#include "TTree.h" + +#include "fastjet/PseudoJet.hh" +#include "fastjet/ClusterSequenceArea.hh" + +#include "include/ProgressBar.h" + +#include "PU14/EventMixer.hh" +#include "PU14/CmdLine.hh" +#include "PU14/PU14.hh" + +#include "include/extraInfo.hh" +#include "include/jetCollection.hh" +#include "include/softDropGroomer.hh" +#include "include/treeWriter.hh" +#include "include/jetMatcher.hh" +#include "include/Angularity.hh" + +using namespace std; +using namespace fastjet; + +// ./runGluonSplittingAnalysis -hard /Users/mverweij/mnt/eos/project/j/jetquenching/JetWorkshop2017/samples/pythia8/dijet120/PythiaEventsTune14PtHat120_0.pu14 -nev 10 + +int main (int argc, char ** argv) { + + auto start_time = chrono::steady_clock::now(); + + CmdLine cmdline(argc,argv); + // inputs read from command line + int nEvent = cmdline.value("-nev",1); // first argument: command line option; second argument: default value + //bool verbose = cmdline.present("-verbose"); + + cout << "will run on " << nEvent << " events" << endl; + + // Uncomment to silence fastjet banner + ClusterSequence::set_fastjet_banner_stream(NULL); + + //to write info to root tree + treeWriter trw("jetTree"); + + //Jet definition + double R = 0.4; + double ghostRapMax = 6.0; + double ghost_area = 0.005; + int active_area_repeats = 1; + GhostedAreaSpec ghost_spec(ghostRapMax, active_area_repeats, ghost_area); + AreaDefinition area_def = AreaDefinition(active_area,ghost_spec); + JetDefinition jet_def(antikt_algorithm, R); + + double jetRapMax = 3.0; + Selector jet_selector = SelectorAbsRapMax(jetRapMax); + + Angularity width(1.,1.,R); + Angularity pTD(0.,2.,R); + + ProgressBar Bar(cout, nEvent); + Bar.SetStyle(-1); + + EventMixer mixer(&cmdline); //the mixing machinery from PU14 workshop + + // loop over events + int iev = 0; + unsigned int entryDiv = (nEvent > 200) ? nEvent / 200 : 1; + while ( mixer.next_event() && iev < nEvent ) + { + // increment event number + iev++; + + Bar.Update(iev); + Bar.PrintWithMod(entryDiv); + + vector particlesMergedAll = mixer.particles(); + + vector eventWeight; + eventWeight.push_back(mixer.hard_weight()); + eventWeight.push_back(mixer.pu_weight()); + + // extract hard partons that initiated the jets + fastjet::Selector parton_selector = SelectorVertexNumber(-1); + vector partons = parton_selector(particlesMergedAll); + + // extract hard partons from first splitting + fastjet::Selector parton_selector_split = SelectorVertexNumber(-2); + vector partonsFirstSplit = parton_selector_split(particlesMergedAll); + + // select final state particles from hard event only + fastjet::Selector sig_selector = SelectorVertexNumber(0); + vector particlesSig = sig_selector(particlesMergedAll); + + // select final state particles from background event only + fastjet::Selector bkg_selector = SelectorVertexNumber(1); + vector particlesBkg = bkg_selector(particlesMergedAll); + + vector particlesMerged = particlesBkg; + particlesMerged.insert( particlesMerged.end(), particlesSig.begin(), particlesSig.end() ); + + //--------------------------------------------------------------------------- + // look at first splitting of hard partons + //--------------------------------------------------------------------------- + std::vector drsplit; + std::vector drmom1; + std::vector drmom2; + int id = 0; + for(int ip = 0; ip widthSig; widthSig.reserve(jetCollectionSig.getJet().size()); + vector pTDSig; pTDSig.reserve(jetCollectionSig.getJet().size()); + for(PseudoJet jet : jetCollectionSig.getJet()) { + widthSig.push_back(width.result(jet)); + pTDSig.push_back(pTD.result(jet)); + } + jetCollectionSig.addVector("widthSig", widthSig); + jetCollectionSig.addVector("pTDSig", pTDSig); + + //--------------------------------------------------------------------------- + // Groom the jets + //--------------------------------------------------------------------------- + + //SoftDrop grooming classic for signal jets (zcut=0.1, beta=0) + softDropGroomer sdgSigBeta00Z01(0.1, 0.0, R); + jetCollection jetCollectionSigSDBeta00Z01(sdgSigBeta00Z01.doGrooming(jetCollectionSig)); + jetCollectionSigSDBeta00Z01.addVector("zgSigSDBeta00Z01", sdgSigBeta00Z01.getZgs()); + jetCollectionSigSDBeta00Z01.addVector("ndropSigSDBeta00Z01", sdgSigBeta00Z01.getNDroppedSubjets()); + jetCollectionSigSDBeta00Z01.addVector("dr12SigSDBeta00Z01", sdgSigBeta00Z01.getDR12()); + + + //--------------------------------------------------------------------------- + // write tree + //--------------------------------------------------------------------------- + + //Give variable we want to write out to treeWriter. + //Only vectors of the types 'jetCollection', and 'double', 'int', 'PseudoJet' are supported + + trw.addCollection("eventWeight", eventWeight); + trw.addPartonCollection("partons", partons); + trw.addPartonCollection("partonsFirstSplit", partonsFirstSplit); + trw.addDoubleCollection("drSplit", drsplit); + trw.addDoubleCollection("drmom1", drmom1); + trw.addDoubleCollection("drmom2", drmom2); + + trw.addCollection("sigJet", jetCollectionSig); + trw.addCollection("sigJetSDBeta00Z01", jetCollectionSigSDBeta00Z01); + + trw.fillTree(); + + }//event loop + + Bar.Update(nEvent); + Bar.Print(); + Bar.PrintLine(); + + TTree *trOut = trw.getTree(); + + TFile *fout = new TFile(cmdline.value("-output", "JetToyHIResultGluonSplittingnalysis.root").c_str(), "RECREATE"); + trOut->Write(); + fout->Write(); + fout->Close(); + + double time_in_seconds = chrono::duration_cast + (chrono::steady_clock::now() - start_time).count() / 1000.0; + cout << "runFromFile: " << time_in_seconds << endl; +}