-
Notifications
You must be signed in to change notification settings - Fork 24
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Showing
5 changed files
with
250 additions
and
14 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1 +1 @@ | ||
| /soft/pythia8235 | ||
| /Users/mverweij/soft/pythia8235 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,208 @@ | ||
| #include <iostream> | ||
| #include <chrono> | ||
|
|
||
| #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" | ||
| #include "include/csSubtractor.hh" | ||
|
|
||
| using namespace std; | ||
| using namespace fastjet; | ||
|
|
||
| // ./runConstituentSubtraction -hard samples/PythiaEventsTune14PtHat120.pu14 -pileup samples/ThermalEventsMult7000PtAv1.20_1.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<int>("-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<PseudoJet> particlesMergedAll = mixer.particles(); | ||
|
|
||
| vector<double> 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<PseudoJet> partons = parton_selector(particlesMergedAll); | ||
|
|
||
| // select final state particles from hard event only | ||
| //vector<PseudoJet> particlesBkg, particlesSig; | ||
| //SelectorIsHard().sift(particlesMerged, particlesSig, particlesBkg); // this sifts the full event into two vectors of PseudoJet, one for the hard event, one for the underlying event | ||
|
|
||
| fastjet::Selector sig_selector = SelectorVertexNumber(0); | ||
| vector<PseudoJet> particlesSig = sig_selector(particlesMergedAll); | ||
|
|
||
| fastjet::Selector bkg_selector = SelectorVertexNumber(1); | ||
| vector<PseudoJet> particlesBkg = bkg_selector(particlesMergedAll); | ||
|
|
||
| vector<PseudoJet> particlesMerged = particlesBkg; | ||
| particlesMerged.insert( particlesMerged.end(), particlesSig.begin(), particlesSig.end() ); | ||
|
|
||
| std::cout << "#merged: " << particlesMerged.size() << " signal: " << particlesSig.size() << " bkg: " << particlesBkg.size() << std::endl; | ||
|
|
||
| //--------------------------------------------------------------------------- | ||
| // jet clustering of signal jets | ||
| //--------------------------------------------------------------------------- | ||
|
|
||
| fastjet::ClusterSequenceArea csSig(particlesSig, jet_def, area_def); | ||
| jetCollection jetCollectionSig(sorted_by_pt(jet_selector(csSig.inclusive_jets(10.)))); | ||
|
|
||
| //calculate some angularities | ||
| vector<double> widthSig; widthSig.reserve(jetCollectionSig.getJet().size()); | ||
| vector<double> 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); | ||
|
|
||
|
|
||
| //--------------------------------------------------------------------------- | ||
| // jet clustering of signal+background jets | ||
| //--------------------------------------------------------------------------- | ||
|
|
||
| fastjet::ClusterSequenceArea csRaw(particlesMerged, jet_def, area_def); | ||
| jetCollection jetCollectionRaw(sorted_by_pt(jet_selector(csRaw.inclusive_jets(10.)))); | ||
|
|
||
|
|
||
| //--------------------------------------------------------------------------- | ||
| // background subtraction | ||
| //--------------------------------------------------------------------------- | ||
|
|
||
| //run jet-by-jet constituent subtraction on mixed (hard+UE) event | ||
| csSubtractor csSub(R, 0., -1, 0.005,ghostRapMax,jetRapMax); | ||
| csSub.setInputParticles(particlesMerged); | ||
| jetCollection jetCollectionCS(csSub.doSubtraction()); | ||
|
|
||
| //Background densities used by constituent subtraction | ||
| std::vector<double> rho; | ||
| std::vector<double> rhom; | ||
| rho.push_back(csSub.getRho()); | ||
| rhom.push_back(csSub.getRhoM()); | ||
|
|
||
| //match CS jets to signal jets | ||
| jetMatcher jmCS(R); | ||
| jmCS.setBaseJets(jetCollectionCS); | ||
| jmCS.setTagJets(jetCollectionSig); | ||
| jmCS.matchJets(); | ||
|
|
||
| jmCS.reorderedToTag(jetCollectionCS); | ||
|
|
||
| //match Raw(=unsubtracted) jets to signal jets | ||
| jetMatcher jmRaw(R); | ||
| jmRaw.setBaseJets(jetCollectionRaw); | ||
| jmRaw.setTagJets(jetCollectionSig); | ||
| jmRaw.matchJets(); | ||
|
|
||
| jmRaw.reorderedToTag(jetCollectionRaw); | ||
|
|
||
| //--------------------------------------------------------------------------- | ||
| // 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.addCollection("csRho", rho); | ||
| trw.addCollection("csRhom", rhom); | ||
|
|
||
| trw.addPartonCollection("partons", partons); | ||
|
|
||
| trw.addCollection("sigJet", jetCollectionSig); | ||
| trw.addCollection("sigJetSDBeta00Z01", jetCollectionSigSDBeta00Z01); | ||
|
|
||
| trw.addCollection("csJet", jetCollectionCS); | ||
| trw.addCollection("rawJet", jetCollectionRaw); | ||
|
|
||
| trw.fillTree(); | ||
|
|
||
| }//event loop | ||
|
|
||
| Bar.Update(nEvent); | ||
| Bar.Print(); | ||
| Bar.PrintLine(); | ||
|
|
||
| TTree *trOut = trw.getTree(); | ||
|
|
||
| TFile *fout = new TFile(cmdline.value<string>("-output", "JetToyHIResultConstituentSubtraction.root").c_str(), "RECREATE"); | ||
| trOut->Write(); | ||
| fout->Write(); | ||
| fout->Close(); | ||
|
|
||
| double time_in_seconds = chrono::duration_cast<chrono::milliseconds> | ||
| (chrono::steady_clock::now() - start_time).count() / 1000.0; | ||
| cout << "runFromFile: " << time_in_seconds << endl; | ||
| } |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters