Skip to content

Instantly share code, notes, and snippets.

@lucasc896
Last active August 29, 2015 14:19
Show Gist options
  • Select an option

  • Save lucasc896/ad37712e9c17c4bd7292 to your computer and use it in GitHub Desktop.

Select an option

Save lucasc896/ad37712e9c17c4bd7292 to your computer and use it in GitHub Desktop.
double get_jet_weight(double isr_pt){
// function to calculate an individual ISR jet's weight
double qcut = 44.;
if(isr_pt<qcut){
return (1. + (2.63e-5)*(isr_pt-qcut) + (1.2e-4)*pow((isr_pt-qcut), 2));
}
else if(isr_pt == qcut){
return 1.;
}
else{
return (1. + 0.44*(1 - exp((-3.34e-3)*(isr_pt - qcut))));
}
// should never get here...!
return 1.;
}
double get_jetsum_weight(double jet_sum, double jet1, double jet2){
//function to calculate the vectorially summed jets weight
double qcut = 44.;
if((jet1 < qcut) || (jet2 < qcut)){
return 1.;
}
//calculate the min value: min(pT, 300)
double ptval = jet_sum < 300. ? jet_sum : 300.;
return (1. - 0.02 + (-4.e-3)*ptval + (7.e-5)*pow(ptval, 2));
}
bool Apply(double * weight){
PolarLorentzV genPtVect(0.,0.,0.,0.);
std::vector<Event::GenObject> gISR;
std::vector<int> stopMums;
// first get stop mother IDs
for( std::vector<Event::GenObject>::const_iterator genObj = mEv->GenParticles().begin(); genObj != mEv->GenParticles().end(); ++genObj ){
if (genObj->GetIndex() > 30) break;
if (genObj->GetStatus() != 3) continue;
if (fabs(genObj->GetID()) == 1000006) stopMums.push_back((*genObj).GetMotherID());
if(stopMums.size() == 2) break; // break once two stop mums have been found
}
// loop again to find all gen particles with mother's the same ID as the stop's mothers
for( std::vector<Event::GenObject>::const_iterator genObj = mEv->GenParticles().begin(); genObj != mEv->GenParticles().end(); ++genObj ){
if (genObj->GetIndex() > 30) break;
if (genObj->GetStatus() != 3) continue;
// don't consider the stops
if (fabs(genObj->GetID()) == 1000006) continue;
// if a gen particle has the same mother as one of the stop's, store it
if(std::find(stopMums.begin(), stopMums.end(), genObj->GetMotherID()) != stopMums.end()){
gISR.push_back((*genObj));
}
}
// don't apply weight to events with other than two ISR particles
if(gISR.size() != 2) return true;
// loop over found particles and get Vectorial Pt Sum
for(UInt_t i=0; i<gISR.size(); i++){
genPtVect += gISR.at(i);
}
// calculate weights according to:
// https://www.dropbox.com/s/gk0agjbtryg669a/15-01-tbt-mikulec.pdf?dl=0
double w_pt1 = get_jet_weight(gISR.at(0).Pt());
double w_pt2 = get_jet_weight(gISR.at(1).Pt());
double w_sum = get_jetsum_weight(genPtVect.Pt(), gISR.at(0).Pt(), gISR.at(1).Pt());
// apply weights to the global weight (*weight)
*weight *= (w_pt1 * w_pt2 * w_sum);
return true;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment