/** *\Class DMAnalysisTreeMaker * * Produces analysis trees from edm-ntuples adding extra variables for resolved and unresolved tops * For custom systematics scenarios * * \Author A. Orso M. Iorio * * *\version $Id: * * */ #include "FWCore/ServiceRegistry/interface/Service.h" #include "FWCore/Framework/interface/EDAnalyzer.h" #include "FWCore/Utilities/interface/InputTag.h" #include "DataFormats/Common/interface/Handle.h" #include "FWCore/Framework/interface/Event.h" #include "FWCore/ParameterSet/interface/ParameterSet.h" #include "FWCore/Utilities/interface/EDMException.h" #include "CommonTools/UtilAlgos/interface/TFileService.h" #include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h" #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h" #include "CondFormats/JetMETObjects/interface/FactorizedJetCorrector.h" #include "CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h" #include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h" #include "DataFormats/Math/interface/LorentzVector.h" #include "DataFormats/Math/interface/deltaPhi.h" #include "DataFormats/Math/interface/deltaR.h" #include "PhysicsTools/Utilities/interface/LumiReWeighting.h" #include #include "CondFormats/BTauObjects/interface/BTagCalibration.h" #include "CondFormats/BTauObjects/interface/BTagCalibrationReader.h" #include "TFile.h" #include "TTree.h" #include "TMath.h" #include #include //using namespace reco; using namespace edm; using namespace std; namespace LHAPDF { void initPDFSet(int nset, const std::string &filename, int member = 0); int numberPDF(int nset); void usePDFMember(int nset, int member); double xfx(int nset, double x, double Q, int fl); double getXmin(int nset, int member); double getXmax(int nset, int member); double getQ2min(int nset, int member); double getQ2max(int nset, int member); void extrapolate(bool extrapolate = true); } class DMAnalysisTreeMaker : public edm::EDAnalyzer { public: explicit DMAnalysisTreeMaker( const edm::ParameterSet & ); private: virtual void beginJob() { std::string distr = "pileUp" + dataPUFile_ + ".root"; LumiWeights_ = edm::LumiReWeighting(distr,"DataPileupHistogram_69mbMinBias.root",std::string("pileup"),std::string("pileup")); LumiWeightsUp_ = edm::LumiReWeighting(distr,"DataPileupHistogram_69mbMinBias_up.root",std::string("pileup"),std::string("pileup")); LumiWeightsDown_ = edm::LumiReWeighting(distr,"DataPileupHistogram_69mbMinBias_down.root",std::string("pileup"),std::string("pileup")); // LumiWeights_ = edm::LumiReWeighting(distr,"PUdata_19468.3.root",std::string("pileup"),std::string("pileup")); } virtual void analyze(const edm::Event &, const edm::EventSetup & ); vector additionalVariables(string); string makeName(string label,string pref,string var); string makeBranchName(string label, string pref, string var); void initializePdf(string centralpdf,string variationpdf); void getEventPdf(); bool getEventTriggers(); void getEventLHEWeights(); bool getMETFilters(); double jetUncertainty(double pt, double eta, string syst); // double smearPt(double pt, double genpt, double eta, string syst); double smear(double pt, double genpt, double eta, string syst); double resolSF(double eta, string syst); double getScaleFactor(double pt, double eta, double partonFlavour, string syst); double muonSF(bool isdata, float pt, float eta, int syst); //------------------ Soureek adding pile-up Info ------------------------------ void getPUSF(); //----------------------------------------------------------------------------- bool isInVector(std::vector v, std::string s); std::vector physObjects; std::vector variablesFloat, variablesInt, singleFloat, singleInt; edm::EDGetTokenT< LHEEventProduct > t_lhes_; edm::EDGetTokenT< GenEventInfoProduct > t_genprod_; edm::EDGetTokenT< std::vector > t_triggerNames_; edm::EDGetTokenT< std::vector > t_triggerBits_; edm::EDGetTokenT< std::vector > t_triggerPrescales_; edm::EDGetTokenT< unsigned int > t_lumiBlock_; edm::EDGetTokenT< unsigned int > t_runNumber_; edm::EDGetTokenT< ULong64_t > t_eventNumber_; edm::EDGetTokenT< bool > t_HBHEFilter_; edm::EDGetTokenT< std::vector > t_metNames_; edm::EDGetTokenT< std::vector > t_metBits_; edm::EDGetTokenT< std::vector > t_pvZ_,t_pvChi2_,t_pvRho_; edm::EDGetTokenT< std::vector > t_pvNdof_; edm::EDGetTokenT t_Rho_; edm::EDGetTokenT t_ntrpu_; //==================================== // edm::LumiReWeighting LumiWeights_, LumiWeightsUp_, LumiWeightsDown_; TTree * treesBase; map trees; std::vector names; std::vector systematics; map< string , float[100] > vfloats_values; map< string , int[100] > vints_values; map< string , string > obs_to_obj; map< string , string > obj_to_pref; map< string , float > float_values; map< string , int > int_values; map< string , int > sizes; map< string , bool > got_label; map< string , int > max_instances; map< int, int > subj_jet_map; map > > h_floats; map > > h_ints; map > h_float; map >h_int; map > > t_floats; map > > t_ints; map > t_float; map >t_int; string mu_label, ele_label, jets_label, met_label, metnohf_label, jetsnohf_label; //MC info: edm::ParameterSet channelInfo; std::string channel; double crossSection, originalEvents; bool useLHEWeights, useLHE, useTriggers,cutOnTriggers, useMETFilters, addPV; bool addLHAPDFWeights; string centralPdfSet,variationPdfSet; std::vector leptonicTriggers, hadronicTriggers, metFilters; int maxPdf, maxWeights; edm::Handle genprod; edm::Handle lhes; edm::Handle > triggerBits; edm::Handle > triggerNames; edm::Handle > triggerPrescales; edm::Handle > metBits; edm::Handle > metNames; edm::Handle HBHE; edm::Handle lumiBlock; edm::Handle runNumber; edm::Handle eventNumber; // float pdf_weights[140]; //float lhe_weights[20]; // std::string lhe_weights_id[20]; edm::Handle > pvZ,pvChi2,pvRho; edm::Handle > pvNdof; int nPV; //----------------------------- Soureek adding for PU info ------------------------- bool doPU_; std::string dataPUFile_; edm::Handle npv, ntrpu; edm::Handle > pubx, puNInt; int nTruePU; edm::LumiReWeighting LumiWeights_, LumiWeightsUp_, LumiWeightsDown_; float puWeight, puWeightUp, puWeightDown; //-------------------------------------------------------------------------------------- //JEC info bool changeJECs = false; bool isData, useMETNoHF; edm::Handle rho; double Rho; std::vector jetScanCuts; std::vector jecPars; JetCorrectorParameters *jecParsL1, *jecParsL2, *jecParsL3, *jecParsL2L3Residuals; JetCorrectionUncertainty *jecUnc; FactorizedJetCorrector *jecCorr; bool isFirstEvent; class BTagWeight { private: int minTags; int maxTags; public: struct JetInfo { JetInfo(float mceff, float datasf) : eff(mceff), sf(datasf) {} float eff; float sf; }; BTagWeight(): minTags(0), maxTags(0) { ; } BTagWeight(int jmin, int jmax) : minTags(jmin) , maxTags(jmax) {} bool filter(int t); float weight(vector jets, int tags); float weightWithVeto(vector jetsTags, int tags, vector jetsVetoes, int vetoes); }; vector jsfscsvt, jsfscsvt_b_tag_up, jsfscsvt_b_tag_down, jsfscsvt_mistag_up, jsfscsvt_mistag_down; vector jsfscsvm, jsfscsvm_b_tag_up, jsfscsvm_b_tag_down, jsfscsvm_mistag_up, jsfscsvm_mistag_down; vector jsfscsvl, jsfscsvl_b_tag_up, jsfscsvl_b_tag_down, jsfscsvl_mistag_up, jsfscsvl_mistag_down; BTagWeight b_csvt_0_tags= BTagWeight(0,0), b_csvt_1_tag= BTagWeight(1,1), b_csvt_2_tags= BTagWeight(2,2); double b_weight_csvt_0_tags, b_weight_csvt_1_tag, b_weight_csvt_2_tags; double b_weight_csvt_0_tags_mistag_up, b_weight_csvt_1_tag_mistag_up, b_weight_csvt_2_tags_mistag_up; double b_weight_csvt_0_tags_mistag_down, b_weight_csvt_1_tag_mistag_down, b_weight_csvt_2_tags_mistag_down; double b_weight_csvt_0_tags_b_tag_down, b_weight_csvt_1_tag_b_tag_down, b_weight_csvt_2_tags_b_tag_down; double b_weight_csvt_0_tags_b_tag_up, b_weight_csvt_1_tag_b_tag_up, b_weight_csvt_2_tags_b_tag_up; BTagWeight b_csvm_0_tags= BTagWeight(0,0), b_csvm_1_tag= BTagWeight(1,1), b_csvm_2_tags= BTagWeight(2,2); double b_weight_csvm_0_tags, b_weight_csvm_1_tag, b_weight_csvm_2_tags; double b_weight_csvm_0_tags_mistag_up, b_weight_csvm_1_tag_mistag_up, b_weight_csvm_2_tags_mistag_up; double b_weight_csvm_0_tags_mistag_down, b_weight_csvm_1_tag_mistag_down, b_weight_csvm_2_tags_mistag_down; double b_weight_csvm_0_tags_b_tag_down, b_weight_csvm_1_tag_b_tag_down, b_weight_csvm_2_tags_b_tag_down; double b_weight_csvm_0_tags_b_tag_up, b_weight_csvm_1_tag_b_tag_up, b_weight_csvm_2_tags_b_tag_up; BTagWeight b_csvl_0_tags= BTagWeight(0,0), b_csvl_1_tag= BTagWeight(1,1), b_csvl_2_tags= BTagWeight(2,2); double b_weight_csvl_0_tags_mistag_up, b_weight_csvl_1_tag_mistag_up, b_weight_csvl_2_tags_mistag_up; double b_weight_csvl_0_tags_mistag_down, b_weight_csvl_1_tag_mistag_down, b_weight_csvl_2_tags_mistag_down; double b_weight_csvl_0_tags_b_tag_down, b_weight_csvl_1_tag_b_tag_down, b_weight_csvl_2_tags_b_tag_down; double b_weight_csvl_0_tags_b_tag_up, b_weight_csvl_1_tag_b_tag_up, b_weight_csvl_2_tags_b_tag_up; double b_weight_csvl_0_tags, b_weight_csvl_1_tag, b_weight_csvl_2_tags; double MCTagEfficiency(string algo, int flavor, double pt); double TagScaleFactor(string algo, int flavor, string syst,double pt); /* // read SF from csv file BTagCalibration btagsf_calib = BTagCalibration("","CSVv2_Aug24.csv"); // central BTagCalibrationReader reader_csvl = BTagCalibrationReader(&btagsf_calib, BTagEntry::OP_LOOSE, "comb", "central"); BTagCalibrationReader reader_csvm = BTagCalibrationReader(&btagsf_calib, BTagEntry::OP_MEDIUM, "comb", "central"); BTagCalibrationReader reader_csvt = BTagCalibrationReader(&btagsf_calib, BTagEntry::OP_TIGHT, "comb", "central"); // btag up BTagCalibrationReader reader_csvlup = BTagCalibrationReader(&btagsf_calib, BTagEntry::OP_LOOSE, "comb", "up"); BTagCalibrationReader reader_csvmup = BTagCalibrationReader(&btagsf_calib, BTagEntry::OP_MEDIUM, "comb", "up"); BTagCalibrationReader reader_csvtup = BTagCalibrationReader(&btagsf_calib, BTagEntry::OP_TIGHT, "comb", "up"); // btag down BTagCalibrationReader reader_csvldown = BTagCalibrationReader(&btagsf_calib, BTagEntry::OP_LOOSE, "comb", "down"); BTagCalibrationReader reader_csvmdown = BTagCalibrationReader(&btagsf_calib, BTagEntry::OP_MEDIUM, "comb", "down"); BTagCalibrationReader reader_csvtdown = BTagCalibrationReader(&btagsf_calib, BTagEntry::OP_TIGHT, "comb", "down"); */ // bool doBTagSF=true; }; DMAnalysisTreeMaker::DMAnalysisTreeMaker(const edm::ParameterSet& iConfig){ mu_label = iConfig.getParameter("muLabel"); ele_label = iConfig.getParameter("eleLabel"); jets_label = iConfig.getParameter("jetsLabel"); met_label = iConfig.getParameter("metLabel"); metnohf_label = iConfig.getParameter("metnohfLabel"); jetsnohf_label = iConfig.getParameter("jetsnohfLabel"); physObjects = iConfig.template getParameter >("physicsObjects"); channelInfo = iConfig.getParameter("channelInfo"); // The physics of the channel, e.g. the cross section, #original events, etc. channel = channelInfo.getParameter("channel"); crossSection = channelInfo.getParameter("crossSection"); originalEvents = channelInfo.getParameter("originalEvents"); useLHEWeights = channelInfo.getUntrackedParameter("useLHEWeights",false); //useLHE = channelInfo.getUntrackedParameter("useLHE",false); useLHE = channelInfo.getUntrackedParameter("useLHE",false); //addLHAPDFWeights = channelInfo.getUntrackedParameter("addLHAPDFWeights",false); addLHAPDFWeights = channelInfo.getUntrackedParameter("addLHAPDFWeights",true); if( useLHE ){ edm::InputTag lhes_ = iConfig.getParameter( "lhes" ); t_lhes_ = consumes< LHEEventProduct >( lhes_ ); } if( addLHAPDFWeights ){ edm::InputTag genprod_ = iConfig.getParameter( "genprod" ); t_genprod_ = consumes( genprod_ ); } useTriggers = iConfig.getUntrackedParameter("useTriggers",true); cutOnTriggers = iConfig.getUntrackedParameter("cutOnTriggers",true); edm::InputTag lumiBlock_ = iConfig.getParameter("lumiBlock"); t_lumiBlock_ = consumes< unsigned int >( lumiBlock_ ); edm::InputTag runNumber_ = iConfig.getParameter("runNumber"); t_runNumber_ = consumes< unsigned int >( runNumber_ ); edm::InputTag eventNumber_ = iConfig.getParameter("eventNumber"); t_eventNumber_ = consumes< ULong64_t >( eventNumber_ ); if(useTriggers){ edm::InputTag triggerBits_ = iConfig.getParameter("triggerBits"); t_triggerBits_ = consumes< std::vector >( triggerBits_ ); edm::InputTag triggerNames_ = iConfig.getParameter("triggerNames"); t_triggerNames_ = consumes< std::vector >( triggerNames_ ); edm::InputTag triggerPrescales_ = iConfig.getParameter("triggerPrescales"); t_triggerPrescales_ = consumes< std::vector >( triggerPrescales_ ); leptonicTriggers= channelInfo.getParameter >("leptonicTriggers"); hadronicTriggers= channelInfo.getParameter >("hadronicTriggers"); } useMETFilters = iConfig.getUntrackedParameter("useMETFilters",true); if(useMETFilters){ metFilters = channelInfo.getParameter >("metFilters"); edm::InputTag metBits_ = iConfig.getParameter("metBits"); t_metBits_ = consumes< std::vector >( metBits_ ); edm::InputTag metNames_ = iConfig.getParameter("metNames"); t_metNames_ = consumes< std::vector >( metNames_ ); edm::InputTag HBHEFilter_ = iConfig.getParameter("HBHEFilter"); t_HBHEFilter_ = consumes< bool >( HBHEFilter_ ); } addPV = iConfig.getUntrackedParameter("addPV",true); changeJECs = iConfig.getUntrackedParameter("changeJECs",false); isData = iConfig.getUntrackedParameter("isData",false); useMETNoHF = iConfig.getUntrackedParameter("useMETNoHF",false); if( changeJECs ) t_Rho_ = consumes( edm::InputTag( "fixedGridRhoFastjetAll" ) ) ; if(addPV || changeJECs){ edm::InputTag pvZ_ = iConfig.getParameter("vertexZ"); t_pvZ_ = consumes< std::vector >( pvZ_ ); edm::InputTag pvChi2_ = iConfig.getParameter("vertexChi2"); t_pvChi2_ = consumes< std::vector >( pvChi2_ ); edm::InputTag pvRho_ = iConfig.getParameter("vertexRho"); t_pvRho_ = consumes< std::vector >( pvRho_ ); edm::InputTag pvNdof_ = iConfig.getParameter("vertexNdof"); t_pvNdof_ = consumes< std::vector< int > >( pvNdof_ ); } //---------------- Soureek adding PU info ----------------------------------------- doPU_=iConfig.getParameter("doPU"); dataPUFile_=iConfig.getParameter("dataPUFile"); if( doPU_ ) t_ntrpu_ = consumes< int >( edm::InputTag( "eventUserData","puNtrueInt" ) ); //--------------------------------------------------------------------------------- if(useLHEWeights){ maxWeights = channelInfo.getUntrackedParameter("maxWeights",9);//Usually we do have 9 weights for the scales, might vary depending on the lhe } if(addLHAPDFWeights){ centralPdfSet = channelInfo.getUntrackedParameter("pdfSet","CT10"); //centralPdfSet = channelInfo.getUntrackedParameter("pdfSet","NNPDF"); variationPdfSet = channelInfo.getUntrackedParameter("pdfSet","CT10"); //variationPdfSet = channelInfo.getUntrackedParameter("pdfSet","NNPDF"); initializePdf(centralPdfSet,variationPdfSet); } systematics = iConfig.getParameter >("systematics"); jetScanCuts = iConfig.getParameter >("jetScanCuts"); std::vector::const_iterator itPsets = physObjects.begin(); bool addNominal=false; for (size_t s = 0; s fs; TFileDirectory DMTrees;// = fs->mkdir( "systematics_trees" ); if(addNominal){ //Service fs; DMTrees = fs->mkdir( "systematics_trees" ); } // treesBase = new TTree("TreeBase", "TreeBase"); trees["noSyst"] = new TTree((channel+"__noSyst").c_str(),(channel+"__noSyst").c_str()); for (;itPsets!=physObjects.end();++itPsets){ int maxI = itPsets->getUntrackedParameter< int >("maxInstances",10); variablesFloat = itPsets->template getParameter >("variablesF"); variablesInt = itPsets->template getParameter >("variablesI"); singleFloat = itPsets->template getParameter >("singleF"); singleInt = itPsets->template getParameter >("singleI"); string namelabel = itPsets->getParameter< string >("label"); string nameprefix = itPsets->getParameter< string >("prefix"); bool saveBaseVariables = itPsets->getUntrackedParameter("saveBaseVariables",true); std::vector toSave= itPsets->getParameter >("toSave"); std::vector::const_iterator itF = variablesFloat.begin(); std::vector::const_iterator itI = variablesInt.begin(); std::vector::const_iterator itsF = singleFloat.begin(); std::vector::const_iterator itsI = singleInt.begin(); stringstream max_instance_str; max_instance_str<instance()+"_"+itF->label(); string nameinstance=itF->instance(); string nameshort=itF->instance(); //if(nameobs!=itF->label())cout<<" warning! label not matching the module! check if members of pset are consistent. If intentional , ignore this warning"; //cout << "nameobs "<< nameobs<< " name " << name <<" nameshort " <instance()<< " is in vector? "<< isInVector(toSave,itF->instance())<< endl; string nametobranch = makeBranchName(namelabel,prefix,nameinstance); name = nametobranch; nameshort = nametobranch; // cout << " name for the branch is now: "<< nametobranch<instance())) trees["noSyst"]->Branch(nameshort.c_str(), &vfloats_values[name],(nameshort+"["+max_instance_str.str()+"]/F").c_str()); if(saveBaseVariables|| isInVector(toSave,itF->instance())) trees["noSyst"]->Branch(nameshort.c_str(), &vfloats_values[name],(nameshort+"["+max_instance_str.str()+"]/F").c_str()); names.push_back(name); obs_to_obj[name] = nameobs; obj_to_pref[nameobs] = prefix; t_floats[ name ] = consumes< std::vector >( *itF ); cout << " branching name "<< name<< " for obs "<< nameobs << " instance "<< nameinstance << endl; } for (;itI != variablesInt.end();++itI){ string name=itI->instance()+"_"+itI->label(); string nameshort=itF->instance(); string nametobranch = makeBranchName(namelabel,prefix,nameshort); name = nametobranch; nameshort = nametobranch; if(saveBaseVariables|| isInVector(toSave,itI->instance())) trees["noSyst"]->Branch(nameshort.c_str(), &vints_values[name],(nameshort+"["+max_instance_str.str()+"]/I").c_str()); names.push_back(name); obs_to_obj[name] = nameobs; obj_to_pref[nameobs] = prefix; t_ints[ name ] = consumes< std::vector >( *itI ); } if (variablesFloat.size()>0){ string nameshortv = namelabel; vector extravars = additionalVariables(nameshortv); for(size_t addv = 0; addv < extravars.size();++addv){ string name = nameshortv+"_"+extravars.at(addv); if (saveBaseVariables || isInVector(toSave, extravars.at(addv)) || isInVector(toSave, "allExtra") )trees["noSyst"]->Branch(name.c_str(), &vfloats_values[name],(name+"["+max_instance_str.str()+"]/F").c_str()); obs_to_obj[name] = nameobs; obj_to_pref[nameobs] = prefix; } } names.push_back(nameobs); cout << "size part: nameobs is "<< nameobs<Branch((nameobs+"_size").c_str(), &sizes[nameobs]); //Initialize single pset objects for (;itsF != singleFloat.end();++itsF){ string name=itsF->instance()+itsF->label(); string nameshort=itsF->instance(); string nametobranch = makeBranchName(namelabel,prefix,nameshort); name = nametobranch; nameshort = nametobranch; t_float[ name ] = consumes< float >( *itsF ); if(saveBaseVariables|| isInVector(toSave,itsF->instance()))trees["noSyst"]->Branch(nameshort.c_str(), &float_values[name]); } for (;itsI != singleInt.end();++itsI){ string name=itsI->instance()+itsI->label(); string nameshort=itsI->instance(); string nametobranch = makeBranchName(namelabel,prefix,nameshort); name = nametobranch; nameshort = nametobranch; t_int[ name ] = consumes< int >( *itsI ); if(saveBaseVariables|| isInVector(toSave,itsI->instance()))trees["noSyst"]->Branch(nameshort.c_str(), &int_values[name]); } } string nameshortv= "Event"; vector extravars = additionalVariables(nameshortv); for(size_t addv = 0; addv < extravars.size();++addv){ string name = nameshortv+"_"+extravars.at(addv); if (name.find("EventNumber") != std::string::npos) trees["noSyst"]->Branch(name.c_str(), &int_values[name],(name+"/I").c_str()); else trees["noSyst"]->Branch(name.c_str(), &float_values[name],(name+"/F").c_str()); } //Prepare the trees cloning all branches and setting the correct names/titles: if(!addNominal){ DMTrees = fs->mkdir( "systematics_trees" ); } for(size_t s=0;s< systematics.size();++s){ std::string syst = systematics.at(s); if(syst=="noSyst")continue; trees[syst]= trees["noSyst"]->CloneTree(); //trees[syst]= treesBase->CloneTree(); trees[syst]->SetName((channel+"__"+syst).c_str()); trees[syst]->SetTitle((channel+"__"+syst).c_str()); } if(isData) jecParsL1 = new JetCorrectorParameters("Fall15_25nsV1_DATA_L1FastJet_AK4PFchs.txt"); else jecParsL1 = new JetCorrectorParameters("Fall15_25nsV1_MC_L1FastJet_AK4PFchs.txt"); jecParsL2 = new JetCorrectorParameters("Fall15_25nsV1_MC_L2Relative_AK4PFchs.txt"); jecParsL3 = new JetCorrectorParameters("Fall15_25nsV1_MC_L3Absolute_AK4PFchs.txt"); jecParsL2L3Residuals = new JetCorrectorParameters("Summer15_25nsV7_DATA_L2L3Residual_AK4PFchs.txt"); jecPars.push_back(*jecParsL1); jecPars.push_back(*jecParsL2); jecPars.push_back(*jecParsL3); if(isData) jecPars.push_back(*jecParsL2L3Residuals); jecCorr = new FactorizedJetCorrector(jecPars); jecUnc = new JetCorrectionUncertainty(*(new JetCorrectorParameters("Summer15_25nsV7_DATA_UncertaintySources_AK4PFchs.txt", "Total"))); // if(addNominal) systematics.push_back("noSyst"); } void DMAnalysisTreeMaker::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) { std::vector::const_iterator itPsets = physObjects.begin(); // if(addLHAPDFWeights){ // event info iEvent.getByToken(t_lumiBlock_,lumiBlock ); iEvent.getByToken(t_runNumber_,runNumber ); iEvent.getByToken(t_eventNumber_,eventNumber ); if(useLHE){ iEvent.getByToken(t_lhes_, lhes); } if(addLHAPDFWeights){ iEvent.getByToken(t_genprod_, genprod); } //Part 0: trigger preselection if(useTriggers){ iEvent.getByToken(t_triggerBits_,triggerBits ); iEvent.getByToken(t_triggerNames_,triggerNames ); iEvent.getByToken(t_triggerPrescales_,triggerPrescales ); bool triggerOr = getEventTriggers(); if(cutOnTriggers && !triggerOr) return; } if(useMETFilters){ iEvent.getByToken(t_metBits_,metBits ); iEvent.getByToken(t_metNames_,metNames ); iEvent.getByToken(t_HBHEFilter_ ,HBHE); if(isFirstEvent){ for(size_t bt = 0; bt < metNames->size();++bt){ std::string tname = metNames->at(bt); cout << "test tname "<< tname << " passes "<< metBits->at(bt)<< endl; } //cout << "test tname "<< tname <size(); } /*std::cout << " initTriggers "<< endl; for(size_t bt = 0; bt < triggerNames->size();++bt){ std::string tname = triggerNames->at(bt); float bit = triggerBits->at(bt); int presc = triggerPrescales->at(bt); std::cout << "name "<< tname << " bit "<< bit << " prescale "<template getParameter >("variablesF"); variablesInt = itPsets->template getParameter >("variablesI"); singleFloat = itPsets->template getParameter >("singleF"); singleInt = itPsets->template getParameter >("singleI"); std::vector::const_iterator itF = variablesFloat.begin(); std::vector::const_iterator itI = variablesInt.begin(); std::vector::const_iterator itsF = singleFloat.begin(); std::vector::const_iterator itsI = singleInt.begin(); string namelabel = itPsets->getParameter< string >("label"); string nameprefix = itPsets->getParameter< string >("prefix"); size_t maxInstance=(size_t)max_instances[namelabel]; //Vectors of floats for (;itF != variablesFloat.end();++itF){ string varname=itF->instance(); string name = makeBranchName(namelabel,nameprefix,varname); //string namelabel; float tmp =1.0; iEvent.getByToken(t_floats[name] ,h_floats[name]); // cout << "name "<< name <size()){tmp = h_floats[name]->at(fi);} else { tmp = -9999.; } // cout << " setting name "<< name<< " at instance "<< fi <<" to value "<< tmp <size(); // cout<< " size for "<< namelabel <<" is then "<< sizes[namelabel]<instance(); string name = makeBranchName(namelabel,nameprefix,varname); int tmp = 1; iEvent.getByToken(t_ints[name] ,h_ints[name]); for (size_t fi = 0;fi < maxInstance;++fi){ if(fi size()){tmp = h_ints[name]->at(fi);} else { tmp = -9999.; } vints_values[name][fi]=tmp; } } // std::cout << " checkpoint ints"<instance(); string name = makeBranchName(namelabel,nameprefix,varname); iEvent.getByToken(t_float[name],h_float[name]); float_values[name]=*h_float[name]; } for (;itsI != singleInt.end();++itsI){ string varname=itsI->instance(); string name = makeBranchName(namelabel,nameprefix,varname); iEvent.getByToken(t_int[name],h_int[name]); int_values[name]=*h_int[name]; } // std::cout << " checkpoint singles"< electrons; vector muons; vector leptons; vector loosemuons; vector jets; vector muons_t; for (size_t s = 0; s< systematics.size();++s){ int ncsvl_tags=0,ncsvt_tags=0,ncsvm_tags=0;//, njets_tottag=0; getEventTriggers(); electrons.clear(); muons.clear(); loosemuons.clear(); jets.clear(); string syst = systematics.at(s); jsfscsvt.clear(); jsfscsvt_b_tag_up.clear(); jsfscsvt_b_tag_down.clear(); jsfscsvt_mistag_up.clear(); jsfscsvt_mistag_down.clear(); jsfscsvm.clear(); jsfscsvm_b_tag_up.clear(); jsfscsvm_b_tag_down.clear(); jsfscsvm_mistag_up.clear(); jsfscsvm_mistag_down.clear(); jsfscsvl.clear(); jsfscsvl_b_tag_up.clear(); jsfscsvl_b_tag_down.clear(); jsfscsvl_mistag_up.clear(); jsfscsvl_mistag_down.clear(); //---------------- Soureek Adding PU Info ------------------------------ if(doPU_){ iEvent.getByToken(t_ntrpu_,ntrpu); nTruePU=*ntrpu; getPUSF(); } //std::cout<<"Check for PU re-weighting 2"<0 && pt> 22 && abs(eta) < 2.1 && iso <0.06){ ++float_values["Event_nTightMuons"]; TLorentzVector muon; muon.SetPtEtaPhiE(pt, eta, phi, energy); muons.push_back(muon); muons_t.push_back(muon); leptons.push_back(muon); mu_sf *= muonSF(isData,pt,eta,0); mu_sf_up *= muonSF(isData,pt,eta,1); mu_sf_down *= muonSF(isData,pt,eta,-1); } if(isLoose>0 && pt> 10 && abs(eta) < 2.1 && iso<0.2){ ++float_values["Event_nLooseMuons"]; } if(isSoft>0 && pt> 30 && abs(eta) < 2.4){ ++float_values["Event_nSoftMuons"]; } if(isLoose>0 && pt > 15){ TLorentzVector muon; muon.SetPtEtaPhiE(pt, eta, phi, energy); loosemuons.push_back(muon); } } float_values["Event_mu_eff"]=mu_sf; float_values["Event_mu_eff_up"]=mu_sf_up; float_values["Event_mu_eff_down"]=mu_sf_down; /************************** Electrons: **************************/ for(int el = 0;el < max_instances[ele_label] ;++el){ string pref = obj_to_pref[ele_label]; float pt = vfloats_values[makeName(ele_label,pref,"Pt")][el]; // if(pt<0)continue; float isTight = vfloats_values[makeName(ele_label,pref,"isTight")][el]; float isLoose = vfloats_values[makeName(ele_label,pref,"isLoose")][el]; float isMedium = vfloats_values[makeName(ele_label,pref,"isMedium")][el]; float isVeto = vfloats_values[makeName(ele_label,pref,"isVeto")][el]; //float vidTight = vfloats_values[makeName(ele_label,pref,"vidTight")][el]; //float vidLoose = vfloats_values[makeName(ele_label,pref,"vidLoose")][el]; //float vidMedium = vfloats_values[makeName(ele_label,pref,"vidMedium")][el]; //float vidVeto = vfloats_values[makeName(ele_label,pref,"vidVeto")][el]; //float vidHEEP = vfloats_values[makeName(ele_label,pref,"vidHEEP")][el]; float eta = vfloats_values[makeName(ele_label,pref,"Eta")][el]; float scEta = vfloats_values[makeName(ele_label,pref,"scEta")][el]; float phi = vfloats_values[makeName(ele_label,pref,"Phi")][el]; float energy = vfloats_values[makeName(ele_label,pref,"E")][el]; float iso = vfloats_values[makeName(ele_label,pref,"Iso03")][el]; //cout << "######" << endl; //std::cout << " electron #"<0.0 && iso < 0.074355 ; } //is barrel electron if (fabs(scEta)>1.479 && fabs(scEta)<2.5){ passesTightCuts = isTight >0.0 && iso < 0.090185 ; } if(pt> 30 && abs(eta) < 2.5 && passesTightCuts){ TLorentzVector ele; ele.SetPtEtaPhiE(pt, eta, phi, energy); double minDR=999; double deltaR = 999; for (size_t m = 0; m < (size_t)loosemuons.size(); ++m){ deltaR = ele.DeltaR(loosemuons[m]); minDR = min(minDR, deltaR); } if(!loosemuons.size()) minDR=999; if(minDR>0.1){ electrons.push_back(ele); leptons.push_back(ele); ++float_values["Event_nTightElectrons"]; } else {passesDRmu = false;} } if(isLoose>0 && pt> 30 && abs(eta) < 2.5){ ++float_values["Event_nLooseElectrons"]; } if(isMedium>0 && pt> 30 && abs(eta) < 2.5){ ++float_values["Event_nMediumElectrons"]; } if(isVeto>0 && pt> 20 && abs(eta) < 2.5){ ++float_values["Event_nVetoElectrons"]; } vfloats_values[ele_label+"_PassesDRmu"][el]=(float)passesDRmu; } /************************** MET: **************************/ // cout << " met "<0){ // ptCorr = smearPt(pt, genpt, eta, syst); // float unc = jetUncertainty(ptCorr,eta,syst); // //cout << "genpt? "<< genpt <<" pt ? "<< pt<<" ptcorr? "<< ptCorr<<"unc? "<< unc<0){ TLorentzVector jetUncorr; jetUncorr.SetPtEtaPhiE(pt,eta,phi,energy); // cout << "jet "<< j <<" standard pt "<setJetEta(jetUncorr.Eta()); jecCorr->setJetPt(jetUncorr.Pt()); jecCorr->setJetA(area); jecCorr->setRho(Rho); jecCorr->setNPV(nPV); double recorr = jecCorr->getCorrection(); jetUncorr = jetUncorr *recorr; // cout << " recorrection "< 0.935; bool isCSVM = csv > 0.800; bool isCSVL = csv > 0.460; vfloats_values[jets_label+"_IsCSVT"][j]=isCSVT; vfloats_values[jets_label+"_IsCSVM"][j]=isCSVM; vfloats_values[jets_label+"_IsCSVL"][j]=isCSVL; float bsf = getScaleFactor(ptCorr,eta,partonFlavour,"noSyst"); float bsfup = getScaleFactor(ptCorr,eta,partonFlavour,"up"); float bsfdown = getScaleFactor(ptCorr,eta,partonFlavour,"down"); vfloats_values[jets_label+"_BSF"][j]=bsf; vfloats_values[jets_label+"_BSFUp"][j]=bsfup; vfloats_values[jets_label+"_BSFDown"][j]=bsfdown; //***** JET ID **** bool passesID = true; if(!(jecscale*energy > 0))passesID = false; else{ double enZero=energy*jecscale; float nDau = vfloats_values[makeName(jets_label,pref,"numberOfDaughters")][j]; float mu_frac = vfloats_values[makeName(jets_label,pref,"MuonEnergy")][j]/enZero; float chMulti = vfloats_values[makeName(jets_label,pref,"chargedMultiplicity")][j]; float chHadEnFrac = vfloats_values[makeName(jets_label,pref,"chargedHadronEnergy")][j]/enZero; float chEmEnFrac = vfloats_values[makeName(jets_label,pref,"chargedEmEnergy")][j]/enZero; float neuEmEnFrac = vfloats_values[makeName(jets_label,pref,"neutralEmEnergy")][j]/enZero; float neuHadEnFrac = vfloats_values[makeName(jets_label,pref,"neutralHadronEnergy")][j]/enZero; float neuMulti = vfloats_values[makeName(jets_label,pref,"neutralMultiplicity")][j]; // recommended at 8 TeV //passesID = (nDau >1.0 && fabs(eta) < 4.7 && (fabs(eta)>=2.4 ||(chHadEnFrac>0.0 && chMulti>0 && chEmEnFrac<0.99) ) && neuEmEnFrac<0.99 && neuHadEnFrac <0.99 && mu_frac < 0.8); // temp fix for HF issues //passesID = (nDau >1.0 && fabs(eta) < 4.7 && (fabs(eta)>=2.4 ||(chHadEnFrac>0.0 && chMulti>0 && chEmEnFrac<0.99 && neuEmEnFrac<0.99 && neuHadEnFrac <0.99) ) && mu_frac < 0.8); // new recommendations from 14/08 if (fabs(eta)<=3.0) passesID = (nDau >1.0 && (fabs(eta)>=2.4 ||(chHadEnFrac>0.0 && chMulti>0 && chEmEnFrac<0.99) ) && neuEmEnFrac<0.99 && neuHadEnFrac <0.99); else passesID = (neuMulti > 10 && neuEmEnFrac < 0.9 ); vfloats_values[jets_label+"_JetID_numberOfDaughters"][j]=nDau; vfloats_values[jets_label+"_JetID_muonEnergyFraction"][j]=mu_frac; vfloats_values[jets_label+"_JetID_chargedMultiplicity"][j]=chMulti; vfloats_values[jets_label+"_JetID_chargedHadronEnergyFraction"][j]=chHadEnFrac; vfloats_values[jets_label+"_JetID_chargedEmEnergyFraction"][j]=chEmEnFrac; vfloats_values[jets_label+"_JetID_neutralEmEnergyFraction"][j]=neuEmEnFrac; vfloats_values[jets_label+"_JetID_neutralHadronEnergyFraction"][j]=neuHadEnFrac; vfloats_values[jets_label+"_JetID_neutralMultiplicity"][j]=neuMulti; /* cout << "### JET ID ###" << endl; cout << "ndau " << nDau << endl; cout << "mu_frac " << mu_frac << endl; cout << "chMulti " << chMulti << endl; cout << "chHadEnFrac " << chHadEnFrac << endl; cout << "chEmEnFrac " << chEmEnFrac << endl; cout << "neuEmEnFrac " << neuEmEnFrac << endl; cout << "neuHadEnFrac " << neuHadEnFrac << endl; cout << "passes loose id " << int(passesID) << endl; */ } vfloats_values[jets_label+"_PassesID"][j]=(float)passesID; //remove overlap with tight electrons/muons double minDR=9999; double minDRtight=9999; double minDRThrEl=0.3; double minDRThrMu=0.3; bool passesDR=true; bool passesDRtight=true; for (size_t e = 0; e < (size_t)electrons.size(); ++e){ minDR = min(minDR,deltaR(math::PtEtaPhiELorentzVector(electrons.at(e).Pt(),electrons.at(e).Eta(),electrons.at(e).Phi(),electrons.at(e).Energy() ) ,math::PtEtaPhiELorentzVector(ptCorr, eta, phi, energyCorr))); if(minDR jetval && fabs(eta) < 4.7); if(!passesID || !passesCut || !passesDR) continue; if(ji==0){ vfloats_values[jets_label+"_IsTight"][j]=1.0; TLorentzVector jet; jet.SetPtEtaPhiE(ptCorr, eta, phi, energyCorr); jets.push_back(jet); //double MCTagEFFiciency(int flavor); // double TagScaleFactor(int flavor, string syst); // jsfscsvt_b_tag_up; // jsfscsvt_b_tag_down, // jsfscsvt_mistag_up, // jsfscsvt_mistag_down; // partonFlavour = vfloats_values[makeName(jets_label,pref,"PartonFlavour")][j]; } // cout << " syst "<< syst<< " jet "<< j << " pt "<< ptCorr <<"cut "<< jetScanCuts.at(ji)<< " extra jet with pt "<< ptCorr<< "eventNJets before is" << float_values["Event_nJets"+j_n.str()]<< " csv "<< csv<< " isCSVM? "<< isCSVM< 40 && fabs(eta) < 2.4) { float_values["Event_nCSVTJets"+j_n.str()]+=1.0; ncsvt_tags +=1; } //if(isCSVL && passesCut && fabs(eta) < 2.4) { if(isCSVL && passesID && passesDRtight && ptCorr > 40 && fabs(eta) < 2.4) { ncsvl_tags +=1; } //if(isCSVM && passesCut && fabs(eta) < 2.4) { if(isCSVM && passesID && passesDRtight && ptCorr > 40 && fabs(eta) < 2.4) { float_values["Event_nCSVMJets"+j_n.str()]+=1.0; if(ji==0){ ncsvm_tags +=1; // cout << " jet "<< j<< " isBJET "<0){ TLorentzVector jetUncorr; jetUncorr.SetPtEtaPhiE(pt,eta,phi,energy); // cout << "jet "<< j <<" standard pt "< 10){ DUnclusteredMETPx+=jetUncorr.Pt()*cos(phi); DUnclusteredMETPy+=jetUncorr.Pt()*sin(phi); } // DUnclusteredMETPx=0.1*jetUncorr.Pt()*cos(phi); // DUnclusteredMETPy=0.1*jetUncorr.Pt()*sin(phi); if(changeJECs){ jecCorr->setJetEta(jetUncorr.Eta()); jecCorr->setJetPt(jetUncorr.Pt()); jecCorr->setJetA(area); jecCorr->setRho(Rho); jecCorr->setNPV(nPV); double recorr = jecCorr->getCorrection(); jetUncorr = jetUncorr *recorr; // cout << " recorrection "<0)metphiCorr = atan(metPy/metPx)+3.141592; if(metPy<0)metphiCorr = atan(metPy/metPx)-3.141592; } else metphiCorr = (atan(metPy/metPx)); vfloats_values[met_label+"_CorrPhi"][0]=metphiCorr; //BTagging part if(doBTagSF){ //CSVT //0 tags b_weight_csvt_0_tags = b_csvt_0_tags.weight(jsfscsvt, ncsvt_tags); b_weight_csvt_0_tags_mistag_up = b_csvt_0_tags.weight(jsfscsvt_mistag_up, ncsvt_tags); b_weight_csvt_0_tags_mistag_down = b_csvt_0_tags.weight(jsfscsvt_mistag_down, ncsvt_tags); b_weight_csvt_0_tags_b_tag_up = b_csvt_0_tags.weight(jsfscsvt_b_tag_up, ncsvt_tags); b_weight_csvt_0_tags_b_tag_down = b_csvt_0_tags.weight(jsfscsvt_b_tag_down, ncsvt_tags); //1 tag b_weight_csvt_1_tag = b_csvt_1_tag.weight(jsfscsvt, ncsvt_tags); b_weight_csvt_1_tag_mistag_up = b_csvt_1_tag.weight(jsfscsvt_mistag_up, ncsvt_tags); b_weight_csvt_1_tag_mistag_down = b_csvt_1_tag.weight(jsfscsvt_mistag_down, ncsvt_tags); b_weight_csvt_1_tag_b_tag_up = b_csvt_1_tag.weight(jsfscsvt_b_tag_up, ncsvt_tags); b_weight_csvt_1_tag_b_tag_down = b_csvt_1_tag.weight(jsfscsvt_b_tag_down, ncsvt_tags); //cout <<"w1t check: is"<< b_weight_csvt_1_tag<hepeup().XWGTUP; float LHEWeightSign = LHE_weight/fabs(LHE_weight); float_values["Event_LHEWeightSign"]=LHEWeightSign; float_values["Event_LHEWeight"]=LHE_weight; } float weightLumi = crossSection/originalEvents; float_values["Event_weight"]=weightLumi; //Part 3: filling the additional variables //Reset event weights/#objects if(useLHEWeights){ getEventLHEWeights(); } if(addLHAPDFWeights){ getEventPdf(); } if(addPV){ float nGoodPV = 0.0; for (size_t v = 0; v < pvZ->size();++v){ bool isGoodPV = ( fabs(pvZ->at(v)) < 24.0 && pvNdof->at(v) > 4.0 && pvRho->at(v) <2.0 ); if (isGoodPV)nGoodPV+=1.0; } float_values["Event_nGoodPV"]=(float)(nGoodPV); float_values["Event_nPV"]=(float)(nPV); } // if(useMETFilters){ // getMETFilters(); // } float_values["Event_passesHBHE"]=(float)(*HBHE); //technical event information int_values["Event_EventNumber"]=(*eventNumber); //float_values["Event_EventNumber"]=(float)(*eventNumber); float_values["Event_LumiBlock"]=*lumiBlock; float_values["Event_RunNumber"]=*runNumber; // std::cout << " event ntight muons " << float_values["Event_nTightMuons"]<< std::endl; // std::cout << " event ntight electrons " << float_values["Event_nTightElectrons"]<< std::endl; // cout << "before filling jets label "<< jets_label << " maxinstances "<< max_instances[jets_label]<< "size "<< sizes[jets_label]<Fill(); //Reset event weights/#objects string nameshortv= "Event"; vector extravars = additionalVariables(nameshortv); for(size_t addv = 0; addv < extravars.size();++addv){ string name = nameshortv+"_"+extravars.at(addv); // std::cout << " resetting variable "<< name<< " before "<< float_values[name]; // bool isTriggerVar // if(isTriggerVar) continue; float_values[name]=-1; // std::cout << " after "<< Float_values[name]<Fill(); } string DMAnalysisTreeMaker::makeBranchName(string label, string pref, string var){ string outVar = var; size_t prefPos=outVar.find(pref); size_t prefLength = pref.length(); outVar.replace(prefPos,prefLength,label+"_"); return outVar; } string DMAnalysisTreeMaker::makeName(string label,string pref,string var){ return label+"_"+var; // string outVar = var; // size_t prefPos=outVar.find(pref); // size_t prefLength = pref.length(); // std::cout << " outvar is "<< outVar<< " prefpos "<< prefPos<< " length "<< prefLength<< endl; // std::cout << "it is " << label+"_"+var< DMAnalysisTreeMaker::additionalVariables(string object){ vector addvar; bool ismuon=object.find("muon")!=std::string::npos; bool iselectron=object.find("electron")!=std::string::npos; bool ismet=object.find("met")!=std::string::npos; bool isjet=object.find("jet")!=std::string::npos && object.find("AK4")!=std::string::npos; bool isak8=object.find("jet")!=std::string::npos && object.find("AK8")!=std::string::npos && object.find("sub")==std::string::npos; bool isak8subjet=object.find("jet")!=std::string::npos && object.find("AK8")!=std::string::npos && object.find("sub")!=std::string::npos; bool isevent=object.find("Event")!=std::string::npos; //bool isResolvedTopHad=object.find("resolvedTopHad")!=std::string::npos; //bool isResolvedTopSemiLep=object.find("resolvedTopSemiLep")!=std::string::npos; if(ismuon || iselectron){ //addvar.push_back("SFTrigger"); //addvar.push_back("SFReco"); //addvar.push_back("isQCD"); // addvar.push_back("isTightOffline"); // addvar.push_back("isLooseOffline"); } if(iselectron){ //addvar.push_back("PassesDRmu"); } if(ismet){ // addvar.push_back("Pt"); // addvar.push_back("Eta"); // addvar.push_back("Phi"); // addvar.push_back("E"); addvar.push_back("CorrPt"); addvar.push_back("CorrPhi"); } if(isjet){ addvar.push_back("CorrPt"); // addvar.push_back("CorrEta"); // addvar.push_back("CorrPhi"); addvar.push_back("CorrE"); addvar.push_back("MinDR"); addvar.push_back("IsCSVT"); addvar.push_back("IsCSVM"); addvar.push_back("IsCSVL"); // addvar.push_back("BSF"); // addvar.push_back("BSFUp"); // addvar.push_back("BSFDown"); addvar.push_back("PassesID"); addvar.push_back("PassesDR"); addvar.push_back("PassesDRtight"); addvar.push_back("CorrMass"); addvar.push_back("IsTight"); addvar.push_back("IsLoose"); // addvar.push_back("CorrNJets"); // addvar.push_back("CorrPartonFlavour"); addvar.push_back("JetID_numberOfDaughters"); addvar.push_back("JetID_muonEnergyFraction"); addvar.push_back("JetID_chargedMultiplicity"); addvar.push_back("JetID_chargedHadronEnergyFraction"); addvar.push_back("JetID_chargedEmEnergyFraction"); addvar.push_back("JetID_neutralEmEnergyFraction"); addvar.push_back("JetID_neutralHadronEnergyFraction"); addvar.push_back("JetID_neutralMultiplicity"); } if(isak8){ addvar.push_back("CorrPt"); addvar.push_back("CorrE"); addvar.push_back("isType1"); addvar.push_back("isType2"); addvar.push_back("TopPt"); addvar.push_back("TopEta"); addvar.push_back("TopPhi"); addvar.push_back("TopE"); addvar.push_back("TopMass"); addvar.push_back("TopWMass"); addvar.push_back("nJ"); addvar.push_back("nCSVM"); } if(isak8subjet){ ;// addvar.push_back("CorrPt"); } /* if(isResolvedTopHad ){ addvar.push_back("Pt"); addvar.push_back("Eta"); addvar.push_back("Phi"); addvar.push_back("E"); addvar.push_back("Mass"); addvar.push_back("WMass"); addvar.push_back("BMPhi"); addvar.push_back("WMPhi"); addvar.push_back("TMPhi"); addvar.push_back("WBPhi"); addvar.push_back("IndexB"); addvar.push_back("IndexJ1"); addvar.push_back("IndexJ2"); } if(isResolvedTopSemiLep ){ addvar.push_back("Pt"); addvar.push_back("Eta"); addvar.push_back("Phi"); addvar.push_back("E"); addvar.push_back("Mass"); addvar.push_back("MT"); addvar.push_back("LBMPhi"); addvar.push_back("LMPhi"); addvar.push_back("LBPhi"); addvar.push_back("BMPhi"); addvar.push_back("TMPhi"); addvar.push_back("IndexL"); addvar.push_back("LeptonFlavour"); addvar.push_back("IndexB"); } */ if(isevent){ //addvar.push_back("weight"); addvar.push_back("nTightMuons"); addvar.push_back("nSoftMuons"); addvar.push_back("nLooseMuons"); addvar.push_back("nTightElectrons"); addvar.push_back("nMediumElectrons"); addvar.push_back("nLooseElectrons"); addvar.push_back("nVetoElectrons"); //addvar.push_back("nElectronsSF"); //addvar.push_back("mt"); //addvar.push_back("Mt2w"); //addvar.push_back("category"); //addvar.push_back("nMuonsSF"); //addvar.push_back("nCSVTJets"); //addvar.push_back("nCSVMJets"); //addvar.push_back("nCSVLJets"); //addvar.push_back("nTightJets"); //addvar.push_back("nLooseJets"); //addvar.push_back("nType1TopJets"); //addvar.push_back("nType2TopJets"); addvar.push_back("bWeight0CSVT"); addvar.push_back("bWeight1CSVT"); addvar.push_back("bWeight2CSVT"); addvar.push_back("bWeight0CSVM"); addvar.push_back("bWeight1CSVM"); addvar.push_back("bWeight2CSVM"); addvar.push_back("bWeight0CSVL"); addvar.push_back("bWeight1CSVL"); addvar.push_back("bWeight2CSVL"); addvar.push_back("bWeightMisTagDown0CSVT"); addvar.push_back("bWeightMisTagDown1CSVT"); addvar.push_back("bWeightMisTagDown2CSVT"); addvar.push_back("bWeightMisTagDown0CSVM"); addvar.push_back("bWeightMisTagDown1CSVM"); addvar.push_back("bWeightMisTagDown2CSVM"); addvar.push_back("bWeightMisTagDown0CSVL"); addvar.push_back("bWeightMisTagDown1CSVL"); addvar.push_back("bWeightMisTagDown2CSVL"); addvar.push_back("bWeightMisTagUp0CSVT"); addvar.push_back("bWeightMisTagUp1CSVT"); addvar.push_back("bWeightMisTagUp2CSVT"); addvar.push_back("bWeightMisTagUp0CSVM"); addvar.push_back("bWeightMisTagUp1CSVM"); addvar.push_back("bWeightMisTagUp2CSVM"); addvar.push_back("bWeightMisTagUp0CSVL"); addvar.push_back("bWeightMisTagUp1CSVL"); addvar.push_back("bWeightMisTagUp2CSVL"); addvar.push_back("bWeightBTagUp0CSVT"); addvar.push_back("bWeightBTagUp1CSVT"); addvar.push_back("bWeightBTagUp2CSVT"); addvar.push_back("bWeightBTagUp0CSVM"); addvar.push_back("bWeightBTagUp1CSVM"); addvar.push_back("bWeightBTagUp2CSVM"); addvar.push_back("bWeightBTagUp0CSVL"); addvar.push_back("bWeightBTagUp1CSVL"); addvar.push_back("bWeightBTagUp2CSVL"); addvar.push_back("bWeightBTagDown0CSVT"); addvar.push_back("bWeightBTagDown1CSVT"); addvar.push_back("bWeightBTagDown2CSVT"); addvar.push_back("bWeightBTagDown0CSVM"); addvar.push_back("bWeightBTagDown1CSVM"); addvar.push_back("bWeightBTagDown2CSVM"); addvar.push_back("bWeightBTagDown0CSVL"); addvar.push_back("bWeightBTagDown1CSVL"); addvar.push_back("bWeightBTagDown2CSVL"); addvar.push_back("LHEWeightSign"); addvar.push_back("LHEWeight"); addvar.push_back("EventNumber"); addvar.push_back("LumiBlock"); addvar.push_back("RunNumber"); addvar.push_back("mu_eff"); addvar.push_back("mu_eff_up"); addvar.push_back("mu_eff_down"); if(addPV){ addvar.push_back("nPV"); addvar.push_back("nGoodPV"); } /* for (size_t j = 0; j < (size_t)jetScanCuts.size(); ++j){ stringstream j_n; double jetval = jetScanCuts.at(j); j_n << "Cut" <size();++bt){ std::string tname = metNames->at(bt); // cout << "test tname "<at(bt)>0); float_values["Event_passes"+fname]=metBits->at(bt); } } } float_values["Event_passesMETFilters"]=(float)METFilterAND; return METFilterAND; } bool DMAnalysisTreeMaker::getEventTriggers(){ bool leptonOR=false,hadronOR=false; for(size_t lt =0; lt< leptonicTriggers.size();++lt){ string lname = leptonicTriggers.at(lt); for(size_t bt = 0; bt < triggerNames->size();++bt){ std::string tname = triggerNames->at(bt); if(tname.find(lname)!=std::string::npos){ leptonOR = leptonOR || (triggerBits->at(bt)>0); float_values["Event_passes"+lname]=triggerBits->at(bt); float_values["Event_prescale"+lname]=triggerPrescales->at(bt); } } } for(size_t ht =0; ht< hadronicTriggers.size();++ht){ string hname = hadronicTriggers.at(ht); for(size_t bt = 0; bt < triggerNames->size();++bt){ std::string tname = triggerNames->at(bt); if(tname.find(hname)!=std::string::npos){ hadronOR = hadronOR || (triggerBits->at(bt)>0); float_values["Event_passes"+hname]=triggerBits->at(bt); float_values["Event_prescale"+hname]=triggerPrescales->at(bt); } } } float_values["Event_passesLeptonicTriggers"]=(float)leptonOR; float_values["Event_passesHadronicTriggers"]=(float)hadronOR; return (leptonOR || hadronOR); } void DMAnalysisTreeMaker::getPUSF(){ puWeight=(float) LumiWeights_.weight(nTruePU); puWeightUp = (float) LumiWeightsUp_.weight(nTruePU); puWeightDown = (float) LumiWeightsDown_.weight(nTruePU); // std::cout<<"nTruePU: "<pdf()->id.first == 21 || genprod->pdf()->id.second == 21) return; std::cout << " getting pdf "<pdf()->scalePDF; double x1 = genprod->pdf()->x.first; double x2 = genprod->pdf()->x.second; int id1 = genprod->pdf()->id.first; int id2 = genprod->pdf()->id.second; std::cout << " maxpdf "<< maxPdf << " accessing x1 " << x1<< " id1 " << id1< maxPDFCount ) continue; LHAPDF::usePDFMember(2, p); double xpdf1_new = LHAPDF::xfx(2, x1, scalePDF, id1); double xpdf2_new = LHAPDF::xfx(2, x2, scalePDF, id2); double pweight = xpdf1_new * xpdf2_new / w0; stringstream w_n; w_n << p; cout << "xpdf1 new " << xpdf1_new << endl; cout << "xpdf2 new " << xpdf2_new << endl; cout << "index " << p << " pweight " << pweight << endl; float_values["PDFWeight"+w_n.str()]= pweight; } } void DMAnalysisTreeMaker::getEventLHEWeights(){ // std::cout << " in weight "<weights().size(); // std::cout << "weight size "<< wgtsize<weights().at(i).wgt; // cout << "ww # " << i<< "is "<0) smear = std::max((double)(0.0), (double)(ptCorr + (ptCorr - genpt) * resolScale) / ptCorr); // return ptCorr * smear; // } double DMAnalysisTreeMaker::smear(double pt, double genpt, double eta, string syst){ double resolScale = resolSF(fabs(eta), syst); double smear =1.0; if(genpt>0) smear = std::max((double)(0.0), (double)(pt + (pt - genpt) * resolScale) / pt); return smear; } /* // 8 TeV scale factors double DMAnalysisTreeMaker::resolSF(double eta, string syst) { double fac = 0.; if (syst == "jer__up")fac = 1.; if (syst == "jer__down")fac = -1.; if (eta <= 0.5) return 0.079 + (0.026 * fac); else if ( eta > 0.5 && eta <= 1.1 ) return 0.099 + (0.028 * fac); else if ( eta > 1.1 && eta <= 1.7 ) return 0.121 + (0.029 * fac); else if ( eta > 1.7 && eta <= 2.3 ) return 0.208 + (0.046 * fac); else if ( eta > 2.3 && eta <= 2.8 ) return 0.254 + (0.062 * fac); else if ( eta > 2.8 && eta <= 3.2 ) return 0.395 + (0.063 * fac); else if ( eta > 3.2 && eta <= 5.0 ) return 0.056 + (0.191 * fac); return 0.1; } */ // preliminary 13 TeV scale factors double DMAnalysisTreeMaker::resolSF(double eta, string syst) { double fac = 0.; if (syst == "jer__up")fac = 1.; if (syst == "jer__down")fac = -1.; if (eta <= 0.8) return 0.061 + (0.023 * fac); else if ( eta > 0.8 && eta <= 1.3 ) return 0.088 + (0.029 * fac); else if ( eta > 1.3 && eta <= 1.9 ) return 0.106 + (0.030 * fac); else if ( eta > 1.9 && eta <= 2.5 ) return 0.126 + (0.094 * fac); else if ( eta > 2.5 && eta <= 3.0 ) return 0.343 + (0.123 * fac); else if ( eta > 3.0 && eta <= 3.2 ) return 0.303 + (0.111 * fac); else if ( eta > 3.2 && eta <= 5.0 ) return 0.320 + (0.286 * fac); return 0.1; } double DMAnalysisTreeMaker::jetUncertainty(double ptCorr, double eta, string syst) { if(ptCorr<0)return ptCorr; if(syst == "jes__up" || syst == "jes__down"){ double fac = 1.; if (syst == "jes__down")fac = -1.; jecUnc->setJetEta(eta); jecUnc->setJetPt(ptCorr); double JetCorrection = jecUnc->getUncertainty(true); return JetCorrection*fac; } return 0.0; } double DMAnalysisTreeMaker::getScaleFactor(double ptCorr,double etaCorr,double partonFlavour, string syst){ return 1.0; } bool DMAnalysisTreeMaker::isInVector(std::vector v, std::string s){ for(size_t i = 0;i= minTags && t <= maxTags); } float DMAnalysisTreeMaker::BTagWeight::weight(vector jetTags, int tags) { if (!filter(tags)) { // std::cout << "nThis event should not pass the selection, what is it doing here?" << std::endl; return 0; } int njetTags = jetTags.size(); //cout<< " njettags "<< njetTags<> j) & 0x1) == 1; if (tagged) { ntagged++; mc *= jetTags[j].eff; data *= jetTags[j].eff * jetTags[j].sf; } else { mc *= (1. - jetTags[j].eff); data *= (1. - jetTags[j].eff * jetTags[j].sf); } } if (filter(ntagged)) { //std::cout << mc << " " << data << endl; pMC += mc; pData += data; } } if (pMC == 0) return 0; return pData / pMC; } double DMAnalysisTreeMaker::MCTagEfficiency(string algo, int flavor, double pt){ if (abs(flavor) ==5){ if(algo=="csvt") return 0.38; if(algo=="csvm") return 0.58; if(algo=="csvl") return 0.755; } if (abs(flavor) ==4){ if(algo=="csvt") return 0.015; if(algo=="csvm") return 0.08; if(algo=="csvl") return 0.28; } if (abs(flavor) !=4 && abs(flavor) !=5){ if(algo=="csvt") return 0.0008; if(algo=="csvm") return 0.007; if(algo=="csvl") return 0.079; } return 1.0; } double DMAnalysisTreeMaker::TagScaleFactor(string algo, int flavor, string syst, double pt){ // source (02/11): // https://twiki.cern.ch/twiki/pub/CMS/BtagRecommendation76X/CSVv2_prelim.csv if(algo == "csvt"){ if(syst == "noSyst") { if(abs(flavor)==5){ if (pt >= 30 && pt < 670) return 0.886376*((1.+(0.00250226*pt))/(1.+(0.00193725*pt))); } if(abs(flavor)==4){ if (pt >= 30 && pt < 670) return 0.886376*((1.+(0.00250226*pt))/(1.+(0.00193725*pt))); } if(abs(flavor)!=5 && abs(flavor)!=4){ return 0.992339; } } if(syst == "mistag_up") { if(abs(flavor)==5){ return 1.00; } if(abs(flavor)==4){ return 1.00; } if(abs(flavor)!=5 && abs(flavor)!=4){ return 1.17457; } } if(syst == "mistag_down") { if(abs(flavor)==5){ return 1.00; } if(abs(flavor)==4){ return 1.00; } if(abs(flavor)!=5 && abs(flavor)!=4){ return 0.810103; } } if(syst == "b_tag_up") { if(abs(flavor)==5){ if (pt >= 30 && pt < 50 ) return (0.886376*((1.+(0.00250226*pt))/(1.+(0.00193725*pt))))+0.019803794100880623; if (pt >= 50 && pt < 70 ) return (0.886376*((1.+(0.00250226*pt))/(1.+(0.00193725*pt))))+0.026958625763654709; if (pt >= 70 && pt < 100) return (0.886376*((1.+(0.00250226*pt))/(1.+(0.00193725*pt))))+0.024285079911351204; if (pt >= 100 && pt < 140) return (0.886376*((1.+(0.00250226*pt))/(1.+(0.00193725*pt))))+0.028512096032500267; if (pt >= 140 && pt < 200) return (0.886376*((1.+(0.00250226*pt))/(1.+(0.00193725*pt))))+0.029808893799781799; if (pt >= 200 && pt < 300) return (0.886376*((1.+(0.00250226*pt))/(1.+(0.00193725*pt))))+0.026503190398216248; if (pt >= 300 && pt < 670) return (0.886376*((1.+(0.00250226*pt))/(1.+(0.00193725*pt))))+0.042264193296432495 ; } if(abs(flavor)==4){ if (pt >= 30 && pt < 50 ) return (0.886376*((1.+(0.00250226*pt))/(1.+(0.00193725*pt))))+0.039607588201761246; if (pt >= 50 && pt < 70 ) return (0.886376*((1.+(0.00250226*pt))/(1.+(0.00193725*pt))))+0.053917251527309418; if (pt >= 70 && pt < 100) return (0.886376*((1.+(0.00250226*pt))/(1.+(0.00193725*pt))))+0.048570159822702408; if (pt >= 100 && pt < 140) return (0.886376*((1.+(0.00250226*pt))/(1.+(0.00193725*pt))))+0.057024192065000534; if (pt >= 140 && pt < 200) return (0.886376*((1.+(0.00250226*pt))/(1.+(0.00193725*pt))))+0.059617787599563599; if (pt >= 200 && pt < 300) return (0.886376*((1.+(0.00250226*pt))/(1.+(0.00193725*pt))))+0.053006380796432495; if (pt >= 300 && pt < 670) return (0.886376*((1.+(0.00250226*pt))/(1.+(0.00193725*pt))))+0.08452838659286499; } if(abs(flavor)!=5 && abs(flavor)!=4){ return 1.0; } } if(syst == "b_tag_down") { if(abs(flavor)==5){ if (pt >= 30 && pt < 50 ) return (0.886376*((1.+(0.00250226*pt))/(1.+(0.00193725*pt))))-0.019803794100880623; if (pt >= 50 && pt < 70 ) return (0.886376*((1.+(0.00250226*pt))/(1.+(0.00193725*pt))))-0.026958625763654709; if (pt >= 70 && pt < 100) return (0.886376*((1.+(0.00250226*pt))/(1.+(0.00193725*pt))))-0.024285079911351204; if (pt >= 100 && pt < 140) return (0.886376*((1.+(0.00250226*pt))/(1.+(0.00193725*pt))))-0.028512096032500267; if (pt >= 140 && pt < 200) return (0.886376*((1.+(0.00250226*pt))/(1.+(0.00193725*pt))))-0.029808893799781799; if (pt >= 200 && pt < 300) return (0.886376*((1.+(0.00250226*pt))/(1.+(0.00193725*pt))))-0.026503190398216248; if (pt >= 300 && pt < 670) return (0.886376*((1.+(0.00250226*pt))/(1.+(0.00193725*pt))))-0.042264193296432495; } if(abs(flavor)==4){ if (pt >= 30 && pt < 50 ) return (0.886376*((1.+(0.00250226*pt))/(1.+(0.00193725*pt))))-0.039607588201761246; if (pt >= 50 && pt < 70 ) return (0.886376*((1.+(0.00250226*pt))/(1.+(0.00193725*pt))))-0.053917251527309418; if (pt >= 70 && pt < 100) return (0.886376*((1.+(0.00250226*pt))/(1.+(0.00193725*pt))))-0.048570159822702408; if (pt >= 100 && pt < 140) return (0.886376*((1.+(0.00250226*pt))/(1.+(0.00193725*pt))))-0.057024192065000534; if (pt >= 140 && pt < 200) return (0.886376*((1.+(0.00250226*pt))/(1.+(0.00193725*pt))))-0.059617787599563599; if (pt >= 200 && pt < 300) return (0.886376*((1.+(0.00250226*pt))/(1.+(0.00193725*pt))))-0.053006380796432495; if (pt >= 300 && pt < 670) return (0.886376*((1.+(0.00250226*pt))/(1.+(0.00193725*pt))))-0.08452838659286499; } if(abs(flavor)!=5 && abs(flavor)!=4){ return 1.0; } } } return 1.0; } float DMAnalysisTreeMaker::BTagWeight::weightWithVeto(vector jetsTags, int tags, vector jetsVetoes, int vetoes) {//This function takes into account cases where you have n b-tags and m vetoes, but they have different thresholds. if (!filter(tags)) { // std::cout << "nThis event should not pass the selection, what is it doing here?" << std::endl; return 0; } int njets = jetsTags.size(); if(njets != (int)(jetsVetoes.size()))return 0;//jets tags and vetoes must have same size! int comb = 1 << njets; float pMC = 0; float pData = 0; for (int i = 0; i < comb; i++) { float mc = 1.; float data = 1.; int ntagged = 0; for (int j = 0; j < njets; j++) { bool tagged = ((i >> j) & 0x1) == 1; if (tagged) { ntagged++; mc *= jetsTags[j].eff; data *= jetsTags[j].eff * jetsTags[j].sf; } else { mc *= (1. - jetsVetoes[j].eff); data *= (1. - jetsVetoes[j].eff * jetsVetoes[j].sf); } } if (filter(ntagged)) { // std::cout << mc << " " << data << endl; pMC += mc; pData += data; } } if (pMC == 0) return 0; return pData / pMC; } // muon SF based on Muon POG values from Dec15 double DMAnalysisTreeMaker::muonSF(bool isdata, float pt, float eta, int syst) { if (isdata) return 1; double eff_id = 0; double eff_iso = 0; double eff_trigger = 0; // tight id SF if (pt > 20 && pt <= 25) { if (abs(eta) > 0 && abs(eta) <= 0.9) eff_id = 0.9752116203308105 + 0.0030660638813280626 * syst; if (abs(eta) > 0.9 && abs(eta) <= 1.2) eff_id = 0.9738101959228516 + 0.004502934246978295 * syst; if (abs(eta) > 1.2 && abs(eta) <= 2.1) eff_id = 0.9983288645744324 + 0.002331323348626783 * syst; if (abs(eta) > 2.1 && abs(eta) <= 2.4) eff_id = 0.9877836108207703 + 0.004915740433340289 * syst; } if (pt > 25 && pt <= 30) { if (abs(eta) > 0 && abs(eta) <= 0.9) eff_id = 0.9848297238349915 + 0.0016307213764927449 * syst; if (abs(eta) > 0.9 && abs(eta) <= 1.2) eff_id = 0.978645384311676 + 0.0027064755458685794 * syst; if (abs(eta) > 1.2 && abs(eta) <= 2.1) eff_id = 0.9905462265014648 + 0.001402578599690647 * syst; if (abs(eta) > 2.1 && abs(eta) <= 2.4) eff_id = 0.9802553653717041 + 0.003173276637083633 * syst; } if (pt > 30 && pt <= 40) { if (abs(eta) > 0 && abs(eta) <= 0.9) eff_id = 0.9861794114112854 + 0.0006187187412138267 * syst; if (abs(eta) > 0.9 && abs(eta) <= 1.2) eff_id = 0.9798933267593384 + 0.001057081371390319 * syst; if (abs(eta) > 1.2 && abs(eta) <= 2.1) eff_id = 0.9923668503761292 + 0.0005653311393042486 * syst; if (abs(eta) > 2.1 && abs(eta) <= 2.4) eff_id = 0.9785045385360718 + 0.0015542030446523895 * syst; } if (pt > 40 && pt <= 50) { if (abs(eta) > 0 && abs(eta) <= 0.9) eff_id = 0.987443208694458 + 0.000494159746725046 * syst; if (abs(eta) > 0.9 && abs(eta) <= 1.2) eff_id = 0.980233907699585 + 0.000819615406448897 * syst; if (abs(eta) > 1.2 && abs(eta) <= 2.1) eff_id = 0.9927627444267273 + 0.0004155573807947332 * syst; if (abs(eta) > 2.1 && abs(eta) <= 2.4) eff_id = 0.9778544902801514 + 0.001456799997296391 * syst; } if (pt > 50 && pt <= 60) { if (abs(eta) > 0 && abs(eta) <= 0.9) eff_id = 0.9834294319152832 + 0.0011818999573518245 * syst; if (abs(eta) > 0.9 && abs(eta) <= 1.2) eff_id = 0.9773300886154175 + 0.001955436343316424 * syst; if (abs(eta) > 1.2 && abs(eta) <= 2.1) eff_id = 0.9886322021484375 + 0.0011254961157344963 * syst; if (abs(eta) > 2.1 && abs(eta) <= 2.4) eff_id = 0.9654409885406494 + 0.003709169009223743 * syst; } if (pt > 60 && pt <= 120) { if (abs(eta) > 0 && abs(eta) <= 0.9) eff_id = 0.9863178730010986 + 0.002073330940717176 * syst; if (abs(eta) > 0.9 && abs(eta) <= 1.2) eff_id = 0.9795225858688354 + 0.0035622593553725837 * syst; if (abs(eta) > 1.2 && abs(eta) <= 2.1) eff_id = 0.9950451850891113 + 0.002673833447209764 * syst; if (abs(eta) > 2.1 && abs(eta) <= 2.4) eff_id = 0.9689615368843079 + 0.011084748199568817 * syst; } if (abs(eta) > 2.4 || pt > 120) eff_id = 1; // additional 1% syst. uncertainty eff_id *= 1 + (0.01 * syst); // tight iso SF if (pt > 20 && pt <= 25) { if (abs(eta) > 0 && abs(eta) <= 0.9) eff_iso = 1.0043761730194092 + 0.003959090391076143 * syst; if (abs(eta) > 0.9 && abs(eta) <= 1.2) eff_iso = 1.004357933998108 + 0.006125539530136138 * syst; if (abs(eta) > 1.2 && abs(eta) <= 2.1) eff_iso = 0.9970762133598328 + 0.003109125287470401 * syst; if (abs(eta) > 2.1 && abs(eta) <= 2.4) eff_iso = 0.9957730770111084 + 0.006137193387970902 * syst; } if (pt > 25 && pt <= 30) { if (abs(eta) > 0 && abs(eta) <= 0.9) eff_iso = 0.9995378255844116 + 0.0022512071035640673 * syst; if (abs(eta) > 0.9 && abs(eta) <= 1.2) eff_iso = 1.002331256866455 + 0.004003683572512011 * syst; if (abs(eta) > 1.2 && abs(eta) <= 2.1) eff_iso = 1.0006532669067383 + 0.002067755362435184 * syst; if (abs(eta) > 2.1 && abs(eta) <= 2.4) eff_iso = 0.9939026832580566 + 0.004261971076013437 * syst; } if (pt > 30 && pt <= 40) { if (abs(eta) > 0 && abs(eta) <= 0.9) eff_iso = 1.000901222229004 + 0.0007979481788689052 * syst; if (abs(eta) > 0.9 && abs(eta) <= 1.2) eff_iso = 1.004658579826355 + 0.0014502638048416372 * syst; if (abs(eta) > 1.2 && abs(eta) <= 2.1) eff_iso = 1.0023553371429443 + 0.0008445520691793605 * syst; if (abs(eta) > 2.1 && abs(eta) <= 2.4) eff_iso = 0.997478187084198 + 0.001781225374381486 * syst; } if (pt > 40 && pt <= 50) { if (abs(eta) > 0 && abs(eta) <= 0.9) eff_iso = 0.9986253976821899 + 0.0004518361024064332 * syst; if (abs(eta) > 0.9 && abs(eta) <= 1.2) eff_iso = 1.0013608932495117 + 0.0004888604573095644 * syst; if (abs(eta) > 1.2 && abs(eta) <= 2.1) eff_iso = 0.999933660030365 + 0.0004309914887707696 * syst; if (abs(eta) > 2.1 && abs(eta) <= 2.4) eff_iso = 1.002805233001709 + 0.001100242856214239 * syst; } if (pt > 50 && pt <= 60) { if (abs(eta) > 0 && abs(eta) <= 0.9) eff_iso = 1.0002487897872925 + 0.000772847340102783 * syst; if (abs(eta) > 0.9 && abs(eta) <= 1.2) eff_iso = 0.9986217021942139 + 0.0012396364566794034 * syst; if (abs(eta) > 1.2 && abs(eta) <= 2.1) eff_iso = 1.0002963542938232 + 0.0007614160360063238 * syst; if (abs(eta) > 2.1 && abs(eta) <= 2.4) eff_iso = 1.0043764114379883 + 0.001806526581100641 * syst; } if (pt > 60 && pt <= 120) { if (abs(eta) > 0 && abs(eta) <= 0.9) eff_iso = 0.9986850023269653 + 0.0008907575174433545 * syst; if (abs(eta) > 0.9 && abs(eta) <= 1.2) eff_iso = 1.0054655075073242 + 0.001589130019220112 * syst; if (abs(eta) > 1.2 && abs(eta) <= 2.1) eff_iso = 1.0004935264587402 + 0.0009382223143922724 * syst; if (abs(eta) > 2.1 && abs(eta) <= 2.4) eff_iso = 1.0010104179382324 + 0.0022795762936220253 * syst; } if (abs(eta) > 2.4 || pt > 120) eff_iso = 1; // additional 1% syst. uncertainty eff_iso *= 1 + (0.01 * syst); // IsoMu20 trigger SF /* This is not trivial as there are two sets of SF depending on the HLT menu. Values are the weighted average of both sets: total lumi: 2197.95 pb-1 HLT v4.2: run <= 257819 393.47 pb-1 0.179 HLT v4.3: run > 257819 1804.48 pb-1 0.821 */ if (pt > 22 && pt <= 25) { if (abs(eta) > 0 && abs(eta) <= 0.9) eff_trigger = 0.987956 + 0.0044549 * syst; if (abs(eta) > 0.9 && abs(eta) <= 1.2) eff_trigger = 1.0178482 + 0.0086129 * syst; if (abs(eta) > 1.2 && abs(eta) <= 2.1) eff_trigger = 0.9971603 + 0.0047654 * syst; if (abs(eta) > 2.1 && abs(eta) <= 2.4) eff_trigger = 1.0132546 + 0.0127635 * syst; } if (pt > 25 && pt <= 30) { if (abs(eta) > 0 && abs(eta) <= 0.9) eff_trigger = 0.9953616 + 0.0023685 * syst; if (abs(eta) > 0.9 && abs(eta) <= 1.2) eff_trigger = 1.0096846 + 0.0048012 * syst; if (abs(eta) > 1.2 && abs(eta) <= 2.1) eff_trigger = 0.9969169 + 0.0028327 * syst; if (abs(eta) > 2.1 && abs(eta) <= 2.4) eff_trigger = 1.0136236 + 0.0073907 * syst; } if (pt > 30 && pt <= 40) { if (abs(eta) > 0 && abs(eta) <= 0.9) eff_trigger = 0.9939314 + 0.0003506 * syst; if (abs(eta) > 0.9 && abs(eta) <= 1.2) eff_trigger = 1.0001498 + 0.0020148 * syst; if (abs(eta) > 1.2 && abs(eta) <= 2.1) eff_trigger = 0.98972 + 0.0013253 * syst; if (abs(eta) > 2.1 && abs(eta) <= 2.4) eff_trigger = 1.0178737 + 0.003679 * syst; } if (pt > 40 && pt <= 50) { if (abs(eta) > 0 && abs(eta) <= 0.9) eff_trigger = 0.9952441 + 0.0033988 * syst; if (abs(eta) > 0.9 && abs(eta) <= 1.2) eff_trigger = 0.9939897 + 0.0015432 * syst; if (abs(eta) > 1.2 && abs(eta) <= 2.1) eff_trigger = 0.9892971 + 0.000579 * syst; if (abs(eta) > 2.1 && abs(eta) <= 2.4) eff_trigger = 1.0114995 + 0.0030685 * syst; } if (pt > 50 && pt <= 60) { if (abs(eta) > 0 && abs(eta) <= 0.9) eff_trigger = 0.9939026 + 0.0014432 * syst; if (abs(eta) > 0.9 && abs(eta) <= 1.2) eff_trigger = 0.9977205 + 0.0243324 * syst; if (abs(eta) > 1.2 && abs(eta) <= 2.1) eff_trigger = 0.9889162 + 0.001979 * syst; if (abs(eta) > 2.1 && abs(eta) <= 2.4) eff_trigger = 1.019157 + 0.004705 * syst; } if (pt > 60 && pt <= 120) { if (abs(eta) > 0 && abs(eta) <= 0.9) eff_trigger = 0.9878582 + 0.0018148 * syst; if (abs(eta) > 0.9 && abs(eta) <= 1.2) eff_trigger = 0.9833346 + 0.0039938 * syst; if (abs(eta) > 1.2 && abs(eta) <= 2.1) eff_trigger = 0.9912019 + 0.0026506 * syst; if (abs(eta) > 2.1 && abs(eta) <= 2.4) eff_trigger = 1.0122248 + 0.0082802 * syst; } if (abs(eta) > 2.4 || pt > 120) eff_trigger = 1; // additional 0.5% syst. uncertainty eff_trigger *= 1 + (0.005 * syst); // total SF double eff_total = eff_id * eff_iso * eff_trigger; return eff_total; } //DMAnalysisTreeMaker::~DMAnalysisTreeMaker(const edm::ParameterSet& iConfig) // ------------ method called once each job just after ending the event loop ------------ #include "FWCore/Framework/interface/MakerMacros.h" DEFINE_FWK_MODULE(DMAnalysisTreeMaker);