diff --git a/include/softDropGroomer.hh b/include/softDropGroomer.hh index 901dbc8..7723c90 100644 --- a/include/softDropGroomer.hh +++ b/include/softDropGroomer.hh @@ -12,6 +12,7 @@ #include "fastjet/contrib/SoftDrop.hh" #include "jetCollection.hh" +#include "jewelMatcher.hh" //--------------------------------------------------------------- // Description @@ -51,6 +52,10 @@ public : std::vector> getConstituents() {return constituents_;} std::vector> getConstituents1() {return constituents1_;} std::vector> getConstituents2() {return constituents2_;} + + std::vector doGroomingWithJewelSub(jetCollection &c, std::vector particlesDummy); + std::vector doGroomingWithJewelSub(std::vector v, std::vector particlesDummy); + std::vector doGroomingWithJewelSub(std::vector particlesDummy); }; softDropGroomer::softDropGroomer(double zcut, double beta, double r0) @@ -188,4 +193,107 @@ std::vector softDropGroomer::doGrooming() return fjOutputs_; } +std::vector softDropGroomer::doGroomingWithJewelSub(jetCollection &c, std::vector particlesDummy) +{ + return doGroomingWithJewelSub(c.getJet(),particlesDummy); +} + +std::vector softDropGroomer::doGroomingWithJewelSub(std::vector v, std::vector particlesDummy) +{ + setInputJets(v); + return doGroomingWithJewelSub(particlesDummy); +} + +std::vector softDropGroomer::doGroomingWithJewelSub(std::vector particlesDummy) +{ + fjOutputs_.reserve(fjInputs_.size()); + zg_.reserve(fjInputs_.size()); + drBranches_.reserve(fjInputs_.size()); + dr12_.reserve(fjInputs_.size()); + constituents_.reserve(fjInputs_.size()); + constituents1_.reserve(fjInputs_.size()); + constituents2_.reserve(fjInputs_.size()); + + for(fastjet::PseudoJet& jet : fjInputs_) { + std::vector particles, ghosts; + fastjet::SelectorIsPureGhost().sift(jet.constituents(), ghosts, particles); + + fastjet::JetDefinition jet_def(fastjet::cambridge_algorithm, fastjet::JetDefinition::max_allowable_R); + fastjet::ClusterSequence cs(particles, jet_def); + + //std::cout << "reclustered with CA" << std::endl; + std::vector tempJets = fastjet::sorted_by_pt(cs.inclusive_jets()); + if(tempJets.size()<1) { + fjOutputs_.push_back(fastjet::PseudoJet(0.,0.,0.,0.)); + zg_.push_back(-1.); + drBranches_.push_back(-1.); + constituents1_.push_back(std::vector()); + constituents2_.push_back(std::vector()); + continue; + } + + fastjet::PseudoJet CurrentJet = tempJets[0]; + fastjet::PseudoJet Part1, Part2; + fastjet::PseudoJet sj1, sj2; + double zg = -1.; + int ndrop = 0; + + // std::cout << "start grooming procedure" << std::endl; + while(CurrentJet.has_parents(Part1, Part2)) { + + if (CurrentJet.pt2() <= 0) break; + + zg = -1.; + + double deltaRsq = Part1.squared_distance(Part2); + double cut = zcut_ * std::pow(deltaRsq / r0_*r0_, 0.5*beta_); + + sj1 = GetCorrectedJet(Part1,particlesDummy); + sj2 = GetCorrectedJet(Part2,particlesDummy); + + if(sj1.pt() + sj2.pt() > 0 && sj1.E()>0. && sj2.E()>0. && sj1.m()>0. && sj2.m()>0.) + zg = min(sj1.pt(), sj2.pt()) / (sj1.pt() + sj2.pt()); + + if(zg sj2.pt()) + CurrentJet = Part1; + else + CurrentJet = Part2; + zg = -1.; + ++ndrop; + } else { + break; + } + } + + //build groomed jet + fastjet::PseudoJet transformedJet; + if(zg>0.) + transformedJet = fastjet::PseudoJet(sj1.px()+sj2.px(),sj1.py()+sj2.py(),sj1.pz()+sj2.pz(),sj1.E()+sj2.E()); + + if ( transformedJet == 0 ) { + fjOutputs_.push_back(fastjet::PseudoJet(0.,0.,0.,0.)); + zg_.push_back(-1.); + drBranches_.push_back(-1.); + constituents1_.push_back(std::vector()); + constituents2_.push_back(std::vector()); + } else { + //storing unsubtracted constituents of the subjets + if(Part1.has_constituents ()) constituents1_.push_back(Part1.constituents()); + if(Part2.has_constituents ()) constituents2_.push_back(Part2.constituents()); + + //get distance between the two subjets + //double deltaR = std::sqrt((sj1.eta() - sj2.eta())*(sj1.eta() - sj2.eta()) + (sj1.delta_phi_to(sj2))*(sj2.delta_phi_to(sj1))); + double deltaR = sj1.delta_R(sj2); + + fjOutputs_.push_back( transformedJet ); //put CA reclusterd jet after softDrop into vector + zg_.push_back(zg); + dr12_.push_back(deltaR); + drBranches_.push_back(ndrop); + } + } //jet loop + return fjOutputs_; +} + + #endif diff --git a/runSDGenVariousJewelSub.cc b/runSDGenVariousJewelSub.cc new file mode 100644 index 0000000..35d2ceb --- /dev/null +++ b/runSDGenVariousJewelSub.cc @@ -0,0 +1,241 @@ +#include +#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/jetCollection.hh" +#include "include/csSubtractor.hh" +#include "include/csSubtractorFullEvent.hh" +#include "include/skSubtractor.hh" +#include "include/softDropGroomer.hh" +#include "include/treeWriter.hh" +#include "include/jetMatcher.hh" +#include "include/randomCones.hh" +#include "include/Angularity.hh" +#include "include/jewelMatcher.hh" + +using namespace std; +using namespace fastjet; + +//./runSDGenVariousJewelSub -hard /eos/project/j/jetquenching/JetWorkshop2017/samples/jewel/DiJet/RecoilOn_0_10/Jewel_0_T_0.pu14 -nev 10 + +int main (int argc, char ** argv) { + + auto start_time = std::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"); + + std::cout << "will run on " << nEvent << " events" << std::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; + fastjet::GhostedAreaSpec ghost_spec(ghostRapMax, active_area_repeats, ghost_area); + fastjet::AreaDefinition area_def = fastjet::AreaDefinition(fastjet::active_area,ghost_spec); + fastjet::JetDefinition jet_def(antikt_algorithm, R); + + double jetRapMax = 3.0; + fastjet::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); + + std::vector particlesMerged = mixer.particles(); + + std::vector eventWeight; + eventWeight.push_back(mixer.hard_weight()); + eventWeight.push_back(mixer.pu_weight()); + + // cluster hard event only + vector particlesDummy, particlesReal; + vector particlesBkg, particlesSig; + SelectorVertexNumber(-1).sift(particlesMerged, particlesDummy, particlesReal); + SelectorVertexNumber(0).sift(particlesReal, particlesSig, particlesBkg); + + for(int i = 0; i < (int)particlesDummy.size(); i++) { + if(particlesDummy[i].perp() < 1e-5 && fabs(particlesDummy[i].pz()) > 2000) { + particlesDummy.erase(particlesDummy.begin() + i); + i = i - 1; + } + } + + //--------------------------------------------------------------------------- + // jet clustering + //--------------------------------------------------------------------------- + + fastjet::ClusterSequenceArea csSig(particlesSig, jet_def, area_def); + jetCollection jetCollectionSig(sorted_by_pt(jet_selector(csSig.inclusive_jets(10.)))); + jetCollection jetCollectionSigJewel(GetCorrectedJets(jetCollectionSig.getJet(), particlesDummy)); + + //--------------------------------------------------------------------------- + // 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()); + + softDropGroomer sdgSigSubBeta00Z01(0.1, 0.0, R); + jetCollection jetCollectionSigSDSubBeta00Z01(sdgSigSubBeta00Z01.doGroomingWithJewelSub(jetCollectionSig,particlesDummy)); + jetCollectionSigSDSubBeta00Z01.addVector("zgSigSDSubBeta00Z01", sdgSigSubBeta00Z01.getZgs()); + jetCollectionSigSDSubBeta00Z01.addVector("ndropSigSDSubBeta00Z01", sdgSigSubBeta00Z01.getNDroppedSubjets()); + jetCollectionSigSDSubBeta00Z01.addVector("dr12SigSDSubBeta00Z01", sdgSigSubBeta00Z01.getDR12()); + + softDropGroomer sdgSigBeta00Z02(0.2, 0.0, R); + jetCollection jetCollectionSigSDBeta00Z02(sdgSigBeta00Z02.doGrooming(jetCollectionSig)); + jetCollectionSigSDBeta00Z02.addVector("zgSigSDBeta00Z02", sdgSigBeta00Z02.getZgs()); + jetCollectionSigSDBeta00Z02.addVector("ndropSigSDBeta00Z02", sdgSigBeta00Z02.getNDroppedSubjets()); + jetCollectionSigSDBeta00Z02.addVector("dr12SigSDBeta00Z02", sdgSigBeta00Z02.getDR12()); + + softDropGroomer sdgSigSubBeta00Z02(0.2, 0.0, R); + jetCollection jetCollectionSigSDSubBeta00Z02(sdgSigSubBeta00Z02.doGroomingWithJewelSub(jetCollectionSig,particlesDummy)); + jetCollectionSigSDSubBeta00Z02.addVector("zgSigSDSubBeta00Z02", sdgSigSubBeta00Z02.getZgs()); + jetCollectionSigSDSubBeta00Z02.addVector("ndropSigSDSubBeta00Z02", sdgSigSubBeta00Z02.getNDroppedSubjets()); + jetCollectionSigSDSubBeta00Z02.addVector("dr12SigSDSubBeta00Z02", sdgSigSubBeta00Z02.getDR12()); + + softDropGroomer sdgSigBeta15Z05(0.5, 1.5, R); + jetCollection jetCollectionSigSDBeta15Z05(sdgSigBeta15Z05.doGrooming(jetCollectionSig)); + jetCollectionSigSDBeta15Z05.addVector("zgSigSDBeta15Z05", sdgSigBeta15Z05.getZgs()); + jetCollectionSigSDBeta15Z05.addVector("ndropSigSDBeta15Z05", sdgSigBeta15Z05.getNDroppedSubjets()); + jetCollectionSigSDBeta15Z05.addVector("dr12SigSDBeta15Z05", sdgSigBeta15Z05.getDR12()); + + softDropGroomer sdgSigSubBeta15Z05(0.5, 1.5, R); + jetCollection jetCollectionSigSDSubBeta15Z05(sdgSigSubBeta15Z05.doGroomingWithJewelSub(jetCollectionSig,particlesDummy)); + jetCollectionSigSDSubBeta15Z05.addVector("zgSigSDSubBeta15Z05", sdgSigSubBeta15Z05.getZgs()); + jetCollectionSigSDSubBeta15Z05.addVector("ndropSigSDSubBeta15Z05", sdgSigSubBeta15Z05.getNDroppedSubjets()); + jetCollectionSigSDSubBeta15Z05.addVector("dr12SigSDSubBeta15Z05", sdgSigSubBeta15Z05.getDR12()); + + softDropGroomer sdgSigBetam1Z01(0.1, -1., R); + jetCollection jetCollectionSigSDBetam1Z01(sdgSigBetam1Z01.doGrooming(jetCollectionSig)); + jetCollectionSigSDBetam1Z01.addVector("zgSigSDBetam1Z01", sdgSigBetam1Z01.getZgs()); + jetCollectionSigSDBetam1Z01.addVector("ndropSigSDBetam1Z01", sdgSigBetam1Z01.getNDroppedSubjets()); + jetCollectionSigSDBetam1Z01.addVector("dr12SigSDBetam1Z01", sdgSigBetam1Z01.getDR12()); + + softDropGroomer sdgSigSubBetam1Z01(0.1, -1., R); + jetCollection jetCollectionSigSDSubBetam1Z01(sdgSigSubBetam1Z01.doGroomingWithJewelSub(jetCollectionSig,particlesDummy)); + jetCollectionSigSDSubBetam1Z01.addVector("zgSigSDSubBetam1Z01", sdgSigSubBetam1Z01.getZgs()); + jetCollectionSigSDSubBetam1Z01.addVector("ndropSigSDSubBetam1Z01", sdgSigSubBetam1Z01.getNDroppedSubjets()); + jetCollectionSigSDSubBetam1Z01.addVector("dr12SigSDSubBetam1Z01", sdgSigSubBetam1Z01.getDR12()); + + softDropGroomer sdgSigBetam1Z02(0.2, -1., R); + jetCollection jetCollectionSigSDBetam1Z02(sdgSigBetam1Z02.doGrooming(jetCollectionSig)); + jetCollectionSigSDBetam1Z02.addVector("zgSigSDBetam1Z02", sdgSigBetam1Z02.getZgs()); + jetCollectionSigSDBetam1Z02.addVector("ndropSigSDBetam1Z02", sdgSigBetam1Z02.getNDroppedSubjets()); + jetCollectionSigSDBetam1Z02.addVector("dr12SigSDBetam1Z02", sdgSigBetam1Z02.getDR12()); + + softDropGroomer sdgSigSubBetam1Z02(0.2, -1.0, R); + jetCollection jetCollectionSigSDSubBetam1Z02(sdgSigSubBetam1Z02.doGroomingWithJewelSub(jetCollectionSig,particlesDummy)); + jetCollectionSigSDSubBetam1Z02.addVector("zgSigSDSubBetam1Z02", sdgSigSubBetam1Z02.getZgs()); + jetCollectionSigSDSubBetam1Z02.addVector("ndropSigSDSubBetam1Z02", sdgSigSubBetam1Z02.getNDroppedSubjets()); + jetCollectionSigSDSubBetam1Z02.addVector("dr12SigSDSubBetam1Z02", sdgSigSubBetam1Z02.getDR12()); + + softDropGroomer sdgSigBetam2Z01(0.1, -2., R); + jetCollection jetCollectionSigSDBetam2Z01(sdgSigBetam2Z01.doGrooming(jetCollectionSig)); + jetCollectionSigSDBetam2Z01.addVector("zgSigSDBetam2Z01", sdgSigBetam2Z01.getZgs()); + jetCollectionSigSDBetam2Z01.addVector("ndropSigSDBetam2Z01", sdgSigBetam2Z01.getNDroppedSubjets()); + jetCollectionSigSDBetam2Z01.addVector("dr12SigSDBetam2Z01", sdgSigBetam2Z01.getDR12()); + + softDropGroomer sdgSigSubBetam2Z01(0.1, -2.0, R); + jetCollection jetCollectionSigSDSubBetam2Z01(sdgSigSubBetam2Z01.doGroomingWithJewelSub(jetCollectionSig,particlesDummy)); + jetCollectionSigSDSubBetam2Z01.addVector("zgSigSDSubBetam2Z01", sdgSigSubBetam2Z01.getZgs()); + jetCollectionSigSDSubBetam2Z01.addVector("ndropSigSDSubBetam2Z01", sdgSigSubBetam2Z01.getNDroppedSubjets()); + jetCollectionSigSDSubBetam2Z01.addVector("dr12SigSDSubBetam2Z01", sdgSigSubBetam2Z01.getDR12()); + + softDropGroomer sdgSigBetam2Z005(0.05, -2., R); + jetCollection jetCollectionSigSDBetam2Z005(sdgSigBetam2Z005.doGrooming(jetCollectionSig)); + jetCollectionSigSDBetam2Z005.addVector("zgSigSDBetam2Z005", sdgSigBetam2Z005.getZgs()); + jetCollectionSigSDBetam2Z005.addVector("ndropSigSDBetam2Z005", sdgSigBetam2Z005.getNDroppedSubjets()); + jetCollectionSigSDBetam2Z005.addVector("dr12SigSDBetam2Z005", sdgSigBetam2Z005.getDR12()); + + softDropGroomer sdgSigSubBetam2Z005(0.05, -2.0, R); + jetCollection jetCollectionSigSDSubBetam2Z005(sdgSigSubBetam2Z005.doGroomingWithJewelSub(jetCollectionSig,particlesDummy)); + jetCollectionSigSDSubBetam2Z005.addVector("zgSigSDSubBetam2Z005", sdgSigSubBetam2Z005.getZgs()); + jetCollectionSigSDSubBetam2Z005.addVector("ndropSigSDSubBetam2Z005", sdgSigSubBetam2Z005.getNDroppedSubjets()); + jetCollectionSigSDSubBetam2Z005.addVector("dr12SigSDSubBetam2Z005", sdgSigSubBetam2Z005.getDR12()); + + //--------------------------------------------------------------------------- + // write tree + //--------------------------------------------------------------------------- + + //Give variable we want to write out to treeWriter. + //Only vectors of the types 'jetCollection', and 'double', 'int', 'fastjet::PseudoJet' are supported + + trw.addCollection("eventWeight", eventWeight); + + trw.addCollection("sigJet", jetCollectionSig); + trw.addCollection("sigJetSDBeta00Z01", jetCollectionSigSDBeta00Z01); + trw.addCollection("sigJetSDBeta00Z02", jetCollectionSigSDBeta00Z02); + trw.addCollection("sigJetSDBeta15Z05", jetCollectionSigSDBeta15Z05); + trw.addCollection("sigJetSDBetam1Z01", jetCollectionSigSDBetam1Z01); + trw.addCollection("sigJetSDBetam1Z02", jetCollectionSigSDBetam1Z02); + trw.addCollection("sigJetSDBetam2Z01", jetCollectionSigSDBetam2Z01); + trw.addCollection("sigJetSDBetam2Z005", jetCollectionSigSDBetam2Z005); + + trw.addCollection("sigJetSub", jetCollectionSigJewel); + trw.addCollection("sigJetSDBeta00Z01Sub", jetCollectionSigSDSubBeta00Z01); + trw.addCollection("sigJetSDBeta00Z02Sub", jetCollectionSigSDSubBeta00Z02); + trw.addCollection("sigJetSDBeta15Z05Sub", jetCollectionSigSDSubBeta15Z05); + trw.addCollection("sigJetSDBetam1Z01Sub", jetCollectionSigSDSubBetam1Z01); + trw.addCollection("sigJetSDBetam1Z02Sub", jetCollectionSigSDSubBetam1Z02); + trw.addCollection("sigJetSDBetam2Z01Sub", jetCollectionSigSDSubBetam2Z01); + trw.addCollection("sigJetSDBetam2Z005Sub", jetCollectionSigSDSubBetam2Z005); + + trw.fillTree(); + + }//event loop + + Bar.Update(nEvent); + Bar.Print(); + Bar.PrintLine(); + + TTree *trOut = trw.getTree(); + + TFile *fout = new TFile("JetToyHIResultSDGenSub.root","RECREATE"); + trOut->Write(); + fout->Write(); + fout->Close(); + + double time_in_seconds = std::chrono::duration_cast + (std::chrono::steady_clock::now() - start_time).count() / 1000.0; + std::cout << "runFromFile: " << time_in_seconds << std::endl; +}