| 1 |
// #############################################
|
| 2 |
// Note : this is a snippet to re-weighting pdf.
|
| 3 |
// Not a stand along code (For Madgraph)
|
| 4 |
// #############################################
|
| 5 |
|
| 6 |
// Force unsetting of the LHAPATH variable to avoid potential problems
|
| 7 |
gSystem->Setenv("LHAPATH","<Your Local LHAPDFLIB Dir>/share/lhapdf/PDFsets/");
|
| 8 |
LHAPDF::setVerbosity(LHAPDF::SILENT);
|
| 9 |
LHAPDF::initPDFSet("cteq65", LHAPDF::LHGRID, 0); // Type: LHGRID(preferred) or LHPDF
|
| 10 |
LHAPDF::getDescription();
|
| 11 |
int nmembers = LHAPDF::numberPDF(1);
|
| 12 |
|
| 13 |
// ....
|
| 14 |
|
| 15 |
void top_iks::Loop(const TString outName, const TString jetName, const int flavor) {
|
| 16 |
|
| 17 |
// Get the collection
|
| 18 |
TopGenPdfInfo pdfstuff;
|
| 19 |
bool runPDF = false;
|
| 20 |
TopGenPdfInfo pdfstuff = ntuple->muons;
|
| 21 |
if (pdfstuff) runPDF = true;
|
| 22 |
|
| 23 |
// Reweighting
|
| 24 |
std::vector<double> weights;
|
| 25 |
if (debug) std::cout<<"runPDF before Top Reconstruction = "<<runPDF<<std::endl;
|
| 26 |
if(runPDF){
|
| 27 |
float Q = pdfstuff->pdf()->scalePDF;
|
| 28 |
int id1 = (int)pdfstuff->pdf()->id.first;
|
| 29 |
double x1 = pdfstuff->pdf()->x.first;
|
| 30 |
int id2 = (int)pdfstuff->pdf()->id.second;
|
| 31 |
double x2 = pdfstuff->pdf()->x.second;
|
| 32 |
|
| 33 |
LHAPDF::usePDFMember(1,0);
|
| 34 |
double pdf1 = LHAPDF::xfx(x1, Q, id1)/x1;
|
| 35 |
double pdf2 = LHAPDF::xfx(x2, Q, id2)/x2;
|
| 36 |
|
| 37 |
if(debug){
|
| 38 |
std::cout<<"Q = "<<Q<<std::endl;
|
| 39 |
std::cout<<"(id1,x1,pdf1) = "<<id1<<","<<x1<<","<<pdf1<<std::endl;
|
| 40 |
std::cout<<"(id2,x2,pdf2) = "<<id2<<","<<x2<<","<<pdf2<<std::endl;
|
| 41 |
}
|
| 42 |
|
| 43 |
if (pdf1 != 0 && pdf2 != 0){
|
| 44 |
// Put PDF weights in the event
|
| 45 |
for (unsigned int i=0; i<=nmembers; ++i) {//2n+1
|
| 46 |
LHAPDF::usePDFMember(1,i);
|
| 47 |
double newpdf1 = LHAPDF::xfx(x1, Q, id1)/x1;
|
| 48 |
double newpdf2 = LHAPDF::xfx(x2, Q, id2)/x2;
|
| 49 |
if(debug){
|
| 50 |
std::cout<<"newpdf1 = "<<newpdf1<<"; ";
|
| 51 |
std::cout<<"newpdf2 = "<<newpdf2<<"; ";
|
| 52 |
std::cout<<"weights["<<i<<"] = "<<newpdf1/pdf1*newpdf2/pdf2<<std::endl;
|
| 53 |
}
|
| 54 |
weights.push_back(newpdf1/pdf1*newpdf2/pdf2);
|
| 55 |
}
|
| 56 |
}else{
|
| 57 |
std::cout<<"weights diverge, skip this event!"<<std::endl;
|
| 58 |
runPDF = false;
|
| 59 |
}
|
| 60 |
}
|
| 61 |
|
| 62 |
|
| 63 |
// Do your event selections
|
| 64 |
|
| 65 |
|
| 66 |
hs["M3_3jet"]->Fill(hadTopP4.M());
|
| 67 |
// Example of how to apply the weight
|
| 68 |
if (runPDF){
|
| 69 |
for(int i=1; i <= nmembers; i++){//2n, not taking into account weights[0] since it's 1.
|
| 70 |
TString suffix = "PDF_";
|
| 71 |
suffix+=i;
|
| 72 |
double ww = weights[i];
|
| 73 |
if(debug){
|
| 74 |
std::cout<<"Weights = "<< ww <<"; ";
|
| 75 |
std::cout<<TString("HadronicTop_mass_M3_PDF_")+suffix<<std::endl;
|
| 76 |
}
|
| 77 |
hs[TString("M3_3jet")+suffix]->Fill(hadTopP4.M(),ww);
|
| 78 |
}
|
| 79 |
}
|
| 80 |
}
|
