|
[Rivet] missingEt variable in WFinder constructorInês Ochoa Ines.Ochoa at cern.chMon Feb 2 12:30:15 GMT 2015
Hi Andy, I attached a simpler version of the routine I'm using. I run it with Rivet 2.1.1 and I see that: 1. the nuPt variable I print out has no cut on it 2. changing the value of the MET argument does not change the number of events with reconstructed W's above a given pT Thanks! Inês On 02/02/15 12:02, Andy Buckley wrote: > Hi again Ines, > > Can you send your example analysis code so we can easily re-run & check? > Which version of Rivet are you using? I should check that some of the > recent interface changes did not accidentally break a MET cut, but this > will be far easier with an example code that I can just run out of the > box in the latest Rivet. > > Cheers, > Andy > > > On 27/01/15 11:06, Inês Ochoa wrote: >> Hi Frank, >> >> The WFinder is defined as: >> >> WFinder wfinder(-6.0, 6.0, 20.0*GeV, PID::MUON, 0.0*GeV, 1000.0*GeV, >> 20.0*GeV, 0.1); >> >> I can see that the documentation on WFinder now has some modifications >> on the list of arguments, but this is consistent with previous >> versions... I've also checked that the constituentLeptons()[0] pT does >> have a cut at 20 GeV. >> >> Cheers, >> Inês >> >> On 27/01/15 10:56, Frank Siegert wrote: >>> Hi Inês, >>> >>> The missingET argument should apply a pT cut to the (in)visible >>> momentum in your inputfs (which should include the neutrino). So your >>> observation would definitely not be expected, and I don't currently >>> see a reason why this would go wrong. Could you post a code snippet of >>> how your WFinder is defined in your analysis, or attach the full >>> analysis code? >>> >>> Cheers, >>> Frank >>> >>> >>> >>> >>> On 27 January 2015 at 11:44, Inês Ochoa <Ines.Ochoa at cern.ch> wrote: >>>> Dear Rivet developers, >>>> >>>> I am using Rivet 2.1.1 and selecting a W candidate with WFinder. I am >>>> applying a cut on the missing Et (as one of the arguments of >>>> WFinder()) but >>>> when accessing the neutrino with constituentNeutrinos()[0] I see that >>>> there >>>> is no cut on its Pt. In fact, changing the MET cut on the WFinder >>>> constructor seems to have no impact on the resulting number of events. >>>> >>>> Is this a known problem? Should that argument be ignored and the cut >>>> applied >>>> on the constituentNeutrino object instead? >>>> >>>> Thanks! >>>> Inês >>>> _______________________________________________ >>>> Rivet mailing list >>>> Rivet at projects.hepforge.org >>>> https://www.hepforge.org/lists/listinfo/rivet >> _______________________________________________ >> Rivet mailing list >> Rivet at projects.hepforge.org >> https://www.hepforge.org/lists/listinfo/rivet > > -------------- next part -------------- #include "Rivet/Tools/ParticleIdUtils.hh" #include "Rivet/ParticleName.hh" #include "Rivet/Analysis.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/WFinder.hh" #include "YODA/WriterAIDA.h" #include "Rivet/AnalysisLoader.hh" #include "fastjet/tools/MassDropTagger.hh" #include "fastjet/tools/Filter.hh" #include <TH1F.h> #include <TH2F.h> #include <TFile.h> #include <TTree.h> #include <TBranch.h> namespace Rivet { class exampleAnalysis : public Analysis { public: /// @name Constructors etc. //@{ /// Constructor exampleAnalysis() : Analysis("exampleAnalysis"){} //@} public: /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { // WFinder: etaRanges,pTmin,pid,m2_min,m2_max,missingET,dRmax WFinder wfinder(-6.0, 6.0, 20.0*GeV, PID::MUON, 0.0*GeV, 1000.0*GeV, 20.0*GeV, 0.1); //original addProjection(wfinder, "WFinder"); VetoedFinalState veto; veto.addVetoOnThisFinalState(wfinder); const FastJets ca(veto, FastJets::CAM, 1.2); const FastJets akt(veto, FastJets::ANTIKT, 0.4); //ca.useInvisibles(true); //akt.useInvisibles(true); addProjection(ca, "CA12"); addProjection(akt, "AKT04"); setNeedsCrossSection(true); //TREE STUFF _treeFileName = "rivetTree.root"; _treeFile = new TFile(_treeFileName, "recreate"); _myTree = new TTree("myTree", "myTree"); //control histograms _h_weights = new TH1F("_h_weights", "MC weights", 10, -5.0, 5.0); _h_xsec = new TH1F("_h_xsec", "x-sec check", 1, 0, 1.0); _h_PtW = new TH1F("_h_PtW", "Pt(W) > 200 GeV", 200, 0.0, 2000.0); _h_precut_PtW = new TH1F("_h_precut_PtW", "Pt(W) precut", 200, 0.0, 2000.0); _h_EtaW = new TH1F("_h_EtaW", "EtaW, Pt(W) > 200 GeV", 1000, -7.0, 7.0); _h_precut_EtaW = new TH1F("_h_precut_EtaW", "EtaW, precut", 1000, -7.0, 7.0); _h_nW = new TH1F("_h_nW", "nW", 20, 0.0, 20.0); //general tree variables _myTree->Branch("weight", &var_weight, "weight/D"); // _myTree->Branch("eventNumber", &var_eventNumber, "eventNumber/I"); _myTree->Branch("WpT", &var_WpT, "WpT/D"); _myTree->Branch("WEta", &var_WEta, "WEta/D"); _myTree->Branch("leptonPt", &var_leptonPt, "leptonPt/D"); _myTree->Branch("leptonEta", &var_leptonEta, "leptonEta/D"); _myTree->Branch("nuPt", &var_nuPt, "nuPt/D"); _myTree->Branch("nuEta", &var_nuEta, "nuEta/D"); _myTree->Branch("alphaS", &var_alphaS, "alphaS/D"); _myTree->Branch("id1", &var_id1, "id1/I"); _myTree->Branch("id2", &var_id2, "id2/I"); _myTree->Branch("x1", &var_x1, "x1/D"); _myTree->Branch("x2", &var_x2, "x2/D"); _myTree->Branch("Q", &var_Q, "Q/D"); //anti-kT variables _myTree->Branch("JetPt1", &var_JetPt1, "JetPt1/D"); _myTree->Branch("JetPt2", &var_JetPt2, "JetPt2/D"); _myTree->Branch("JetPt3", &var_JetPt3, "JetPt3/D"); _myTree->Branch("JetEta1", &var_JetEta1, "JetEta1/D"); _myTree->Branch("JetEta2", &var_JetEta2, "JetEta2/D"); _myTree->Branch("JetEta3", &var_JetEta3, "JetEta3/D"); _myTree->Branch("JetPhi1", &var_JetPhi1, "JetPhi1/D"); _myTree->Branch("JetPhi2", &var_JetPhi2, "JetPhi2/D"); _myTree->Branch("JetPhi3", &var_JetPhi3, "JetPhi3/D"); //CA split/filt variables _myTree->Branch("FatJetPt", &var_FatJetPt, "FatJetPt/D"); _myTree->Branch("FatJetEta", &var_FatJetEta, "FatJetEta/D"); _myTree->Branch("FatJetRap", &var_FatJetRap, "FatJetRap/D"); _myTree->Branch("FatJetPhi", &var_FatJetPhi, "FatJetPhi/D"); _myTree->Branch("FatJet2Pt", &var_FatJet2Pt, "FatJet2Pt/D"); _myTree->Branch("FatJet2Eta", &var_FatJet2Eta, "FatJet2Eta/D"); _myTree->Branch("FatJet2Rap", &var_FatJet2Rap, "FatJet2Rap/D"); _myTree->Branch("FatJet2Phi", &var_FatJet2Phi, "FatJet2Phi/D"); _myTree->Branch("Rfilt", &var_Rfilt, "Rfilt/D"); _myTree->Branch("nSj", &var_nSj, "nSj/I"); _myTree->Branch("SubJetPt1", &var_SubJetPt1, "SubJetPt1/D"); _myTree->Branch("SubJetEta1", &var_SubJetEta1, "SubJetEta1/D"); _myTree->Branch("SubJetPhi1", &var_SubJetPhi1, "SubJetPhi1/D"); _myTree->Branch("SubJetPt2", &var_SubJetPt2, "SubJetPt2/D"); _myTree->Branch("SubJetEta2", &var_SubJetEta2, "SubJetEta2/D"); _myTree->Branch("SubJetPhi2", &var_SubJetPhi2, "SubJetPhi2/D"); _myTree->Branch("SubJetPt3", &var_SubJetPt3, "SubJetPt3/D"); _myTree->Branch("SubJetEta3", &var_SubJetEta3, "SubJetEta3/D"); _myTree->Branch("SubJetPhi3", &var_SubJetPhi3, "SubJetPhi3/D"); //event counters nev = 0; sumw = 0; //event counters for cross-checks - number of events passing each selection //njs = 0; //njs_alt = 0; //nakt = 0; } /// Perform the per-event analysis void analyze(const Event& event) { double weight = event.weight(); var_weight = weight; //var_eventNumber = event.genEvent()->event_number(); _h_weights->Fill(weight); _h_xsec->Fill(0.0,weight); using namespace fastjet; sumw = sumw + weight; //----W stuff----// const WFinder& wfinder = applyProjection<WFinder>(event, "WFinder"); if ( wfinder.bosons().size() != 1 ) vetoEvent; FourMomentum wmom(0.0,0.0,0.0,0.0); wmom = wfinder.bosons().front().momentum(); float ptW = wmom.pT(); float EtaW = wmom.eta(); _h_precut_PtW->Fill(ptW/GeV, weight); _h_precut_EtaW->Fill(EtaW, weight); _h_nW->Fill(wfinder.bosons().size(), weight); //pT>200 GeV cut if ( ptW < 200*GeV ) vetoEvent; _h_EtaW->Fill(EtaW, weight); _h_PtW->Fill(ptW/GeV, weight); nev = nev + weight; //look at lepton const FourMomentum& l = wfinder.constituentLeptons()[0].momentum(); var_leptonPt = l.pT(); var_leptonEta = l.eta(); //look at MET const FourMomentum& nu = wfinder.constituentNeutrinos()[0].momentum(); var_nuPt = nu.pT(); var_nuEta = nu.eta(); //----PDF stuff----// HepMC::PdfInfo pdfi = *(event.genEvent()->pdf_info()); var_Q = pdfi.scalePDF(); var_x1 = pdfi.x1(); var_x2 = pdfi.x2(); var_id1 = pdfi.id1(); var_id2 = pdfi.id2(); var_alphaS = event.genEvent()->alphaQCD(); //----Jet Stuff---// const FastJets& ca = applyProjection<FastJets>(event, "CA12"); const FastJets& akt = applyProjection<FastJets>(event, "AKT04"); std::vector<Jet> akt_signal; std::vector<Jet> ca_total; std::vector<PseudoJet> ca_signal, casj_signal; std::vector<PseudoJet> ca_alt_signal, casj_alt_signal; foreach (const Jet& ajet, akt.jetsByPt(20*GeV)) { akt_signal.push_back(ajet); } foreach (const Jet& cajet, ca.jetsByPt(150*GeV)) { ca_total.push_back(cajet); } //Define the mass drop tagger to be used MassDropTagger mdtagger(0.67, 0.09); //default values // Split and filtering, place loser cut on CA jets before grooming // then take split/filtered ones only if pT>150 GeV PseudoJets cajets = ca.pseudoJetsByPt(20*GeV); //20 is a bit arbitrary here foreach (const PseudoJet pjet, cajets) { //MD PseudoJet mdjet = mdtagger(pjet); if ( mdjet == 0 ) continue; var_mu = mdjet.structure_of<MassDropTagger>().mu(); var_y = mdjet.structure_of<MassDropTagger>().y(); //Filtering vector<PseudoJet> mdjet_pieces = mdjet.pieces(); var_Rfilt = min(0.3, 0.5 * mdjet_pieces[0].delta_R(mdjet_pieces[1])); PseudoJet filt_jet = Filter(var_Rfilt, SelectorNHardest(3))(mdjet); if ( filt_jet.perp() < 150*GeV ) continue; //Subjet collection vector<PseudoJet> subjets = filt_jet.pieces(); var_nSj = subjets.size(); ca_signal.push_back(filt_jet); foreach ( const PseudoJet sj, subjets) { casj_signal.push_back(sj); } } //start analysis cuts here - consider all possibilities and apply specific cuts in the plot routine if ( akt_signal.size() > 0 || ca_signal.size() > 0 || ca_alt_signal.size() > 0 ) { // nakt = nakt + weight; // njs = njs + weight; // njs_alt = njs_alt + weight; var_WpT = ptW; var_WEta = EtaW; foreach ( const Jet& akt, akt_signal ) { var_nAktJets++; //no eta cut } foreach ( const PseudoJet ca, ca_signal ) { var_nCAJets++; //no eta cut } var_HTakt = wmom.perp(); var_HTca = wmom.perp(); var_AltHTca = wmom.perp(); var_nBTagJets = var_nBMatchJets = var_nBMatchSubJets = var_AltnBMatchSubJets = 0; //------akt jets-----// if ( akt_signal.size() > 0 ) { var_JetPt1 = akt_signal[0].perp(); var_JetEta1 = akt_signal[0].eta(); var_JetPhi1 = akt_signal[0].phi(); var_JetB1 = bTag(akt_signal[0]); if ( akt_signal.size() > 1 ) { var_JetPt2 = akt_signal[1].perp(); var_JetEta2 = akt_signal[1].eta(); var_JetPhi2 = akt_signal[1].phi(); double dEtajj = deltaEta(var_JetEta1, var_JetEta2); double dPhijj = deltaPhi(var_JetPhi1, var_JetPhi2); var_dRjj = sqrt(dEtajj*dEtajj + dPhijj*dPhijj); FourMomentum HiggsCand_1(akt_signal[0].momentum()+akt_signal[1].momentum()); var_Mjj = HiggsCand_1.mass(); var_JetB2 = bTag(akt_signal[1]); var_JetB1alt = bMatch(event, akt_signal[0], akt_signal[1]).first; //new var_JetB2alt = bMatch(event, akt_signal[0], akt_signal[1]).second; //new var_JetC1alt = cMatch(event, akt_signal[0], akt_signal[1]).first; //new var_JetC2alt = cMatch(event, akt_signal[0], akt_signal[1]).second; //new var_dEtaWjj = deltaEta(HiggsCand_1.eta(), wmom.eta()); FourMomentum WH_1(HiggsCand_1 + wmom); var_MWjj = WH_1.mass(); } else { var_JetPt2 = -999.; var_JetEta2 = -999.; var_JetPhi2 = -999.; var_dRjj = -999.; var_Mjj = -999.; var_JetB2 = false; var_JetB2alt = false; var_JetB1alt = false; //new var_JetC2alt = false; var_JetC1alt = false; //new var_JetTau1alt = false; var_JetTau2alt = false; } if ( akt_signal.size() > 2 ) { var_JetPt3 = akt_signal[2].perp(); var_JetEta3 = akt_signal[2].eta(); var_JetPhi3 = akt_signal[2].phi(); } else { var_JetPt3 = -999.; var_JetEta3 = -999.; var_JetPhi3 = -999.; } } else { var_JetPt1 = -999.; var_JetEta1 = -999.; var_JetPhi1 = -999.; var_JetB1 = false; } //-----ca jets-----// if ( ca_signal.size() > 0 ) { var_FatJetPt = ca_signal[0].perp(); var_FatJetEta = ca_signal[0].eta(); var_FatJetRap = ca_signal[0].rap(); var_FatJetPhi = ca_signal[0].phi(); FourMomentum HiggsCand_3(ca_signal[0].E(), ca_signal[0].px(), ca_signal[0].py(), ca_signal[0].pz()); var_Mj = HiggsCand_3.mass(); if ( ca_signal.size() > 1 ) { var_FatJet2Pt = ca_signal[1].perp(); var_FatJet2Eta = ca_signal[1].eta(); var_FatJet2Rap = ca_signal[1].rap(); var_FatJet2Phi = ca_signal[1].phi(); FourMomentum test(ca_signal[1].E(), ca_signal[1].px(), ca_signal[1].py(), ca_signal[1].pz()); } else { var_FatJet2Pt = -999.; var_FatJet2Eta = -999.; var_FatJet2Rap = -999.; var_FatJet2Phi = -999.; } var_SubJetPt1 = casj_signal[0].perp(); var_SubJetEta1 = casj_signal[0].eta(); var_SubJetPhi1 = casj_signal[0].phi(); var_SubJetPt2 = casj_signal[1].perp(); var_SubJetEta2 = casj_signal[1].eta(); var_SubJetPhi2 = casj_signal[1].phi(); if ( var_nSj > 2 ) { var_SubJetPt3 = casj_signal[2].perp(); var_SubJetEta3 = casj_signal[2].eta(); var_SubJetPhi3 = casj_signal[2].phi(); } else { var_SubJetPt3 = -999.; var_SubJetEta3 = -999.; var_SubJetPhi3 = -999.; } var_SubJetB1 = bTagPseudoJet(event, casj_signal[0], casj_signal[1], 0.3).first; //new var_SubJetB2 = bTagPseudoJet(event, casj_signal[0], casj_signal[1], 0.3).second; //new var_dPhiWJ = deltaPhi(var_FatJetPhi, wmom.phi()); var_dEtaWJ = deltaEta(var_FatJetEta, wmom.eta()); FourMomentum WH_2(HiggsCand_3 + wmom); var_MWJ = WH_2.mass(); } else { var_FatJetPt = -999.; var_FatJetEta = -999.; var_FatJetRap = -999.; var_FatJetPhi = -999.; var_Mj = -999.; var_SubJetPt1 = -999.; var_SubJetEta1 = -999.; var_SubJetPhi1 = -999.; var_SubJetPt2 = -999.; var_SubJetEta2 = -999.; var_SubJetPhi2 = -999.; var_SubJetB1 = false; var_SubJetB2 = false; var_SubJetC1 = false; var_SubJetC2 = false; var_SubJetTau1 = false; var_SubJetTau2 = false; var_SubJetVarB1 = false; var_SubJetVarB2 = false; var_dRCentral = -999.; var_dPhiWJ = -999.; } _myTree->Fill(); } } //analyse bool bTag(const Jet& jet) { foreach (const Particle& p, jet.particles()) { const PdgId pid = p.pdgId(); //if (abs(pid) == PID::BQUARK) return true; //new if (PID::isHadron(pid) && PID::hasBottom(pid) && p.momentum().pT()>5.0*GeV) return true; HepMC::GenVertex* gv = p.genParticle()->production_vertex(); if (gv) { foreach (const GenParticle* pi, Rivet::particles(gv, HepMC::ancestors)) { const PdgId pid2 = pi->pdg_id(); if (PID::isHadron(pid2) && PID::hasBottom(pid2) && pi->momentum().perp()>5.0*GeV) return true; } } } return false; } pair<bool, bool> bMatch(const Event& event, const Jet& jet1, const Jet& jet2) { //new -> look at 2 jets ParticleVector bquarks; for ( GenEvent::particle_const_iterator p=event.genEvent()->particles_begin(); p!=event.genEvent()->particles_end(); ++p ) { const PdgId pid = (*p)->pdg_id(); //if ( abs(pid) == PID::BQUARK ) isb = true; //new if ( PID::isHadron(pid) && PID::hasBottom(pid) && (*p)->momentum().perp()>5.0*GeV ) bquarks.push_back(Particle(**p)); } if ( bquarks.size() == 0 ) return pair<bool, bool>(false, false); bool isb1 = false, isb2 = false; foreach ( const Particle& b, bquarks ) { double dEta1 = deltaEta(jet1.eta(),b.eta()); double dPhi1 = deltaPhi(jet1.phi(),b.phi()); double dEta2 = deltaEta(jet2.eta(),b.eta()); double dPhi2 = deltaPhi(jet2.phi(),b.phi()); double dRmatch1 = sqrt(dEta1*dEta1 + dPhi1*dPhi1); double dRmatch2 = sqrt(dEta2*dEta2 + dPhi2*dPhi2); //closest sjet method if ( dRmatch1 < 0.4 && dRmatch2 > 0.4 ) isb1 = true; else if ( dRmatch1 > 0.4 && dRmatch2 < 0.4 ) isb2 = true; else if ( dRmatch1 < 0.4 && dRmatch2 < 0.4 ) { if ( dRmatch1 < dRmatch2 ) { if ( !isb1 ) isb1 = true; else isb2 = true; } else if ( dRmatch2 < dRmatch1 ) { if ( !isb2 ) isb2 = true; else isb1 = true; } } } return make_pair(isb1, isb2); } pair<bool, bool> bTagPseudoJet(const Event& event, const PseudoJet& pjet1, PseudoJet& pjet2, double rMatch) { ParticleVector bquarks; for ( GenEvent::particle_const_iterator p=event.genEvent()->particles_begin(); p!=event.genEvent()->particles_end(); ++p ) { const PdgId pid = (*p)->pdg_id(); //if ( abs(pid) == PID::BQUARK ) isb = true; //new if ( PID::isHadron(pid) && PID::hasBottom(pid) && (*p)->momentum().perp()>5.0*GeV ) bquarks.push_back(Particle(**p)); } if ( bquarks.size() == 0 ) return pair<bool, bool>(false, false); bool isb1 = false, isb2 = false; foreach ( const Particle& b, bquarks ) { double dEta1 = deltaEta(pjet1.eta(),b.eta()); double dPhi1 = deltaPhi(pjet1.phi(),b.phi()); double dEta2 = deltaEta(pjet2.eta(),b.eta()); double dPhi2 = deltaPhi(pjet2.phi(),b.phi()); double dRmatch1 = sqrt(dEta1*dEta1 + dPhi1*dPhi1); double dRmatch2 = sqrt(dEta2*dEta2 + dPhi2*dPhi2); if ( dRmatch1 < rMatch && dRmatch2 > rMatch ) isb1 = true; else if ( dRmatch1 > rMatch && dRmatch2 < rMatch ) isb2 = true; else if ( dRmatch1 < rMatch && dRmatch2 < rMatch ) { if ( dRmatch1 < dRmatch2 ) { if ( !isb1 ) isb1 = true; else isb2 = true; } else if ( dRmatch2 < dRmatch1 ) { if ( !isb2 ) isb2 = true; else isb1 = true; } } } return make_pair(isb1, isb2); } pair<bool, bool> cMatch(const Event& event, const Jet& jet1, const Jet& jet2) { //new -> look at 2 jets ParticleVector cquarks; for ( GenEvent::particle_const_iterator p=event.genEvent()->particles_begin(); p!=event.genEvent()->particles_end(); ++p ) { const PdgId pid = (*p)->pdg_id(); if ( PID::isHadron(pid) && PID::hasCharm(pid) && !PID::hasBottom(pid) && (*p)->momentum().perp()>5.0*GeV ) cquarks.push_back(Particle(**p)); } if ( cquarks.size() == 0 ) return pair<bool, bool>(false, false); bool isc1 = false, isc2 = false; foreach ( const Particle& c, cquarks ) { double dEta1 = deltaEta(jet1.eta(),c.eta()); double dPhi1 = deltaPhi(jet1.phi(),c.phi()); double dEta2 = deltaEta(jet2.eta(),c.eta()); double dPhi2 = deltaPhi(jet2.phi(),c.phi()); double dRmatch1 = sqrt(dEta1*dEta1 + dPhi1*dPhi1); double dRmatch2 = sqrt(dEta2*dEta2 + dPhi2*dPhi2); if ( dRmatch1 < 0.4 && dRmatch2 > 0.4 ) isc1 = true; else if ( dRmatch1 > 0.4 && dRmatch2 < 0.4 ) isc2 = true; else if ( dRmatch1 < 0.4 && dRmatch2 < 0.4 ) { if ( dRmatch1 < dRmatch2 ) { if ( !isc1 ) isc1 = true; else isc2 = true; } else if ( dRmatch2 < dRmatch1 ) { if ( !isc2 ) isc2 = true; else isc1 = true; } } } return make_pair(isc1, isc2); } pair<bool, bool> cTagPseudoJet(const Event& event, const PseudoJet& pjet1, PseudoJet& pjet2, double rMatch) { ParticleVector cquarks; for ( GenEvent::particle_const_iterator p=event.genEvent()->particles_begin(); p!=event.genEvent()->particles_end(); ++p ) { const PdgId pid = (*p)->pdg_id(); if ( PID::isHadron(pid) && PID::hasCharm(pid) && !PID::hasBottom(pid) && (*p)->momentum().perp()>5.0*GeV ) cquarks.push_back(Particle(**p)); } if ( cquarks.size() == 0 ) return pair<bool, bool>(false, false); bool isc1 = false, isc2 = false; foreach ( const Particle& c, cquarks ) { double dEta1 = deltaEta(pjet1.eta(),c.eta()); double dPhi1 = deltaPhi(pjet1.phi(),c.phi()); double dEta2 = deltaEta(pjet2.eta(),c.eta()); double dPhi2 = deltaPhi(pjet2.phi(),c.phi()); double dRmatch1 = sqrt(dEta1*dEta1 + dPhi1*dPhi1); double dRmatch2 = sqrt(dEta2*dEta2 + dPhi2*dPhi2); if ( dRmatch1 < rMatch && dRmatch2 > rMatch ) isc1 = true; else if ( dRmatch1 > rMatch && dRmatch2 < rMatch ) isc2 = true; else if ( dRmatch1 < rMatch && dRmatch2 < rMatch ) { if ( dRmatch1 < dRmatch2 ) { if ( !isc1 ) isc1 = true; else isc2 = true; } else if ( dRmatch2 < dRmatch1 ) { if ( !isc2 ) isc2 = true; else isc1 = true; } } } return make_pair(isc1, isc2); } //taus!! pair<bool, bool> tauMatch(const Event& event, const Jet& jet1, const Jet& jet2) { //new -> look at 2 jets ParticleVector taulep; for ( GenEvent::particle_const_iterator p=event.genEvent()->particles_begin(); p!=event.genEvent()->particles_end(); ++p ) { const PdgId pid = (*p)->pdg_id(); if ( abs(pid) == 15 ) taulep.push_back(Particle(**p)); //applying no pt cut on taus... } if ( taulep.size() == 0 ) return pair<bool, bool>(false, false); bool istau1 = false, istau2 = false; foreach ( const Particle& t, taulep ) { double dEta1 = deltaEta(jet1.eta(),t.eta()); double dPhi1 = deltaPhi(jet1.phi(),t.phi()); double dEta2 = deltaEta(jet2.eta(),t.eta()); double dPhi2 = deltaPhi(jet2.phi(),t.phi()); double dRmatch1 = sqrt(dEta1*dEta1 + dPhi1*dPhi1); double dRmatch2 = sqrt(dEta2*dEta2 + dPhi2*dPhi2); if ( dRmatch1 < 0.4 && dRmatch2 > 0.4 ) istau1 = true; else if ( dRmatch1 > 0.4 && dRmatch2 < 0.4 ) istau2 = true; else if ( dRmatch1 < 0.4 && dRmatch2 < 0.4 ) { if ( dRmatch1 < dRmatch2 ) { if ( !istau1 ) istau1 = true; else istau2 = true; } else if ( dRmatch2 < dRmatch1 ) { if ( !istau2 ) istau2 = true; else istau1 = true; } } } return make_pair(istau1, istau2); } pair<bool, bool> tauMatchPseudoJet(const Event& event, const PseudoJet& pjet1, PseudoJet& pjet2, double rMatch) { ParticleVector taulep; for ( GenEvent::particle_const_iterator p=event.genEvent()->particles_begin(); p!=event.genEvent()->particles_end(); ++p ) { const PdgId pid = (*p)->pdg_id(); if ( abs(pid) == 15 ) taulep.push_back(Particle(**p)); //applying no pt cut on taus... } if ( taulep.size() == 0 ) return pair<bool, bool>(false, false); bool istau1 = false, istau2 = false; foreach ( const Particle& t, taulep ) { double dEta1 = deltaEta(pjet1.eta(),t.eta()); double dPhi1 = deltaPhi(pjet1.phi(),t.phi()); double dEta2 = deltaEta(pjet2.eta(),t.eta()); double dPhi2 = deltaPhi(pjet2.phi(),t.phi()); double dRmatch1 = sqrt(dEta1*dEta1 + dPhi1*dPhi1); double dRmatch2 = sqrt(dEta2*dEta2 + dPhi2*dPhi2); if ( dRmatch1 < rMatch && dRmatch2 > rMatch ) istau1 = true; else if ( dRmatch1 > rMatch && dRmatch2 < rMatch ) istau2 = true; else if ( dRmatch1 < rMatch && dRmatch2 < rMatch ) { if ( dRmatch1 < dRmatch2 ) { if ( !istau1 ) istau1 = true; else istau2 = true; } else if ( dRmatch2 < dRmatch1 ) { if ( !istau2 ) istau2 = true; else istau1 = true; } } } return make_pair(istau1, istau2); } /// Normalise histograms etc., after the run void finalize() { cout << "This is the high-WpT version of the analysis!" << endl; cout << "Weighted number of events with pT(W)>200GeV: " << nev << endl; cout << "Total weighted number of events: " << sumw << endl; _treeFile->cd(); _treeFile->WriteTObject(_h_precut_PtW); _treeFile->WriteTObject(_h_PtW); _treeFile->WriteTObject(_h_precut_EtaW); _treeFile->WriteTObject(_h_EtaW); _treeFile->WriteTObject(_h_nW); _treeFile->WriteTObject(_h_weights); _treeFile->WriteTObject(_h_xsec); _myTree->Write(); _treeFile->Close(); } //@} private: // Data members like post-cuts event weight counters go here //NTUPLE VARIABLES double var_weight, var_WpT, var_WEta, var_leptonPt, var_leptonEta, var_nuPt, var_nuEta; int var_nAktJets, var_nCAJets, var_AltnCAJets, var_nBMatchJets, var_nBTagJets, var_nBMatchSubJets, var_AltnBMatchSubJets, var_nTotCAJets; double var_JetPt1, var_JetEta1, var_JetPhi1; double var_JetPt2, var_JetEta2, var_JetPhi2; double var_JetPt3, var_JetEta3, var_JetPhi3; bool var_JetB1, var_JetB2, var_JetB1alt, var_JetB2alt, var_JetB1double, var_JetC1alt, var_JetC2alt, var_JetTau1alt, var_JetTau2alt; double var_dRjj; double var_Mj, var_Mjj; double var_SubJetPt1, var_SubJetEta1, var_SubJetPhi1; double var_SubJetPt2, var_SubJetEta2, var_SubJetPhi2; double var_SubJetPt3, var_SubJetEta3, var_SubJetPhi3; bool var_SubJetB1, var_SubJetB2, var_SubJetVarB1, var_SubJetVarB2, var_SubJetC1, var_SubJetC2, var_SubJetTau1, var_SubJetTau2; double var_FatJetPt, var_FatJetEta, var_FatJetPhi, var_dRCentral, var_Rfilt, var_dPhiWJ, var_mu, var_y, var_FatJetRap; double var_FatJet2Pt, var_FatJet2Eta, var_FatJet2Phi, var_FatJet2Rap; int var_nSj; double var_Q, var_x1, var_x2, var_alphaS; int var_id1, var_id2; double var_dEtaWjj, var_dEtaWJ, var_AltdEtaWJ; double var_MWjj, var_MWJ, var_AltMWJ; double var_AltHTca, var_HTca, var_HTakt; private: float nev; float sumw; // float njs, nakt, njs_alt; TH1F* _h_weights; TH1F* _h_xsec; TH1F* _h_PtW; TH1F* _h_precut_PtW; TH1F* _h_EtaW; TH1F* _h_nW; TH1F* _h_precut_EtaW; const char* _treeFileName; TFile* _treeFile; TTree* _myTree; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(exampleAnalysis); }
More information about the Rivet mailing list |