Three programs/utilities are discussed below: -------------------------------------------- //############################################################################ CASE 1: YOU WANT TO DETERMINE THE X1 AND X2 VALUES USED IN PYTHIA IN THE GENERATION STEP //############################################################################ - Disclaimer: I have not checked the result against the true KINE generated values. Provided that no bugs are present, both should give the same answer, apart from negligible ad-hoc Pt components, introduced in Pythia to "simulate" relative parton momenta inside the proton. - The idea is to reconstruct the "backward showering" logic implemented in Pythia. This is described in the section "INITIAL STATE SHOWERS", subsection "Transverse evolution" of the manual. There are successive rotations and boosts applied to create additional partons in the initial state. At the end, both protons are transformed back into the true center-of-mass frame of the detector, with the "domino" effect of changing the original Pz/E values used in the hard collision and creating non-zero Pt components. There is also some Q2 ordering involved in the backward evolution. - In most cases, the calculation is not far from the naive x=abs(Pz)/E expectation. It may differ by factors of order two only when the final generated PT of the parton is very high. The Q2 value for the PDF has to be known a priori. In qqbar->Z, it is evident that it is MZ**2. - Below you will find a "method" to be used in ExRootAnalysis. However, it can be easily adapted to any analysis structure if you have access to all Pythia partons up to the hard scattering process. "gen1" and "gen2" are the links to the two partons just before the hard scattering. The usage is illustrated for the qqbar->Z case. The relevant "method" is "FindPDFVariables". //>>>>>>>>>> BEGIN FILE >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> // // The example should run as: // $ root.exe // > gSystem->Load("libExRootAnalysisReader"); // > .x PDFexample.C // bool FindPDFVariables(TClonesArray* branchGen , TRootGenParticle* gen1, TRootGenParticle* gen2); void PDFexample(){ TChain cms("Analysis"); cms.Add("sig/test_*.root"); ExRootTreeReader* tree = new ExRootTreeReader(&cms); //printf("Entries in stream: %d\n", tree->GetEntries()); TClonesArray* branchGen = tree->UseBranch("Gen"); for (int entry=0; entry<10; entry++) { tree->ReadEntry(entry); int ngen = branchGen->GetEntries(); TRootGenParticle* gen1 = 0; TRootGenParticle* gen2 = 0; for (int ig=0; igAt(ig); if (gen->PID!=23) continue; if (gen->M1<0) continue; if (gen->M2<0) continue; TRootGenParticle* gtmp1 = (TRootGenParticle*)branchGen->At(gen->M1); if (gtmp1->PID==23) continue; TRootGenParticle* gtmp2 = (TRootGenParticle*)branchGen->At(gen->M2); if (gtmp2->PID==23) continue; gen1 = gtmp1; gen2 = gtmp2; break; } if (gen1<=0 || gen2<=0) return -1; TArrayD* x = FindPDFVariables(branchGen,gen1,gen2); if (x) { TRootGenParticle* gbeam = (TRootGenParticle*)branchGen->At(0); printf("\nPt/E values >>> Pt1/E1 %.6f Pt2/E1 %.6f\n" , fabs(gen1->PT)/gen1->E , fabs(gen2->PT)/gen2->E); printf("NAIVE x values >>> x1 %.6f x2 %.6f\n" , fabs(gen1->Pz)/gbeam->E , fabs(gen2->Pz)/gbeam->E); printf("TRUE x values >>> x1 %.6f x2 %.6f\n" , x->At(0), x->At(1)); } } } //>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> double Mass2(TRootGenParticle* ge1,TRootGenParticle* ge2){ double E = ge1->E + ge2->E; double Px = ge1->Px + ge2->Px; double Py = ge1->Py + ge2->Py; double Pz = ge1->Pz + ge2->Pz; return E*E - Px*Px - Py*Py - Pz*Pz; } //>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> double Q2(TRootGenParticle* ge1){ double E = ge1->E; double Px = ge1->Px; double Py = ge1->Py; double Pz = ge1->Pz; return fabs(E*E - Px*Px - Py*Py - Pz*Pz); } //>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> TArrayD* FindPDFVariables(TClonesArray* branchGen, TRootGenParticle* gen1, TRootGenParticle* gen2){ // gen1 and gen2 are the links to the partons // before the hard scattering process if (!gen1) return NULL; if (!gen2) return NULL; // gbeam1 and gbeam2 are the links to the two initial protons TRootGenParticle* gbeam1 = branchGen->At(0); TRootGenParticle* gbeam2 = branchGen->At(1); static TArrayD x(2); x[0] = 1.; x[1] = 1.; TRootGenParticle* g1 = gen1; TRootGenParticle* g2 = gen2; TRootGenParticle* gprev1 = NULL; TRootGenParticle* gprev2 = NULL; while (1) { if (g1==gbeam1 && g2==gbeam2) break; if (g1==gprev1 && g2==gprev2) { printf("WHAT? Backward link lost!\n"); return NULL; } gprev1 = g1; gprev2 = g2; double sold = Mass2(g1,g2); if (Q2(g1) > Q2(g2)) { if (g1==gbeam1) continue; //printf("M1 %d\n", g1->M1); if (g1->M1<0) continue; g1 = (TRootGenParticle*)branchGen->At(g1->M1); x[0] *= sold/Mass2(g1,g2); //printf("x1_step %.3f\n", sold/Mass2(g1,g2)); } else { if (g2==gbeam2) continue; if (g2->M1<0) continue; //printf("M2 %d\n", g2->M1); g2 = (TRootGenParticle*)branchGen->At(g2->M1); x[1] *= sold/Mass2(g1,g2); //printf("x2_step %.3f\n", sold/Mass2(g1,g2)); } } return &x; } //>>>>>>>>>> END FILE >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> //############################################################################ CASE 2: YOU KNOW ALREADY X1, X2 and Q**2 (or good approximations to them) AND WANT TO CALCULATE THE REWEIGHTING FACTORS IN ORDER TO GO FROM THE CENTRAL SET OF A PDF TO THE SET NUMBER "I" (OF THE SAME PDF //############################################################################ a) Copy the latest LHAPDF package version into your machine. For instance: > cd MY_LHAPDF_DIR > wget http://www.cedar.ac.uk/downloads/lhapdf-4.2.tar.gz > tar xvfvz lhapdf-4.2.tar.gz > cd lhapdf-4.2 > ./configure; make b) Get my private "ccwrap" version ( a modified version from the original one by Stefan Gieseke at http://hepforge.cedar.ac.uk/lhapdf/cc): > wget http://cern.ch/alcaraz/ccwrap.tgz > tar xvfz ccwrap.tgz > cd ccwrap > make (make sure before that your ROOT stuff is correctly setup) c) Make a link to the LHAPDF stuff in yous analysis directory: > cd MY_ANALYSIS_DIR > ln -s MY_LHAPDF_DIR/lhapdf-4.2 LHAPDF /// Here it is an (almost) imaginary Root example showing how the reweighting step should be implemented: //>>>>>>>>>> BEGIN FILE >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> // ROOTSYS=/afs/cern.ch/cms/external/lcg/external/root/3.10.02/slc3_ia32_gcc323/root // // The example should run as: // $ root.exe rewtest.C // If you want to compile it first, then proceed as follows. // First uncomment the include line for "LHAPDFWRap.h" and then: // $ root.exe // > gSystem->Load("LHAPDF/src/libLHAPDFWrap.so"); // > .x rewtest.C+ #include "TSystem.h" #include "TH1F.h" #include "TRandom.h" //#include "LHAPDF/ccwrap/LHAPDFWrap.h" int rewtest(){ gSystem->Load("LHAPDF/src/libLHAPDFWrap.so"); char pdfset[80]; sprintf(pdfset, "%s", "cteq61"); int subseta = 0; int subsetb = 30; char pdfname[80]; sprintf(pdfname, "LHAPDF/PDFsets/%s.LHgrid", pdfset); FILE* ftest = fopen(pdfname, "r"); if (!ftest) { sprintf(pdfname, "LHAPDF/PDFsets/%s.LHpdf", pdfset); ftest = fopen(pdfname, "r"); } if (!ftest) { printf("NO files found for PDF: %s\n", pdfset); return -1; } fclose(ftest); const int MAXEVENTS = 1000; const double XMIN = 0.0001; const double XMAX = 0.5; TH1F* hbefore = new TH1F ("hbefore", "x histogram before", 50, XMIN, XMAX); TH1F* hafter = new TH1F ("hafter", "x histogram after", 50, XMIN, XMAX); // Your loop on events: LHAPDFWrap pdf(pdfname, 0); pdf.setSilent(); for (int event=0; eventRndm(); double x2 = XMIN + (XMAX-XMIN)*gRandom->Rndm(); // Find the parton flavour (we will assume uubar here) // Gluon is id=0; quark flavours follow the PDG numbering scheme int ifl1 = 2; int ifl2 = -2; pdf.initPDF(subseta); double w1a = pdf.xfx(x1, Qval, ifl1); double w2a = pdf.xfx(x2, Qval, ifl2); //printf("w1a %f w2a %f\n", w1a, w2a); pdf.initPDF(subsetb); double w1b = pdf.xfx(x1, Qval, ifl1); double w2b = pdf.xfx(x2, Qval, ifl2); //printf("w1b %f w2b %f\n", w1b, w2b); double weight = (w1b*w2b)/(w1a*w2a); //printf("x1 %f x2 %f, event weight: %f\n", x1, x2, weight); hbefore->Fill(x1); hbefore->Fill(x2); hafter->Fill(x1, weight); hafter->Fill(x2, weight); } hbefore->SetLineStyle(1); hafter->SetLineStyle(2); hafter->Draw(); hbefore->Draw("same"); return 1; } //>>>>>>>>>> END FILE >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> //>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> //############################################################################ CASE 3: YOU KNOW ALREADY X1, X2 and Q**2 (or good aproximations to them) AND WANT TO CALCULATE THE REWEIGHTING FACTORS IN ORDER TO GO FROM THE CENTRAL SET OF A PDF TO THE CENTRAL SET OF ANOTHER ONE. //############################################################################ You can of course do it but, unfortunately, you need to reinitialize the two sets of PDFS FOR EACH EVENT. This usually takes an enormous amount of time. Here below you can find an example. In general, I advise not to do this, but to proceed in three steps: a) Store the x*pdf values in a vector in memory in binning form. I have tabulated them in the past using 5000 bins for x between 0 and 1, for all 13 flavours (6 quarks, 6 antiquarks and the gluon) and for a given Q2 (which is fixed by the nature of the process under study). This means a vector of 5000*13 components. b) Calculate the 5000*13 ratios of the x*pdf values of the final set over the x*pdf values of the original PDF set. Then write them into a file (or into 13 histograms of 5000 entries, if you prefer). c) Read those "w(i,flavour)" 5000*13 weights back and reweight event by event by the product of w(i1,flavour1)*w(i2,flavour2), being i1, i2 the bins corresponding to the observed x1 and x2 values; and being flavour1 and flavour2 the corresponding parton flavours. //>>>>>>>>>> BEGIN FILE >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> // WARNING!!! THIS EXAMPLE IS EXTREMELY TIME EXPENSIVE. // // ROOTSYS=/afs/cern.ch/cms/external/lcg/external/root/3.10.02/slc3_ia32_gcc323/root // // The example should run as: // $ root.exe rewtest.C // IF you want to compile it first, then proceed as follows. // First uncomment the include line for "LHAPDFWRap.h" and then: // $ root.exe // > gSystem->Load("LHAPDF/src/libLHAPDFWrap.so"); // > .x rewtest2.C+ #include "TSystem.h" #include "TH1F.h" #include "TRandom.h" //#include "LHAPDF/ccwrap/LHAPDFWrap.h" int rewtest2(){ gSystem->Load("LHAPDF/src/libLHAPDFWrap.so"); char pdfseta[80]; sprintf(pdfseta, "%s", "cteq61"); char pdfsetb[80]; sprintf(pdfsetb, "%s", "MRST2001E"); int subseta = 0; int subsetb = 0; char pdfnamea[80]; sprintf(pdfnamea, "LHAPDF/PDFsets/%s.LHgrid", pdfseta); FILE* ftest = fopen(pdfnamea, "r"); if (!ftest) { sprintf(pdfnamea, "LHAPDF/PDFsets/%s.LHpdf", pdfseta); ftest = fopen(pdfnamea, "r"); } if (!ftest) { printf("NO files found for PDF: %s\n", pdfseta); return -1; } fclose(ftest); char pdfnameb[80]; sprintf(pdfnameb, "LHAPDF/PDFsets/%s.LHgrid", pdfsetb); ftest = fopen(pdfnameb, "r"); if (!ftest) { sprintf(pdfnameb, "LHAPDF/PDFsets/%s.LHpdf", pdfsetb); ftest = fopen(pdfnameb, "r"); } if (!ftest) { printf("NO files found for PDF: %s\n", pdfsetb); return -1; } fclose(ftest); const int MAXEVENTS = 10; const double XMIN = 0.0001; const double XMAX = 0.5; TH1F* hbefore = new TH1F ("hbefore", "x histogram before", 50, XMIN, XMAX); TH1F* hafter = new TH1F ("hafter", "x histogram after", 50, XMIN, XMAX); // Your loop on events: LHAPDFWrap pdf(pdfnamea, 0); pdf.setSilent(); for (int event=0; eventRndm(); double x2 = XMIN + (XMAX-XMIN)*gRandom->Rndm(); // Find the parton flavour (we will assume uubar here) // Gluon is id=0; quark flavours follow the PDG numbering scheme int ifl1 = 2; int ifl2 = -2; pdf = LHAPDFWrap(pdfnamea, subseta); pdf.initPDF(subseta); double w1a = pdf.xfx(x1, Qval, ifl1); double w2a = pdf.xfx(x2, Qval, ifl2); //printf("w1a %f w2a %f\n", w1a, w2a); pdf = LHAPDFWrap(pdfnameb, subsetb); pdf.initPDF(subsetb); double w1b = pdf.xfx(x1, Qval, ifl1); double w2b = pdf.xfx(x2, Qval, ifl2); //printf("w1b %f w2b %f\n", w1b, w2b); double weight = (w1b*w2b)/(w1a*w2a); //printf("x1 %f x2 %f, event weight: %f\n", x1, x2, weight); hbefore->Fill(x1); hbefore->Fill(x2); hafter->Fill(x1, weight); hafter->Fill(x2, weight); } hbefore->SetLineStyle(1); hafter->SetLineStyle(2); hafter->Draw(); hbefore->Draw("same"); return 1; } //>>>>>>>>>> END FILE >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>