/************************************************************************ * WGRAD_HistModule * * ---------------- * * * * D. Waters * * 21.8.99 * *************************************************************************/ //----------------------- // This Class's Header -- //----------------------- #include "Electroweak/WGRAD_HistModule.hh" //------------- // C Headers -- //------------- #include //------------------------------- // Collaborating Class Headers -- //--------------------------- #include "BaBar/Cdf.hh" #ifdef USE_CDFEDM2 #include "Edm/ConstHandle.hh" #include "Edm/EventRecord.hh" #include "SimulationObjects/HEPG_StorableBank.hh" #else #include "Banks/HEPG_Bank.hh" #include "Trybos/TRY_Abstract_Record.hh" #endif // USE_CDFEDM2 #include "CLHEP/Matrix/Vector.h" #include "CLHEP/Vector/LorentzVector.h" #include "CLHEP/Random/Ranlux64Engine.h" #include "CLHEP/Random/RandGamma.h" #include "CLHEP/Random/RandGauss.h" #include "Electroweak/HighPtElectrons.hh" #include "HepTuple/HepHist1D.h" // ------------------------------------------------- // -- Declarations of variables defined elsewhere -- // ------------------------------------------------- //----------------------------------------------------------------------- // Local Macros, Typedefs, Structures, Unions and Forward Declarations -- //----------------------------------------------------------------------- static const char rcsid[] = " "; //---------------- // Constructors -- //---------------- WGRAD_HistModule::WGRAD_HistModule(const char* const theName, const char* const theDescription ) : HepHistModule( theName, theDescription ), _randomEngine(0),_gammaGenerator(0),_gaussGenerator(0) { // Set up the CLHEP random number generating machinery : _randomEngine = new Ranlux64Engine(280273); // Initialize the generators with the random number engine. Send as // a reference to tell the generators to leave memory cleanup to us. //----------------------------------------------------------------------// // NOTE : _gammaGenerator->shoot(X,1.0) produces the same distribution // // as CERNLIB routine RNGAMA(X). // // ---------------------------------------------------------------------// if (_randomEngine) { _gammaGenerator = new RandGamma(*_randomEngine); _gaussGenerator = new RandGauss(*_randomEngine); } } //-------------- // Destructor -- //-------------- WGRAD_HistModule::~WGRAD_HistModule( ) { if (_gammaGenerator) { delete _gammaGenerator; _gammaGenerator = NULL; } if (_gaussGenerator) { delete _gaussGenerator; _gaussGenerator = NULL; } if (_randomEngine) { delete _randomEngine; _randomEngine = NULL; } } //-------------- // Operations -- //-------------- AppResult WGRAD_HistModule::bookHistograms( ) { HepFileManager* manager = fileManager( ); assert(0 != manager); _Z_pt_fit = &manager->hist1D("Z pT fit (GeV)", 500 , 0.0, 50.0, 100); assert(0 != _Z_pt_fit); _W_pt = &manager->hist1D("W pT (GeV)", 320 , 0.0, 160.0, 200); assert(0 != _W_pt); _W_recoil = &manager->hist1D("W recoil (GeV)", 100 , 0.0, 40.0, 201); assert(0 != _W_recoil); _W_recoilReWeighted = &manager->hist1D("W recoil RW (GeV)", 100 , 0.0, 40.0, 211); assert(0 != _W_recoilReWeighted); _W_tmass = &manager->hist1D("W Mt (GeV)", 200 , 0.0, 200.0, 202); assert(0 != _W_tmass); _W_pt_resum = &manager->hist1D("W pT resum (GeV)", 500 , 0.0, 50.0, 300); assert(0 != _W_pt_resum); // ------------------------------------------------------------------------------- // // Dave Waters 24.8.99 // // Fill a histogram with a fit to the measured Z pT distribution provided by // // Mark Lancaster. // // --------------------------------------------------------------------------------// // Set the number of bins and upper limit on pT : int nbins = _Z_pt_fit->nBins(); double ptmax = _Z_pt_fit->max(); double ptmin = _Z_pt_fit->min(); // Make the original Z pT distribution from the parameterisation up to 50 GeV : double binwidth = (ptmax-ptmin)/nbins; for (int i = 0; i < nbins; i++) { double pt = ptmin+(i+0.5)*binwidth; double x = pt/50.0; _Z_pt_orig[i] = zpt(x); _Z_pt_fit->accumulate(pt,zpt(x)); } // Fill weighting arrays : _pythia_weights_ktdefault[0] = 0.448492; // 0 -> 0.5 GeV _pythia_weights_ktdefault[1] = 0.333937; // 0.5 -> 1.0 GeV _pythia_weights_ktdefault[2] = 0.304173; // ... etc _pythia_weights_ktdefault[3] = 0.294817; _pythia_weights_ktdefault[4] = 0.305336; _pythia_weights_ktdefault[5] = 0.335484; _pythia_weights_ktdefault[6] = 0.354203; _pythia_weights_ktdefault[7] = 0.385636; _pythia_weights_ktdefault[8] = 0.411423; _pythia_weights_ktdefault[9] = 0.457582; _pythia_weights_ktdefault[10] = 0.453598; _pythia_weights_ktdefault[11] = 0.488375; _pythia_weights_ktdefault[12] = 0.504292; _pythia_weights_ktdefault[13] = 0.537012; _pythia_weights_ktdefault[14] = 0.517696; _pythia_weights_ktdefault[15] = 0.523655; _pythia_weights_ktdefault[16] = 0.555729; _pythia_weights_ktdefault[17] = 0.564055; _pythia_weights_ktdefault[18] = 0.554538; _pythia_weights_ktdefault[19] = 0.582682; _pythia_weights_ktdefault[20] = 0.595276; _pythia_weights_ktdefault[21] = 0.591751; _pythia_weights_ktdefault[22] = 0.593324; _pythia_weights_ktdefault[23] = 0.610165; _pythia_weights_ktdefault[24] = 0.566659; _pythia_weights_ktdefault[25] = 0.595331; _pythia_weights_ktdefault[26] = 0.597675; _pythia_weights_ktdefault[27] = 0.551751; _pythia_weights_ktdefault[28] = 0.574142; _pythia_weights_ktdefault[29] = 0.576573; _pythia_weights_ktdefault[30] = 0.55632; _pythia_weights_ktdefault[31] = 0.553323; _pythia_weights_ktdefault[32] = 0.585789; _pythia_weights_ktdefault[33] = 0.567996; _pythia_weights_ktdefault[34] = 0.570269; _pythia_weights_ktdefault[35] = 0.578268; _pythia_weights_ktdefault[36] = 0.565487; _pythia_weights_ktdefault[37] = 0.553053; _pythia_weights_ktdefault[38] = 0.480225; _pythia_weights_ktdefault[39] = 0.522937; return AppResult::OK; } AppResult WGRAD_HistModule::fillHistograms( AbsEvent* anEvent ) { // std::cout << " WGRAD_HistModule::fillHistograms : " << std::endl; double lepton[4],neutrino[4]; bool lepton_found = false; bool neutrino_found = false; // Try negotiating the HEPG bank : #ifdef USE_CDFEDM2 for (EventRecord::ConstIterator hepgIter(anEvent, "HEPG_StorableBank"); hepgIter.is_valid(); ++hepgIter) { ConstHandle hepgCH(hepgIter); const HEPG_StorableBank& hepg = *hepgCH; #else for (TRY_Record_Iter_Same hepgIter(anEvent, "HEPG"); hepgIter.is_valid(); ++hepgIter ) { HEPG_Bank hepg(hepgIter); #endif // USE_CDFEDM2 // Iterate over the particles in the HEPG bank : #ifdef USE_CDFEDM2 for (HEPG_StorableBank::particleIter partIter(hepg); partIter.is_valid(); ++partIter) #else for (HEPG_Bank::particleIter partIter(hepg); partIter.is_valid(); ++partIter) #endif // USE_CDFEDM2 { // Look for the first charged lepton and neutrino in the event with ISTHEP = 1 : if (hepg.istdhep(partIter) == 1) { int absPartCode = abs(hepg.idhep(partIter)); if (!lepton_found && (absPartCode == 11 || absPartCode == 13 || absPartCode == 15)) { // std::cout << " WGRAD_HistModule::fillHistograms : charged lepton found " << std::endl; lepton_found = true; lepton[0] = hepg.Px(partIter); lepton[1] = hepg.Py(partIter); lepton[2] = hepg.Pz(partIter); lepton[3] = hepg.E(partIter); } if (!neutrino_found && (absPartCode == 12 || absPartCode == 14 || absPartCode == 16)) { // std::cout << " WGRAD_HistModule::fillHistograms : neutrino found " << std::endl; neutrino_found = true; neutrino[0] = hepg.Px(partIter); neutrino[1] = hepg.Py(partIter); neutrino[2] = hepg.Pz(partIter); neutrino[3] = hepg.E(partIter); } } } // Make a Lorentz vector by combining the lepton and neutrino 4-momenta : HepLorentzVector W((lepton[0]+neutrino[0]), (lepton[1]+neutrino[1]), (lepton[2]+neutrino[2]), (lepton[3]+neutrino[3])); double W_pt = sqrt(pow(W.px(),2)+pow(W.py(),2)); HepVector transLepton(2),transNeutrino(2); transLepton[0] = lepton[0]; transLepton[1] = lepton[1]; transNeutrino[0] = neutrino[0]; transNeutrino[1] = neutrino[1]; double W_tmass = HighPtElectrons::transverseMass(transLepton,transNeutrino); // Fill some histograms with generator level information : _W_pt->accumulate(W_pt); _W_tmass->accumulate(W_tmass); // Does this event, after applying the recoil model, pass the cuts ? bool acceptEvent = passCuts(W.px(),W.py(),lepton[0],lepton[1]); if (acceptEvent) { _W_recoil->accumulate(_recoilModelRecoil); // Reweight on the basis of the true W pT. This will take Pythia (default kT) -> // correct (i.e. measured from Z pT) W pT distribution. int wptBin = W_pt/0.5 + 1; wptBin = min(40,wptBin); // std::cout << " W_pt, bin = " << W_pt << " , " << wptBin << std::endl; _W_recoilReWeighted->accumulate(_recoilModelRecoil,_pythia_weights_ktdefault[wptBin-1]); // Fill a histogram with the W pT distribution for this event on the basis of // the W rapidity : int nbins = _W_pt_resum->nBins(); double ptmax = _W_pt_resum->max(); double ptmin = _W_pt_resum->min(); double binwidth = (ptmax-ptmin)/nbins; double ptArray[1000],scaledPtArray[1000]; double scaledPtArraySum = 0.; for (int i = 0; i < nbins; i++) { double pt = ptmin+(i+0.5)*binwidth; ptArray[i] = pt; scaledPtArray[i] = _Z_pt_orig[i]* zypt(pt,fabs(W.rapidity()),1)* zwpt(pt,fabs(W.rapidity()),1); scaledPtArraySum += scaledPtArray[i]; // std::cout << " pt, zwpt = " << pt << " , " << zwpt(pt,fabs(W.rapidity()),1) << std::endl; // std::cout << " W rapidity; zypt; zwpt = " << W.rapidity() << " , " // << zypt(pt,W.rapidity(),1) << " , " << zwpt(pt,W.rapidity(),1) << std::endl; } // Fill histogram with normalised result : for (int i = 0; i < nbins; i++) { _W_pt_resum->accumulate(ptArray[i],(scaledPtArray[i]/scaledPtArraySum)); } } } return AppResult::OK; } AppModule* WGRAD_HistModule::clone(const char* cloneName) { return new WGRAD_HistModule(cloneName,"this module is a clone WGRAD_HistModule"); } const char * WGRAD_HistModule::rcsId( ) const { return rcsid; } double WGRAD_HistModule::zpt(double x) { double par[4]; par[0] = 0.51356; par[1] = 0.27675; par[2] = 7.8193; par[3] = 0.41704; double g1,g2,g3,a1,a2; double gg = lgamma(par[3]); // log(Gamma(x)):log(n!)=log( Gamma(n+1) ) g1 = par[3]*exp(gg) ; g2 = pow(x,par[3]) ; g3 = g2/g1 ; a1 = (1.0-par[0])*pow(par[1],(par[3]+1))*exp(-par[1]*x) ; a2 = (0.0+par[0])*pow(par[2],(par[3]+1))*exp(-par[2]*x) ; return g3*(a1+a2); } float WGRAD_HistModule::zypt(float pt, float y, int iflag){ /* --- Given an input pT and Y this function returns a multiplicative --- weight to apply to the event to convert the pT spectrum from one --- for Z events generated at one Y (assumed to be the mean after cuts) --- i.e Y=0.3 to one for any Y. --- The pT spectrum is softer at high-Y. This basically does this. --- i.e. turns an average pT spectrum into one differential in Y. --- The ratio is taken from running Arnold-Kaufmann, to alphs^2**2 --- in b-space with LY parameters except g2=0.6. --- Supports two arguements ----- IFLAG=1 : MRS-R2 ----- IFLAG=2 : MRS-R1 --- The function is based on fitting the ratio in 10 Y bins as a function --- of pT. The coefficients of the pT ratio fit are then fit as a --- function of Y. Mark Lancaster : 31st July 1998. ------------------------------------------------------------------- */ static float p0, p1, p2, p3, p4, y1, y2, y3, pt1, pt2, pt3, pt4, p0_y[4], p1_y[4], p2_y[4], p3_y[4], p4_y[4]; static int first = 0 ; /* -->> Files "z_y_ratio_r*.dat" from running ~/wmass/pt/arnold2/make_y_z_ratio.f */ pt1 = pt; pt2 = pt1 * pt; pt3 = pt2 * pt; pt4 = pt3 * pt; y1 = y; y2 = y1 * y; y3 = y2 * y; /* : coefficients of the Y dependence of the ratio fit... */ if (first == 0 && iflag == 1) { /* MRS-R2 */ printf("\n-- Z P(T) Y REWEIGHTING FROM MRS-R2 --- \n") ; p0_y[0] = .993285; p0_y[1] = .0260119; p0_y[2] = -.00789399; p0_y[3] = .025722; p1_y[0] = .0014923; p1_y[1] = -.00688932; p1_y[2] = .00616331; p1_y[3] = -.00664817; p2_y[0] = -7.14482e-5; p2_y[1] = 3.78316e-4; p2_y[2] = -5.48541e-4; p2_y[3] = 3.38331e-4; p3_y[0] = 1.43036e-6; p3_y[1] = -7.92301e-6; p3_y[2] = 1.32963e-5; p3_y[3] = -7.36784e-6; p4_y[0] = -1.04138e-8; p4_y[1] = 6.25338e-8; p4_y[2] = -1.181e-7; p4_y[3] = 6.01555e-8; } if (first == 0 && iflag == 2) { /* MRS-R1 */ printf("\n-- Z P(T) Y REWEIGHTING FROM MRS-R1 --- \n") ; p0_y[0] = .993844; p0_y[1] = .0243082; p0_y[2] = -.00871364; p0_y[3] = .0238222; p1_y[0] = .00143883; p1_y[1] = -.00684189; p1_y[2] = .00639546; p1_y[3] = -.00653028; p2_y[0] = -7.31236e-5; p2_y[1] = 4.03191e-4; p2_y[2] = -5.89889e-4; p2_y[3] = 3.47032e-4; p3_y[0] = 1.63919e-6; p3_y[1] = -9.49411e-6; p3_y[2] = 1.54553e-5; p3_y[3] = -8.04981e-6; p4_y[0] = -1.38028e-8; p4_y[1] = 8.54909e-8; p4_y[2] = -1.48003e-7; p4_y[3] = 7.03315e-8; } first = 1; p0 = p0_y[0] + p0_y[1] * y1 + p0_y[2] * y2 + p0_y[3] * y3; p1 = p1_y[0] + p1_y[1] * y1 + p1_y[2] * y2 + p1_y[3] * y3; p2 = p2_y[0] + p2_y[1] * y1 + p2_y[2] * y2 + p2_y[3] * y3; p3 = p3_y[0] + p3_y[1] * y1 + p3_y[2] * y2 + p3_y[3] * y3; p4 = p4_y[0] + p4_y[1] * y1 + p4_y[2] * y2 + p4_y[3] * y3; return p0 + p1 * pt1 + p2 * pt2 + p3 * pt3 + p4 * pt4; } float WGRAD_HistModule::zwpt(float pt, float yy, int iflag){ /* --- Given an input pT and Y this function returns a multiplicative --- weight to apply to the event to convert the pT spectrum from one --- for Z events to one for W events. --- The ratio is taken from running Arnold-Kaufmann, to alphas^2**2. --- Supports IFLAG = 1->4 ----- IFLAG=1 : MRS-R2 : b-space : g2 = 0.60 ----- IFLAG=2 : MRS-R2 : b-space : g2 = 0.45 (low one sigma) ----- IFLAG=3 : MRS-R2 : b-space : g2 = 0.75 (high one sigma) ----- IFLAG=4 : MRS-R1 : b-space : g2 = 0.60 (new pdf) --- The function is based on fitting the ratio in 10 Y bins as a function --- of pT. The coefficients of the pT ratio fit are then fit as a --- function of Y. Mark Lancaster : 31st July 1998. ------------------------------------------------------------------- -->> Files "r2_b_g2060.y_p_fits.f" etc from running ~/wmass/pt/arnold2/make_zwpt_ratio.f --- Fits as a function of Y only really valid in reagion y < 2.0 --- So if y > 2.0 set y = 2.0. */ static int first = 0 ; float ret_val,y ; static float p0, p1, p2, p3, p4, y1, y2, y3, pt1, pt2, p0_y[4], p1_y[4], p2_y[4], p3_y[4], p4_y[4]; if (yy > 2.0) y = 2.0 ; if (yy <= 2.0) y = yy ; pt1 = pt; pt2 = pt1 * pt; y1 = y; y2 = y1 * y; y3 = y2 * y; /* : coefficients of the Y dependence of the ratio fit... */ if (first == 0 && iflag == 1) { printf("\n -- Z/W pT REWEIGHT : MRS-R2 : b-space:g2 =0.6 \n"); p0_y[0] = 1.00598; p0_y[1] = -.107536; p0_y[2] = .0494103; p0_y[3] = -.00306264; p1_y[0] = .00480112; p1_y[1] = 1.56684; p1_y[2] = -.17463; p1_y[3] = -.185392; p2_y[0] = 1.08705; p2_y[1] = .853719; p2_y[2] = -.481549; p2_y[3] = .172936; p3_y[0] = -.0808729; p3_y[1] = -.0139689; p3_y[2] = -.0237642; p3_y[3] = 1.07492e-4; p4_y[0] = -.00194144; p4_y[1] = -.00236383; p4_y[2] = .00398231; p4_y[3] = -.0011232; } if (first == 0 && iflag == 2) { printf("\n -- Z/W pT REWEIGHT : MRS-R2 : b-space:g2 =0.45\n"); p0_y[0] = 1.01169; p0_y[1] = -.104554; p0_y[2] = .0512159; p0_y[3] = -.00447775; p1_y[0] = -.0135311; p1_y[1] = 1.48419; p1_y[2] = -.326563; p1_y[3] = -.110707; p2_y[0] = 1.18464; p2_y[1] = .662812; p2_y[2] = -.270765; p2_y[3] = .113823; p3_y[0] = -.0209091; p3_y[1] = -.169685; p3_y[2] = .0705705; p3_y[3] = -.00935475; p4_y[0] = -.00206341; p4_y[1] = -.0021843; p4_y[2] = .00367509; p4_y[3] = -9.99954e-4; } if (first == 0 && iflag == 3) { printf("\n -- Z/W pT REWEIGHT : MRS-R2 : b-space:g2 =0.75\n"); p0_y[0] = 1.00435; p0_y[1] = -.117018; p0_y[2] = .0545736; p0_y[3] = -.00379051; p1_y[0] = -.0119275; p1_y[1] = 1.82437; p1_y[2] = -.283149; p1_y[3] = -.177396; p2_y[0] = .95542; p2_y[1] = 1.13133; p2_y[2] = -.6962; p2_y[3] = .209955; p3_y[0] = -.0605829; p3_y[1] = -.07905; p3_y[2] = .0347798; p3_y[3] = -.0125066; p4_y[0] = -.00190873; p4_y[1] = -.00220052; p4_y[2] = .00379941; p4_y[3] = -.00106369; } if (first == 0 && iflag == 4) { printf("\n -- Z/W pT REWEIGHT : MRS-R1 : b-space:g2 =0.6\n"); p0_y[0] = 1.0023; p0_y[1] = -.106408; p0_y[2] = .0515487; p0_y[3] = -.00363927; p1_y[0] = -.0472909; p1_y[1] = 2.11157; p1_y[2] = -.791024; p1_y[3] = -.0288853; p2_y[0] = 1.15169; p2_y[1] = 1.14896; p2_y[2] = -.759639; p2_y[3] = .253782; p3_y[0] = -.105211; p3_y[1] = -.00112298; p3_y[2] = -.0297693; p3_y[3] = -.00350382; p4_y[0] = -.00173566; p4_y[1] = -.0020258; p4_y[2] = .00369263; p4_y[3] = -.00106164; } first = 1; p0 = p0_y[0] + p0_y[1] * y1 + p0_y[2] * y2 + p0_y[3] * y3; p1 = p1_y[0] + p1_y[1] * y1 + p1_y[2] * y2 + p1_y[3] * y3; p2 = p2_y[0] + p2_y[1] * y1 + p2_y[2] * y2 + p2_y[3] * y3; p3 = p3_y[0] + p3_y[1] * y1 + p3_y[2] * y2 + p3_y[3] * y3; p4 = p4_y[0] + p4_y[1] * y1 + p4_y[2] * y2 + p4_y[3] * y3; if (pt1 <= 6.75) { ret_val = p0 + p1 * exp(-p2 * pt1 - p3 * pt2); } else if (pt1 > 6.75) { ret_val = p0 + p1 * exp(-p2 * 6.75 - p3 * 45.5625) + p4 * (pt1 - 6.75); } return ret_val; } bool WGRAD_HistModule::passCuts(double Wpx, double Wpy, double Mpx, double Mpy) { // Based on some run I fortran from Mark Lancaster. // Simplified - e.g. no PHOTON vector. No thought about // how backgrounds effect the recoil model parameters. // --- SUMET PARAMETERS:From Gamma fit to Z data vs PT -------------------- double SUMET_P1[2],SUMET_P2[2]; SUMET_P1[0] = 1.2852; SUMET_P1[1] = 0.099309; SUMET_P2[0] = 22.421; SUMET_P2[1] = -0.21066; // --- Min Bias resolution vs Sumet ---------------------------------- double MINB_P1,MINB_P2; MINB_P1 = 0.324; MINB_P2 = 0.577; // --- 7 Recoil model parameters (see CDF 4770) // --- Values in 4770 are from fits to Z data. // --- These values are the W-optimised values. double A,B,C,D,E,F,G; // -- U1 response function A = -.50798; B = -.57666; C = -.00051; // -- U1 resolution parameters D = 0.94174; E = 0.01860; // -- U2 resolution parameters F = 0.90071; G = -.00556; double PT = sqrt(pow(Wpx,2)+pow(Wpy,2)); // // --- Form a SUMET for the event (based on fits to Z data) // double SP1 = SUMET_P1[0] + SUMET_P1[1]*PT; double SP2 = SUMET_P2[0] + SUMET_P2[1]*PT; double SUMET = (_gammaGenerator->shoot(_randomEngine, SP1+1.0, 1.0))*SP2; // // --- Generate the min-bias resolution (based on fits to MB data) // double SIGMB = MINB_P1*pow(SUMET,MINB_P2); // --- Response function double MU1 = min(0., A + B*PT + C*pow(PT,2) ); // // --- Resolution functions // double SIG1 = SIGMB * (D + E*PT); double SIG2 = SIGMB * (F + G*PT); double G1 = _gaussGenerator->shoot(_randomEngine); double G2 = _gaussGenerator->shoot(_randomEngine); // // --- Form U1, u2 -- projections of U along VB direction. // double U1 = MU1 + G1*SIG1; double U2 = 0. + G2*SIG2; // // --- Rotate to X,Y -- use true phi for rotation (slightly dodgy) // double PHI = atan2(Wpy,Wpx); double pi2val = 6.2831853; if (PHI < 0.) PHI = PHI + pi2val; double PWSEEN[2]; PWSEEN[0] = -1.0*(U1*cos(PHI) - U2*sin(PHI)); PWSEEN[1] = -1.0*(U1*sin(PHI) + U2*cos(PHI)); // // --- Need to add any photons into U here + get the sign of PHOTONS() correct. // // WARNING : PHOTONS NOT TAKEN INTO ACCOUNT : // ========================================== double PHOTONS[2]; PHOTONS[0] = 0.; PHOTONS[1] = 0.; PWSEEN[0] = PWSEEN[0] - PHOTONS[0]; PWSEEN[1] = PWSEEN[1] - PHOTONS[1]; // // --->> This is what DO are calling W pT !! // double U = sqrt(pow(PWSEEN[0],2) + pow(PWSEEN[1],2)); double ETMISS[2]; // WARNING : MUONS NOT TAKEN INTO ACCOUNT : // ======================================== double PMUON[2];; PMUON[0] = Mpx; PMUON[1] = Mpy; ETMISS[0] = PWSEEN[0] - PMUON[0]; ETMISS[1] = PWSEEN[1] - PMUON[1]; double PTMUON = sqrt(pow(PMUON[0],2)+pow(PMUON[1],2)); double PTMISS = sqrt(pow(ETMISS[0],2)+pow(ETMISS[1],2)); double MT = sqrt(2.*(PTMUON*PTMISS-PMUON[0]*ETMISS[0]-PMUON[1]*ETMISS[1])); // Set the event properties after applying the recoil model : _recoilModelRecoil = U; _recoilModelT_Mass = MT; _recoilModelPtMiss = PTMISS; // // --- Make cuts on U < 20.0, MET > 25.0, ET > 25.0, 65 < MT < 100 // if (U < 20.0 && PTMISS > 25.0 && PTMUON > 25.0 && MT > 65.0 && MT < 100.0) { return true; } else { return false; } }