#include <EmissionVetoHook1.h>

Inheritance diagram for EmissionVetoHook1:

Public Member Functions
bool	canVetoFSREmission ()

bool	canVetoISREmission ()

bool	canVetoMPIEmission ()

bool	canVetoMPIStep ()

bool	doVetoFSREmission (int, const Pythia8::Event &e, int iSys, bool)

bool	doVetoISREmission (int, const Pythia8::Event &e, int iSys)

bool	doVetoMPIEmission (int, const Pythia8::Event &e)

bool	doVetoMPIStep (int nMPI, const Pythia8::Event &e)

	EmissionVetoHook1 (int nFinalIn, bool vetoOnIn, int vetoCountIn, int pThardModeIn, int pTemtModeIn, int emittedModeIn, int pTdefModeIn, bool MPIvetoOnIn, int VerbosityIn)

void	fatalEmissionVeto (string message)

int	numberVetoMPIStep ()

double	pTcalc (const Pythia8::Event &e, int i, int j, int k, int r, int xSRin)

double	pTpowheg (const Pythia8::Event &e, int i, int j, bool FSR)

double	pTpythia (const Pythia8::Event &e, int RadAfterBranch, int EmtAfterBranch, int RecAfterBranch, bool FSR)

	~EmissionVetoHook1 ()

Private Attributes
bool	accepted

int	emittedMode

int	MPIvetoOn

int	nAcceptSeq

int	nFinal

int	nFinalExt

unsigned long int	nFSRveto

unsigned long int	nISRveto

int	pTdefMode

int	pTemtMode

double	pThard

int	pThardMode

double	pTMPI

int	Verbosity

int	vetoCount

int	vetoOn

Detailed Description

Definition at line 3 of file EmissionVetoHook1.h.

Constructor & Destructor Documentation

EmissionVetoHook1::EmissionVetoHook1	(	int	nFinalIn,
		bool	vetoOnIn,
		int	vetoCountIn,
		int	pThardModeIn,
		int	pTemtModeIn,
		int	emittedModeIn,
		int	pTdefModeIn,
		bool	MPIvetoOnIn,
		int	VerbosityIn
	)

inline

Definition at line 8 of file EmissionVetoHook1.h.

                                                                         :
                     nFinalExt(nFinalIn),
                     vetoOn(vetoOnIn), vetoCount(vetoCountIn),
                     pThardMode(pThardModeIn), pTemtMode(pTemtModeIn),
                     emittedMode(emittedModeIn), pTdefMode(pTdefModeIn),
                     MPIvetoOn(MPIvetoOnIn), nISRveto(0), nFSRveto(0),
                     Verbosity(VerbosityIn) {}

EmissionVetoHook1::~EmissionVetoHook1 ( )

inline

Definition at line 17 of file EmissionVetoHook1.h.

References gather_cfg::cout, nFSRveto, and nISRveto.

                       {
     cout << "Number of ISR vetoed = " << nISRveto << endl;
     cout << "Number of FSR vetoed = " << nFSRveto << endl;
   }

Member Function Documentation

bool EmissionVetoHook1::canVetoFSREmission ( )

inline

Definition at line 31 of file EmissionVetoHook1.h.

References vetoOn.

31 { return vetoOn; }

EmissionVetoHook1::vetoOn

int vetoOn

Definition: EmissionVetoHook1.h:49

bool EmissionVetoHook1::canVetoISREmission ( )

inline

Definition at line 28 of file EmissionVetoHook1.h.

References vetoOn.

28 { return vetoOn; }

EmissionVetoHook1::vetoOn

int vetoOn

Definition: EmissionVetoHook1.h:49

bool EmissionVetoHook1::canVetoMPIEmission ( )

inline

Definition at line 34 of file EmissionVetoHook1.h.

References MPIvetoOn.

34 { return MPIvetoOn; }

EmissionVetoHook1::MPIvetoOn

int MPIvetoOn

Definition: EmissionVetoHook1.h:49

bool EmissionVetoHook1::canVetoMPIStep ( )

inline

Definition at line 24 of file EmissionVetoHook1.h.

24 { return true; }

bool EmissionVetoHook1::doVetoFSREmission	(	int	,
		const Pythia8::Event &	e,
		int	iSys,
		bool
	)

Definition at line 384 of file EmissionVetoHook1.cc.

References accepted, gather_cfg::cout, emittedMode, fatalEmissionVeto(), i, j, gen::k, min, nAcceptSeq, nFSRveto, pTcalc(), pTemtMode, pThard, alignCSCRings::r, ntuplemaker::status, swap(), vetoCount, and vetoOn.

                                                                                     {
   // Must be radiation from the hard system
   if (iSys != 0) return false;
 
   // If we already have accepted 'vetoCount' emissions in a row, do nothing
   if (vetoOn && nAcceptSeq >= vetoCount) return false;
 
   // Pythia radiator (before and after), emitted and recoiler (after)
   int iRecAft = e.size() - 1;
   int iEmt    = e.size() - 2;
   int iRadAft = e.size() - 3;
   int iRadBef = e[iEmt].mother1();
   if ( (e[iRecAft].status() != 52 && e[iRecAft].status() != -53) ||
        e[iEmt].status() != 51 || e[iRadAft].status() != 51) {
     e.list();
     fatalEmissionVeto(string("Couldn't find Pythia FSR emission"));
   }
 
   // Behaviour based on pTemtMode:
   //  0 - pT of emitted w.r.t. radiator before
   //  1 - min(pT of emitted w.r.t. all incoming/outgoing)
   //  2 - min(pT of all outgoing w.r.t. all incoming/outgoing)
   int xSR = (pTemtMode == 0) ? 1       : -1;
   int i   = (pTemtMode == 0) ? iRadBef : -1;
   int k   = (pTemtMode == 0) ? iRadAft : -1;
   int r   = (pTemtMode == 0) ? iRecAft : -1;
 
   // When pTemtMode is 0 or 1, iEmt has been selected
   double pTemt = -1.;
   if (pTemtMode == 0 || pTemtMode == 1) {
     // Which parton is emitted, based on emittedMode:
     //  0 - Pythia definition of emitted
     //  1 - Pythia definition of radiated after emission
     //  2 - Random selection of emitted or radiated after emission
     //  3 - Try both emitted and radiated after emission
     int j = iRadAft;
     if (emittedMode == 0 || (emittedMode == 2 && rndmPtr->flat() < 0.5)) j++;
 
     for (int jLoop = 0; jLoop < 2; jLoop++) {
       if      (jLoop == 0) pTemt = pTcalc(e, i, j, k, r, xSR);
       else if (jLoop == 1) pTemt = min(pTemt, pTcalc(e, i, j, k, r, xSR));
   
       // For emittedMode == 3, have tried iRadAft, now try iEmt
       if (emittedMode != 3) break;
       if (k != -1) swap(j, k); else j = iEmt;
     }
 
   // If pTemtMode is 2, then try all final-state partons as emitted
   } else if (pTemtMode == 2) {
     pTemt = pTcalc(e, i, -1, k, r, xSR);
 
   }
 
 #ifdef DBGOUTPUT
   cout << "doVetoFSREmission: pTemt = " << pTemt << endl << endl;
 #endif
 
   // Veto if pTemt > pThard
   if (pTemt > pThard) {
     nAcceptSeq = 0;
     nFSRveto++;
     return true;
   }
 
   // Else mark that an emission has been accepted and continue
   nAcceptSeq++;
   accepted = true;
   return false;
 }

bool EmissionVetoHook1::doVetoISREmission	(	int	,
		const Pythia8::Event &	e,
		int	iSys
	)

Definition at line 333 of file EmissionVetoHook1.cc.

References accepted, gather_cfg::cout, fatalEmissionVeto(), i, j, gen::k, nAcceptSeq, nISRveto, pTcalc(), pTemtMode, pThard, alignCSCRings::r, ntuplemaker::status, vetoCount, and vetoOn.

                                                                               {
   // Must be radiation from the hard system
   if (iSys != 0) return false;
 
   // If we already have accepted 'vetoCount' emissions in a row, do nothing
   if (vetoOn && nAcceptSeq >= vetoCount) return false;
 
   // Pythia radiator after, emitted and recoiler after.
   int iRadAft = -1, iEmt = -1, iRecAft = -1;
   for (int i = e.size() - 1; i > 0; i--) {
     if      (iRadAft == -1 && e[i].status() == -41) iRadAft = i;
     else if (iEmt    == -1 && e[i].status() ==  43) iEmt    = i;
     else if (iRecAft == -1 && e[i].status() == -42) iRecAft = i;
     if (iRadAft != -1 && iEmt != -1 && iRecAft != -1) break;
   }
   if (iRadAft == -1 || iEmt == -1 || iRecAft == -1) {
     e.list();
     fatalEmissionVeto(string("Couldn't find Pythia ISR emission"));
   }
 
   // pTemtMode == 0: pT of emitted w.r.t. radiator
   // pTemtMode == 1: min(pT of emitted w.r.t. all incoming/outgoing)
   // pTemtMode == 2: min(pT of all outgoing w.r.t. all incoming/outgoing)
   int xSR      = (pTemtMode == 0) ? 0       : -1;
   int i        = (pTemtMode == 0) ? iRadAft : -1;
   int j        = (pTemtMode != 2) ? iEmt    : -1;
   int k        = -1;
   int r        = (pTemtMode == 0) ? iRecAft : -1;
   double pTemt = pTcalc(e, i, j, k, r, xSR);
 
 #ifdef DBGOUTPUT
   cout << "doVetoISREmission: pTemt = " << pTemt << endl << endl;
 #endif
 
   // Veto if pTemt > pThard
   if (pTemt > pThard) {
     nAcceptSeq = 0;
     nISRveto++;
     return true;
   }
 
   // Else mark that an emission has been accepted and continue
   nAcceptSeq++;
   accepted = true;
   return false;
 }

bool EmissionVetoHook1::doVetoMPIEmission	(	int	,
		const Pythia8::Event &	e
	)

Definition at line 458 of file EmissionVetoHook1.cc.

References gather_cfg::cout, MPIvetoOn, and pTMPI.

                                                                     {
   if (MPIvetoOn) {
     if (e[e.size() - 1].pT() > pTMPI) {
 #ifdef DBGOUTPUT
       cout << "doVetoMPIEmission: pTnow = " << e[e.size() - 1].pT()
            << ", pTMPI = " << pTMPI << endl << endl;
 #endif
       return true;
     }
   }
   return false;
 }

bool EmissionVetoHook1::doVetoMPIStep	(	int	nMPI,
		const Pythia8::Event &	e
	)

Definition at line 248 of file EmissionVetoHook1.cc.

References abs, accepted, prof2calltree::count, gather_cfg::cout, fatalEmissionVeto(), first, i, prof2calltree::last, MPIvetoOn, nAcceptSeq, nFinal, nFinalExt, pTcalc(), pThard, pThardMode, pTMPI, and Verbosity.

                                                                      {
   // Extra check on nMPI
   if (nMPI > 1) return false;
 
   // Find if there is a POWHEG emission. Go backwards through the
   // event record until there is a non-final particle. Also sum pT and
   // find pT_1 for possible MPI vetoing
   int count = 0, inonfinal = 0;
   double pT1 = 0., pTsum = 0.;
   for (int i = e.size() - 1; i > 0; i--) {
     inonfinal = i;
     if (e[i].isFinal()) {
       count++;
       pT1    = e[i].pT();
       pTsum += e[i].pT();
     } else break;
   }
 
   nFinal = nFinalExt;
 
   if (nFinal < 0) {      // nFinal is not specified from external, try to find out
     int first = -1, myid;
     int last = -1;
     for(int ip = 2; ip < e.size(); ip++) {
       myid = e[ip].id();
       if(abs(myid) < 6 || abs(myid) == 21) continue;
       first = ip;
       break;
     }
     if(first < 0) fatalEmissionVeto(string("signal particles not found"));
     for(int ip = first; ip < e.size(); ip++) {
       myid = e[ip].id();
       if(abs(myid) < 6 || abs(myid) == 21) continue;
       last = ip;
     }
     nFinal = last - inonfinal;
   }
 
   // Extra check that we have the correct final state
   if (count != nFinal && count != nFinal + 1)
     fatalEmissionVeto(string("Wrong number of final state particles in event"));
 
   // Flag if POWHEG radiation present and index
   bool isEmt = (count == nFinal) ? false : true;
   int  iEmt  = (isEmt) ? e.size() - 1 : -1;
 
   // If there is no radiation or if pThardMode is 0 then set pThard to QRen.
   if (!isEmt || pThardMode == 0) {
     pThard = infoPtr->QRen();
       
   // If pThardMode is 1 then the pT of the POWHEG emission is checked against
   // all other incoming and outgoing partons, with the minimal value taken
   } else if (pThardMode == 1) {
     pThard = pTcalc(e, -1, iEmt, -1, -1, -1);
 
   // If pThardMode is 2, then the pT of all final-state partons is checked
   // against all other incoming and outgoing partons, with the minimal value
   // taken
   } else if (pThardMode == 2) {
     pThard = pTcalc(e, -1, -1, -1, -1, -1);
 
   }
 
   // Find MPI veto pT if necessary
   if (MPIvetoOn) {
     pTMPI = (isEmt) ? pTsum / 2. : pT1;
   }
 
   if(Verbosity)
     cout << "doVetoMPIStep: QFac = " << infoPtr->QFac()
          << ", QRen = " << infoPtr->QRen()
          << ", pThard = " << pThard << endl << endl;
 
   // Initialise other variables
   accepted   = false;
   nAcceptSeq = 0;
 
   // Do not veto the event
   return false;
 }

void EmissionVetoHook1::fatalEmissionVeto ( string message )

Definition at line 4 of file EmissionVetoHook1.cc.

References edm::errors::Configuration, and edm::hlt::Exception.

Referenced by doVetoFSREmission(), doVetoISREmission(), and doVetoMPIStep().

                                                         {
   throw edm::Exception(edm::errors::Configuration,"Pythia8Interface")
     << "EmissionVeto: " << message << endl;
 }

int EmissionVetoHook1::numberVetoMPIStep ( )

inline

Definition at line 25 of file EmissionVetoHook1.h.

25 { return 1; }

double EmissionVetoHook1::pTcalc	(	const Pythia8::Event &	e,
		int	i,
		int	j,
		int	k,
		int	r,
		int	xSRin
	)

Definition at line 117 of file EmissionVetoHook1.cc.

References gather_cfg::cout, min, pTdefMode, pTpowheg(), and pTpythia().

Referenced by doVetoFSREmission(), doVetoISREmission(), and doVetoMPIStep().

                                                                                              {
 
   // Loop over ISR and FSR if necessary
   double pTemt = -1., pTnow;
   int xSR1 = (xSRin == -1) ? 0 : xSRin;
   int xSR2 = (xSRin == -1) ? 2 : xSRin + 1;
   for (int xSR = xSR1; xSR < xSR2; xSR++) {
     // FSR flag
     bool FSR = (xSR == 0) ? false : true;
 
     // If all necessary arguments have been given, then directly calculate.
     // POWHEG ISR and FSR, need i and j.
     if ((pTdefMode == 0 || pTdefMode == 1) && i > 0 && j > 0) {
       pTemt = pTpowheg(e, i, j, (pTdefMode == 0) ? false : FSR);
 
     // Pythia ISR, need i, j and r.
     } else if (!FSR && pTdefMode == 2 && i > 0 && j > 0 && r > 0) {
       pTemt = pTpythia(e, i, j, r, FSR);
 
     // Pythia FSR, need k, j and r.
     } else if (FSR && pTdefMode == 2 && j > 0 && k > 0 && r > 0) {
       pTemt = pTpythia(e, k, j, r, FSR);
 
     // Otherwise need to try all possible combintations.
     } else {
       // Start by finding incoming legs to the hard system after
       // branching (radiator after branching, i for ISR).
       // Use partonSystemsPtr to find incoming just prior to the
       // branching and track mothers.
       int iInA = partonSystemsPtr->getInA(0);
       int iInB = partonSystemsPtr->getInB(0);
       while (e[iInA].mother1() != 1) { iInA = e[iInA].mother1(); }
       while (e[iInB].mother1() != 2) { iInB = e[iInB].mother1(); }
 
       // If we do not have j, then try all final-state partons
       int jNow = (j > 0) ? j : 0;
       int jMax = (j > 0) ? j + 1 : e.size();
       for (; jNow < jMax; jNow++) {
 
         // Final-state and coloured jNow or photon only
         if (!e[jNow].isFinal()) continue;
         if (e[jNow].colType() == 0 && e[jNow].id() != 22) continue;
 
         // POWHEG
         if (pTdefMode == 0 || pTdefMode == 1) {
 
           // ISR - only done once as just kinematical pT
           if (!FSR) {
             pTnow = pTpowheg(e, iInA, jNow, (pTdefMode == 0) ? false : FSR);
             if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow);
   
           // FSR - try all outgoing partons from system before branching 
           // as i. Note that for the hard system, there is no 
           // "before branching" information.
           } else {
     
             int outSize = partonSystemsPtr->sizeOut(0);
             for (int iMem = 0; iMem < outSize; iMem++) {
               int iNow = partonSystemsPtr->getOut(0, iMem);
 
               // Coloured only, i != jNow and no carbon copies
               if (iNow == jNow || e[iNow].colType() == 0) continue;
               if (jNow == e[iNow].daughter1() 
                   && jNow == e[iNow].daughter2()) continue;
 
               pTnow = pTpowheg(e, iNow, jNow, (pTdefMode == 0) 
                 ? false : FSR);
               if (pTnow > 0.) pTemt = (pTemt < 0) 
                 ? pTnow : min(pTemt, pTnow);
             } // for (iMem)
   
           } // if (!FSR)
   
         // Pythia
         } else if (pTdefMode == 2) {
   
           // ISR - other incoming as recoiler
           if (!FSR) {
             pTnow = pTpythia(e, iInA, jNow, iInB, FSR);
             if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow);
             pTnow = pTpythia(e, iInB, jNow, iInA, FSR);
             if (pTnow > 0.) pTemt = (pTemt < 0) ? pTnow : min(pTemt, pTnow);
   
           // FSR - try all final-state coloured partons as radiator
           //       after emission (k).
           } else {
             for (int kNow = 0; kNow < e.size(); kNow++) {
               if (kNow == jNow || !e[kNow].isFinal() ||
                   e[kNow].colType() == 0) continue;
   
               // For this kNow, need to have a recoiler.
               // Try two incoming.
               pTnow = pTpythia(e, kNow, jNow, iInA, FSR);
               if (pTnow > 0.) pTemt = (pTemt < 0) 
                 ? pTnow : min(pTemt, pTnow);
               pTnow = pTpythia(e, kNow, jNow, iInB, FSR);
               if (pTnow > 0.) pTemt = (pTemt < 0) 
                 ? pTnow : min(pTemt, pTnow);
 
               // Try all other outgoing.
               for (int rNow = 0; rNow < e.size(); rNow++) {
                 if (rNow == kNow || rNow == jNow ||
                     !e[rNow].isFinal() || e[rNow].colType() == 0) continue;
                 pTnow = pTpythia(e, kNow, jNow, rNow, FSR);
                 if (pTnow > 0.) pTemt = (pTemt < 0) 
                   ? pTnow : min(pTemt, pTnow);
               } // for (rNow)
   
             } // for (kNow)
           } // if (!FSR)
         } // if (pTdefMode)
       } // for (j)
     }
   } // for (xSR)
 
 #ifdef DBGOUTPUT
   cout << "pTcalc: i = " << i << ", j = " << j << ", k = " << k
        << ", r = " << r << ", xSR = " << xSRin
        << ", pTemt = " << pTemt << endl;
 #endif
 
   return pTemt;
 }

double EmissionVetoHook1::pTpowheg	(	const Pythia8::Event &	e,
		int	i,
		int	j,
		bool	FSR
	)

Definition at line 74 of file EmissionVetoHook1.cc.

References gather_cfg::cout, alignCSCRings::e, j, AlCaHLTBitMon_ParallelJobs::p, and mathSSE::sqrt().

Referenced by pTcalc().

                                                                                 {
 
   // pT value for FSR and ISR
   double pTnow = 0.;
   if (FSR) {
     // POWHEG d_ij (in CM frame). Note that the incoming beams have not
     // been updated in the parton systems pointer yet (i.e. prior to any
     // potential recoil).
     int iInA = partonSystemsPtr->getInA(0);
     int iInB = partonSystemsPtr->getInB(0);
     double betaZ = - ( e[iInA].pz() + e[iInB].pz() ) /
                      ( e[iInA].e()  + e[iInB].e()  );
     Pythia8::Vec4 iVecBst(e[i].p()), jVecBst(e[j].p());
     iVecBst.bst(0., 0., betaZ);
     jVecBst.bst(0., 0., betaZ);
     pTnow = sqrt( (iVecBst + jVecBst).m2Calc() *
                   iVecBst.e() * jVecBst.e() /
                   Pythia8::pow2(iVecBst.e() + jVecBst.e()) );
 
   } else {
     // POWHEG pT_ISR is just kinematic pT
     pTnow = e[j].pT();
   }
 
   // Check result
   if (pTnow < 0.) {
     cout << "Warning: pTpowheg was negative" << endl;
     return -1.;
   }
 
 #ifdef DBGOUTPUT
   cout << "pTpowheg: i = " << i << ", j = " << j
        << ", pTnow = " << pTnow << endl;
 #endif
 
    return pTnow;
 }

double EmissionVetoHook1::pTpythia	(	const Pythia8::Event &	e,
		int	RadAfterBranch,
		int	EmtAfterBranch,
		int	RecAfterBranch,
		bool	FSR
	)

Definition at line 17 of file EmissionVetoHook1.cc.

References abs, gather_cfg::cout, mathSSE::sqrt(), and detailsBasic3DVector::z.

Referenced by pTcalc().

                                                                                      {
 
   // Convenient shorthands for later
   Pythia8::Vec4 radVec = e[RadAfterBranch].p();
   Pythia8::Vec4 emtVec = e[EmtAfterBranch].p();
   Pythia8::Vec4 recVec = e[RecAfterBranch].p();
   int  radID  = e[RadAfterBranch].id();
 
   // Calculate virtuality of splitting
   double sign = (FSR) ? 1. : -1.;
   Pythia8::Vec4 Q(radVec + sign * emtVec); 
   double Qsq = sign * Q.m2Calc();
 
   // Mass term of radiator
   double m2Rad = (abs(radID) >= 4 && abs(radID) < 7) ?
                   Pythia8::pow2(particleDataPtr->m0(radID)) : 0.;
 
   // z values for FSR and ISR
   double z, pTnow;
   if (FSR) {
     // Construct 2 -> 3 variables
     Pythia8::Vec4 sum = radVec + recVec + emtVec;
     double m2Dip = sum.m2Calc();
     double x1 = 2. * (sum * radVec) / m2Dip;
     double x3 = 2. * (sum * emtVec) / m2Dip;
     z     = x1 / (x1 + x3);
     pTnow = z * (1. - z);
 
   } else {
     // Construct dipoles before/after splitting
     Pythia8::Vec4 qBR(radVec - emtVec + recVec);
     Pythia8::Vec4 qAR(radVec + recVec);
     z     = qBR.m2Calc() / qAR.m2Calc();
     pTnow = (1. - z);
   }
 
   // Virtuality with correct sign
   pTnow *= (Qsq - sign * m2Rad);
 
   // Can get negative pT for massive splittings
   if (pTnow < 0.) {
     cout << "Warning: pTpythia was negative" << endl;
     return -1.;
   }
 
 #ifdef DBGOUTPUT
   cout << "pTpythia: rad = " << RadAfterBranch << ", emt = "
        << EmtAfterBranch << ", rec = " << RecAfterBranch
        << ", pTnow = " << sqrt(pTnow) << endl;
 #endif
 
   // Return pT
   return sqrt(pTnow);
 }