| 1 |
jengbou |
1.1 |
// #############################################
|
| 2 |
|
|
// Note : this is a snippet to re-weighting pdf.
|
| 3 |
|
|
// Not a stand along code
|
| 4 |
|
|
// #############################################
|
| 5 |
|
|
|
| 6 |
|
|
|
| 7 |
|
|
PDFAnalyzer::PDFAnalyzer(const edm::ParameterSet& iConfig)
|
| 8 |
|
|
{
|
| 9 |
|
|
|
| 10 |
|
|
/// Parameters:
|
| 11 |
|
|
pdfInfoTag = iConfig.getUntrackedParameter<edm::InputTag> ("PdfInfoTag", edm::InputTag("genEventPdfInfo"));
|
| 12 |
|
|
pdfSetName = iConfig.getUntrackedParameter<std::string> ("PdfSetName", "cteq65");
|
| 13 |
|
|
|
| 14 |
|
|
|
| 15 |
|
|
// Force unsetting of the LHAPATH variable to avoid potential problems
|
| 16 |
|
|
gSystem->Setenv("LHAPATH","<Your Local LHAPDFLIB Dir>/share/lhapdf/PDFsets/");
|
| 17 |
|
|
//gSystem->Setenv("LHAPATH","/uscmst1/prod/sw/cms/slc5_ia32_gcc434/external/lhapdf/5.6.0-cms2/share/lhapdf/PDFsets/");
|
| 18 |
|
|
//gSystem->Setenv("LHAPATH","/afs/cern.ch/cms/sw/slc4_ia32_gcc345/external/lhapdf/5.6.0-cms/share/lhapdf/PDFsets/");
|
| 19 |
|
|
LHAPDF::setVerbosity(LHAPDF::SILENT);
|
| 20 |
|
|
LHAPDF::initPDFSet(pdfSetName, LHAPDF::LHGRID, 0); // Type: LHGRID(preferred) or LHPDF
|
| 21 |
|
|
LHAPDF::getDescription();
|
| 22 |
|
|
|
| 23 |
|
|
nmembers = LHAPDF::numberPDF() + 1;
|
| 24 |
|
|
}
|
| 25 |
|
|
|
| 26 |
|
|
// ....
|
| 27 |
|
|
|
| 28 |
|
|
void
|
| 29 |
|
|
PDFAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
|
| 30 |
|
|
{
|
| 31 |
|
|
|
| 32 |
|
|
// Get the collection
|
| 33 |
|
|
edm::Handle<reco::PdfInfo> pdfstuff;
|
| 34 |
|
|
bool runPDF = false;
|
| 35 |
|
|
if ((!iEvent.isRealData()) && iEvent.getByLabel(pdfInfoTag, pdfstuff))
|
| 36 |
|
|
runPDF = true;
|
| 37 |
|
|
|
| 38 |
|
|
|
| 39 |
|
|
// Reweighting
|
| 40 |
|
|
std::vector<double> weights;
|
| 41 |
|
|
if (debug) std::cout<<"runPDF before Top Reconstruction = "<<runPDF<<std::endl;
|
| 42 |
|
|
if(runPDF){
|
| 43 |
|
|
float Q = pdfstuff->scalePDF;
|
| 44 |
|
|
int id1 = (int)pdfstuff->id1;
|
| 45 |
|
|
double x1 = pdfstuff->x1;
|
| 46 |
|
|
double pdf1 = pdfstuff->pdf1;
|
| 47 |
|
|
int id2 = (int)pdfstuff->id2;
|
| 48 |
|
|
double x2 = pdfstuff->x2;
|
| 49 |
|
|
double pdf2 = pdfstuff->pdf2;
|
| 50 |
|
|
if(debug){
|
| 51 |
|
|
std::cout<<"Q = "<<Q<<std::endl;
|
| 52 |
|
|
std::cout<<"(id1,x1,pdf1) = "<<id1<<","<<x1<<","<<pdf1<<std::endl;
|
| 53 |
|
|
std::cout<<"(id2,x2,pdf2) = "<<id2<<","<<x2<<","<<pdf2<<std::endl;
|
| 54 |
|
|
}
|
| 55 |
|
|
|
| 56 |
|
|
if (pdf1 != 0 && pdf2 != 0){
|
| 57 |
|
|
// Put PDF weights in the event
|
| 58 |
|
|
weights.reserve(nmembers);
|
| 59 |
|
|
for (unsigned int i=0; i<nmembers; ++i) {
|
| 60 |
|
|
LHAPDF::usePDFMember(i);
|
| 61 |
|
|
double newpdf1 = LHAPDF::xfx(x1, Q, id1)/x1;
|
| 62 |
|
|
double newpdf2 = LHAPDF::xfx(x2, Q, id2)/x2;
|
| 63 |
|
|
if(debug){
|
| 64 |
|
|
std::cout<<"newpdf1 = "<<newpdf1<<"; ";
|
| 65 |
|
|
std::cout<<"newpdf2 = "<<newpdf2<<"; ";
|
| 66 |
|
|
std::cout<<"weights["<<i<<"] = "<<newpdf1/pdf1*newpdf2/pdf2<<std::endl;
|
| 67 |
|
|
}
|
| 68 |
|
|
weights.push_back(newpdf1/pdf1*newpdf2/pdf2);
|
| 69 |
|
|
}
|
| 70 |
|
|
}else{
|
| 71 |
|
|
std::cout<<"weights diverge, skip this event!"<<std::endl;
|
| 72 |
|
|
runPDF = false;
|
| 73 |
|
|
}
|
| 74 |
|
|
}
|
| 75 |
|
|
|
| 76 |
|
|
|
| 77 |
|
|
// Do your event selections
|
| 78 |
|
|
|
| 79 |
|
|
|
| 80 |
|
|
// Example of how to apply the weight
|
| 81 |
|
|
if (runPDF){
|
| 82 |
|
|
for(int i=0; i < nmembers; i++){
|
| 83 |
|
|
TString suffix = "PDF_";
|
| 84 |
|
|
suffix+=i;
|
| 85 |
|
|
double ww = weights[i]/weights[0];
|
| 86 |
|
|
if(debug){
|
| 87 |
|
|
std::cout<<"Weights = "<< ww <<"; ";
|
| 88 |
|
|
std::cout<<TString("HadronicTop_mass_M3_")+suffix<<std::endl;
|
| 89 |
|
|
}
|
| 90 |
|
|
h1[TString("HadronicTop_mass_M3_")+suffix]->Fill(hadTopP4.M(),ww);
|
| 91 |
|
|
}
|
| 92 |
|
|
}
|
| 93 |
|
|
}
|
