#include <PH_HEPEVT_Interface.h>

Static Public Member Functions
static int	set (PhotosBranch *branch)

static void	get ()

static void	prepare ()

static void	check_ME_channel ()

static void	complete ()

static void	clear ()

Static Public Attributes
static int	decay_idx =0

static int	ME_channel =0

Static Private Member Functions
static void	add_particle (int i, PhotosParticle *particle, int first_mother, int last_mother, int first_daughter, int last_daughter)

Static Private Attributes
static std::vector < PhotosParticle * >	m_particle_list

Detailed Description

Definition at line 40 of file PH_HEPEVT_Interface.h.

Member Function Documentation

void Photospp::PH_HEPEVT_Interface::add_particle	(	int	i,
		PhotosParticle *	particle,
		int	first_mother,
		int	last_mother,
		int	first_daughter,
		int	last_daughter
	)

staticprivate

Add single particle to HEPEVT

Definition at line 45 of file PH_HEPEVT_Interface.cxx.

References Photospp::PhotosParticle::getE(), Photospp::PhotosParticle::getMass(), Photospp::PhotosParticle::getPdgID(), Photospp::PhotosParticle::getPx(), Photospp::PhotosParticle::getPy(), Photospp::PhotosParticle::getPz(), Photospp::PhotosParticle::getStatus(), Photospp::PhotosParticle::getVirtuality(), i, Photospp::ph_hepevt_, and Photospp::ph_phoqed_.

                                                                              {
 
   if(i>0)
     i--; //account for fortran indicies begining at 1
   else
     Log::Warning()<<"Index given to PH_HEPEVT_Interface::add_particle "
                   <<"is too low (it must be > 0)."<<endl;
 
   //add to our internal list of pointer/index pairs
   m_particle_list.push_back(particle);
 
   //now set the element of PH_HEPEVT
   ph_hepevt_.nevhep=0; //dummy
   //ph_hepevt_.nhep=ph_hepevt_.nhep++; // [R.K. 9/2018] This line may lead to undefined behaviour.
   ph_hepevt_.nhep++;                   // [R.K. 9/2018] I assume the number should be increased by one
   ph_hepevt_.isthep[i]=particle->getStatus();
   ph_hepevt_.idhep[i]=particle->getPdgID();
 
   ph_hepevt_.jmohep[i][0]=first_mother;
   ph_hepevt_.jmohep[i][1]=last_mother;
 
   ph_hepevt_.jdahep[i][0]=first_daughter;
   ph_hepevt_.jdahep[i][1]=last_daughter;
 
   ph_hepevt_.phep[i][0]=particle->getPx();
   ph_hepevt_.phep[i][1]=particle->getPy();
   ph_hepevt_.phep[i][2]=particle->getPz();
   ph_hepevt_.phep[i][3]=particle->getE();
 
   // if massFrom4Vector=true (default) - get sqrt(e^2-p^2)
   // otherwise - get mass from event record
   if(!Photos::massFrom4Vector) ph_hepevt_.phep[i][4]=particle->getMass();
   else                         ph_hepevt_.phep[i][4]=particle->getVirtuality();
 
   int pdgid = abs(particle->getPdgID());
 
   // if 'forceMass' for this PDGID was used - overwrite mass
   if(Photos::forceMassList)
   {
     for(unsigned int j=0;j<Photos::forceMassList->size();j++)
     {
       if(pdgid == abs(Photos::forceMassList->at(j)->first))
       {
         double mass = Photos::forceMassList->at(j)->second;
 
         // when 'forceMass' is used the mass provided is larger than 0.0
         // when 'forceMassFromEventRecord' is used mass is -1.0
         // in this case - get mass from event record
         if(mass<0.0) mass = particle->getMass();
         ph_hepevt_.phep[i][4] = mass;
       }
     }
   }
 
   ph_hepevt_.vhep[i][0]=0;
   ph_hepevt_.vhep[i][1]=0;
   ph_hepevt_.vhep[i][2]=0;
   ph_hepevt_.vhep[i][3]=0;
 
   ph_phoqed_.qedrad[i]=1;
 
 }

void Photospp::PH_HEPEVT_Interface::check_ME_channel ( )

static

Check channel for complete matrix element calculation

Definition at line 419 of file PH_HEPEVT_Interface.cxx.

References i, and Photospp::ph_hepevt_.

 {
         ME_channel=0;
 
 // Check mothers:
 
         if(decay_idx==2)                              return; // Only one mother present
         if(ph_hepevt_.idhep[0]*ph_hepevt_.idhep[1]>0) return; // Mothers have same sign
 
         Log::Debug(900)<<"ME_channel: Mothers PDG:  "<<ph_hepevt_.idhep[0]<<" "<<ph_hepevt_.idhep[1]<<endl;
         if(decay_idx)
                 Log::Debug(900,false)<<"            Intermediate: "<<ph_hepevt_.idhep[decay_idx-1]<<endl;
 
         int              firstDaughter=3;
         if(decay_idx==0) firstDaughter=2; // if no intermediate particle - daughters start at idx 2
 
         // Are mothers in range +/- 1-6; +/- 11-16?
         int mother1 = abs(ph_hepevt_.idhep[0]);
         int mother2 = abs(ph_hepevt_.idhep[1]);
         if( mother1<1 || (mother1>6 && mother1<11) || mother1>16 ) return;
         if( mother2<1 || (mother2>6 && mother2<11) || mother2>16 ) return;
 
 //Check daughters
 
         // Z: check for pairs 11 -11 ; 13 -13 ; 15 -15
         // -------------------------------------------
         int firstPDG =0;
         int secondPDG=0;
         for(int i=firstDaughter; i<ph_hepevt_.nhep;i++)
         {
                 int pdg = abs(ph_hepevt_.idhep[i]);
                 if(pdg==11 || pdg==13 || pdg==15)
                 {
                         if(firstPDG==0) firstPDG=ph_hepevt_.idhep[i];
                         else
                         {
                                 secondPDG=ph_hepevt_.idhep[i];
                                 // Just in case two pairs are genereted - verify that we have a pair with oposite signs
                                 if(firstPDG*secondPDG>0) secondPDG=0;
                                 break;
                         }
                 }
         }
 
         if( ME_channel==0 && firstPDG!=0 && secondPDG!=0 &&
             firstPDG==-secondPDG ) ME_channel=1;
 
         // W: check for pairs 11 -12; -11 12; 13 -14; -13 14; 15 -16; -15 16
         // -----------------------------------------------------------------
         firstPDG =0;
         secondPDG=0;
         for(int i=firstDaughter; i<ph_hepevt_.nhep;i++)
         {
                 int pdg = abs(ph_hepevt_.idhep[i]);
                 if(pdg>=11 && pdg<=16)
                 {
                         if(firstPDG==0) firstPDG=ph_hepevt_.idhep[i];
                         else
                         {
                                 secondPDG=ph_hepevt_.idhep[i];
                                 // Just in case two pairs are genereted - verify that we have a pair with oposite signs
                                 if(firstPDG*secondPDG>0) secondPDG=0;
                                 break;
                         }
                 }
         }
 
         firstPDG =abs(firstPDG);
         secondPDG=abs(secondPDG);
 
         if(  ME_channel==0 && firstPDG!=0 && secondPDG!=0 &&
            ( ( firstPDG==11 && secondPDG==12 ) || (firstPDG == 12 && secondPDG == 11) ||
              ( firstPDG==13 && secondPDG==14 ) || (firstPDG == 14 && secondPDG == 13) ||
              ( firstPDG==15 && secondPDG==16 ) || (firstPDG == 16 && secondPDG == 15)
            )
           ) ME_channel=2;
 
         Log::Debug(901)<<"ME_channel: Found ME_channel: "<<ME_channel<<endl;
 
 // Check intermediate particle (if exists):
 
         // Verify that intermediate particle PDG matches ME_channel found
         if(ME_channel>0 && decay_idx)
         {
                 int pdg=ph_hepevt_.idhep[decay_idx-1];
 
                 if(ME_channel==1 && !(pdg==22 || pdg==23) ) ME_channel=0; //gamma/Z
                 if(ME_channel==2 && !(pdg==24 || pdg==-24)) ME_channel=0; //W+/W-
 
                 if(ME_channel==0)
                         Log::Debug(901,false)<<"            but set to 0: wrong intermediate particle: "<<pdg<<endl;
         }
 
 // Check flags
 
         switch(ME_channel)
         {
                 case  0: break; // Ok - no channel found
                 case  1: if(!Photos::meCorrectionWtForZ) ME_channel=0; break;
                 case  2: if(!Photos::meCorrectionWtForW) ME_channel=0; break;
                 default: Log::Error()<<"Unimplemented ME channel: "<<ME_channel<<endl; break;
         }
         Log::Debug(902)<<"ME_channel: Finally, after flag check, ME_channel is: "<<ME_channel<<endl;
 }

void Photospp::PH_HEPEVT_Interface::clear ( void )

static

Clear HEPEVT

for(int i=0; i < NMXHEP; i++){

ph_hepevt_.isthep[i]=0; ph_hepevt_.idhep[i]=0;

for(int j=0; j<2; j++){ ph_hepevt_.jmohep[i][j]=0; ph_hepevt_.jdahep[i][j]=0; }

for(int j=0; j<5; j++) ph_hepevt_.phep[i][j]=0;

for(int j=0; j<4; j++) ph_hepevt_.vhep[i][j]=0;

ph_phoqed_.qedrad[i]=0;

}

Definition at line 16 of file PH_HEPEVT_Interface.cxx.

References Photospp::ph_hepevt_.

                                {
 
   m_particle_list.clear();
 
   ph_hepevt_.nevhep=0;
   ph_hepevt_.nhep=0;
 
 }

void Photospp::PH_HEPEVT_Interface::complete ( )

static

Finalize processing

Definition at line 414 of file PH_HEPEVT_Interface.cxx.

Referenced by Photospp::PhotosBranch::process().

 {
 
 }

void Photospp::PH_HEPEVT_Interface::get ( )

static

Update event record with data from HEPEVT

Definition at line 161 of file PH_HEPEVT_Interface.cxx.

Referenced by Photospp::PhotosBranch::process().

                              {
 
   int index = 0;
 
   //if no photons have been added to the event record, do nothing.
   if(ph_hepevt_.nhep == (int) m_particle_list.size())
     return;
 
   //phodmp_();
 
   int  particle_count  = m_particle_list.size();
   int  daughters_start = ph_hepevt_.jmohep[ph_hepevt_.nhep-1][0];
   int  photons         = ph_hepevt_.nhep - m_particle_list.size();
   bool isPhotonCreated = (photons>0);
 
   std::vector<PhotosParticle*> photon_list; // list of added photons
                                             // which need kinematical treatment
                                             // in special case
 
   // we decipher daughters_start from  last entry
   // that is last daughter in  ph_hepevt_
   // another option of this functionality may be
   // ph_hepevt_.jdahep[ ph_hepevt_.jmohep[ph_hepevt_.nhep-1][0]-1][0];
   // Update daughters_start if there are two mothers
   // NOTE: daughters_start is index for C++ arrays, while ph_hepevt_.jmohep
   //       contains indices for Fortran arrays.
   if(ph_hepevt_.jmohep[ph_hepevt_.nhep-1][1]>0)
     daughters_start = ph_hepevt_.jmohep[ph_hepevt_.nhep-1][1];
 
   index = particle_count;
 
   // Add extra photons
   for(;photons>0; photons--, index++){
 
     if(ph_hepevt_.idhep[index]!=PhotosParticle::GAMMA)
       Log::Fatal("PH_HEPEVT_Interface::get(): Extra particle added to the PH_HEPEVT common block in not a photon!",6);
 
     //create a new particle
     PhotosParticle * new_photon;
     new_photon = m_particle_list.at(0)->createNewParticle(ph_hepevt_.idhep[index],
                                                           ph_hepevt_.isthep[index],
                                                           ph_hepevt_.phep[index][4],
                                                           ph_hepevt_.phep[index][0],
                                                           ph_hepevt_.phep[index][1],
                                                           ph_hepevt_.phep[index][2],
                                                           ph_hepevt_.phep[index][3]);
 
     //add into the event record
     //get mother particle of photon
     PhotosParticle * mother =  m_particle_list.at(ph_hepevt_.jmohep[index][0]-1);
     mother->addDaughter(new_photon);
 
     //add to list of photons
     photon_list.push_back(new_photon);
   }
 
   // Before we update particles, we check for special cases
   // At this step, particles are yet unmodified
   // but photons are already in the event record
   bool special=false;
   PhotosParticle *p1 = NULL;
   PhotosParticle *p2 = NULL;
 
   if( isPhotonCreated )
   {
     std::vector<PhotosParticle*> daughters;
 
     // in the following we create list of   daughters,
     // later  we calculate bool special which is true only if all
     // daughters self-decay
     // at peresent warning for  mixed self-decay and not self decay
     // daughters is not printed.
 
     for(int i=daughters_start;i<particle_count;i++)
     {
       PhotosParticle *p = m_particle_list.at(i);
 
       daughters.push_back(p);
     }
 
     // Check if this is a special case
     special = true;
 
     if(daughters.size()==0) special = false;
 
     // special = false if there is a stable particle on the list
     //                 or there is a particle without self-decay
     for(unsigned int i=0;i<daughters.size();i++)
     {
       if(daughters[i]->getStatus()==1)
       {
         special = false;
         break;
       }
 
       // NOTE: We can use 'getDaughters' here, because vertices
       //       of daughters are not being modified by Photos right now
       //       (so there will be no caching)
       std::vector<PhotosParticle*> daughters2 = daughters[i]->getDaughters();
 
       if(daughters2.size()!=1 ||
          daughters2[0]->getPdgID() != daughters[i]->getPdgID() )
       {
         special = false;
         break;
       }
     }
 
     if( special )
     {
       double px1=0.0, py1=0.0, pz1=0.0, e1=0.0;
       double px2=0.0, py2=0.0, pz2=0.0, e2=0.0;
 
       // get sum of 4-momenta of unmodified particles
       for(unsigned int i=0;i<daughters.size();i++)
       {
         // ignore photons
         if(daughters[i]->getPdgID()==22) continue;
 
         px1+=daughters[i]->getPx();
         py1+=daughters[i]->getPy();
         pz1+=daughters[i]->getPz();
         e1 +=daughters[i]->getE();
       }
 
       // get sum of 4-momenta of particles in self-decay vertices
       for(unsigned int i=0;i<daughters.size();i++)
       {
         // ignore photons
         if(daughters[i]->getPdgID()==22) continue;
 
         // since 'allDaughtersSelfDecay()' is true
         // each of these particles has exactly one daughter
         px2 += daughters[i]->getDaughters().at(0)->getPx();
         py2 += daughters[i]->getDaughters().at(0)->getPy();
         pz2 += daughters[i]->getDaughters().at(0)->getPz();
         e2  += daughters[i]->getDaughters().at(0)->getE();
       }
 
       //cout<<"ORIG: "<<px1<<" "<<py1<<" "<<pz1<<" "<<e1<<endl;
       //cout<<"SELF: "<<px2<<" "<<py2<<" "<<pz2<<" "<<e2<<endl;
 
       p1 = m_particle_list.at(0)->createNewParticle(0,-1,0.0,px1,py1,pz1,e1);
       p2 = m_particle_list.at(0)->createNewParticle(0,-2,0.0,px2,py2,pz2,e2);
 
       // Finaly, boost photons to appropriate frame
       for(unsigned int i=0;i<photon_list.size();i++)
       {
         PhotosParticle *boosted = photon_list[i]->createNewParticle( 22, 1,
                                     0.0,
                                     photon_list[i]->getPx(),
                                     photon_list[i]->getPy(),
                                     photon_list[i]->getPz(),
                                     photon_list[i]->getE()   );
 
         boosted->boostToRestFrame(p1);
         boosted->boostFromRestFrame(p2);
 
         photon_list[i]->createSelfDecayVertex(boosted);
 
         delete boosted;
       }
 
       Log::Warning()<<"Hidden interaction, all daughters self decay."
           <<"Potentially over simplified solution applied."<<endl;
     }
   }
 
   //otherwise loop over particles which are already in the
   //event record and modify their 4 momentum
   //4.03.2012: Fix to prevent kinematical trap in vertex of simultaneous:
   //           z-collinear and non-conservation pf E,p for dauthters of grandmothers
   for(index=daughters_start; index < particle_count && index < (int) m_particle_list.size(); index++){
 
     PhotosParticle * particle = m_particle_list.at(index);
 
     if(ph_hepevt_.idhep[index]!=particle->getPdgID())
       Log::Fatal("PH_HEPEVT_Interface::get(): Something is wrong with the PH_HEPEVT common block",5);
 
     // If photons were added - for each daughter create a history entry
     if(isPhotonCreated && Photos::isCreateHistoryEntries)
     {
       particle->createHistoryEntry();
     }
 
     //check to see if this particle's 4-momentum has been modified
     bool   update=false;
 
     // don't update particle if difference lower than THRESHOLD * particle energy (default threshold = 10e-8)
     double threshold = NO_BOOST_THRESHOLD*ph_hepevt_.phep[index][3];
     if( fabs(ph_hepevt_.phep[index][0]-particle->getPx()) > threshold ||
         fabs(ph_hepevt_.phep[index][1]-particle->getPy()) > threshold ||
         fabs(ph_hepevt_.phep[index][2]-particle->getPz()) > threshold ||
         fabs(ph_hepevt_.phep[index][3]-particle->getE())  > threshold    ) update=true;
 
     if(update)
     {
 
       //modify this particle's momentum and it's daughters momentum
       //Steps 1., 2. and 3. must be executed in order.
 
       //1. boost the particles daughters into it's (old) rest frame
       particle->boostDaughtersToRestFrame(particle);
 
       //2. change this particles 4 momentum
       particle->setPx(ph_hepevt_.phep[index][0]);
       particle->setPy(ph_hepevt_.phep[index][1]);
       particle->setPz(ph_hepevt_.phep[index][2]);
       particle->setE(ph_hepevt_.phep[index][3]);
 
       //3. boost the particles daughters back into the lab frame
       particle->boostDaughtersFromRestFrame(particle);
 
       if(special && particle->getDaughters().size()>0){
 
       // Algorithm for special case:
       // a. get self-daughter of 'particle'
       PhotosParticle *particled = particle->getDaughters().at(0);
 
       // b. boost 'particled' daughters to rest frame
       particled->boostDaughtersToRestFrame(particled);
 
       // c. copy four momentum of 'particle' into four momentum of
       //    its self-daughter 'particled'
 
       particled->setPx( particle->getPx() );
       particled->setPy( particle->getPy() );
       particled->setPz( particle->getPz() );
       particled->setE ( particle->getE()  );
 
       // d. boost self daughter to rest-frame of <e1>
       //    boost self daughter from rest-frame of <e2>
 
       particled->boostToRestFrame(p1);
       particled->boostFromRestFrame(p2);
 
       // e. boost the 'particled' daughters back into the lab frame
       particled->boostDaughtersFromRestFrame(particled);
       }
 
     }
   }
 
   // cleanup
   if(p1) delete p1;
   if(p2) delete p2;
 }

void Photospp::PH_HEPEVT_Interface::prepare ( )

static

Prepare particles for processing

Definition at line 409 of file PH_HEPEVT_Interface.cxx.

Referenced by Photospp::PhotosBranch::process().

 {
         check_ME_channel();
 }

int Photospp::PH_HEPEVT_Interface::set ( PhotosBranch * branch )

static

Convert PhotosBranch to HEPEVT

Definition at line 110 of file PH_HEPEVT_Interface.cxx.

References Photospp::PhotosBranch::getDaughters(), Photospp::PhotosBranch::getDecayingParticle(), Photospp::PhotosBranch::getMothers(), i, and idx.

Referenced by Photospp::PhotosBranch::process().

 {
         PH_HEPEVT_Interface::clear();
         int idx=1;
 
         //get mothers
         vector<PhotosParticle *> mothers = branch->getMothers();
         int nmothers=mothers.size();
 
         //check if mid-particle exist
         decay_idx=0;
         PhotosParticle *decay_particle = branch->getDecayingParticle();
         if(decay_particle) decay_idx=nmothers+1;
 
         //get daughters
         vector<PhotosParticle *> daughters = branch->getDaughters();
         int ndaughters=daughters.size();
 
         for(int i=0;i<nmothers;i++)
         {
                 if(decay_idx)
                         add_particle(idx++,mothers.at(i),
                                      0,0, //mothers
                                      decay_idx,decay_idx); //daughters
                 else
                         add_particle(idx++,mothers.at(i),
                                      0,0, //mothers
                                      nmothers+1,nmothers+ndaughters); //daughters
         }
 
         if(decay_particle)
                 add_particle(idx++,decay_particle,
                              1,nmothers, //mothers
                              nmothers+2,nmothers+1+ndaughters); //daughters
 
         for(int i=0;i<ndaughters;i++)
         {
                 if(decay_idx)
                         add_particle(idx++,daughters.at(i),
                                      decay_idx,decay_idx, //mothers
                                      0,0); //daughters
                 else
                         add_particle(idx++,daughters.at(i),
                                      1,nmothers, //mothers
                                      0,0); //daughters
         }
         //Log::RedirectOutput( phodmp_ , Log::Debug(1000) );
         Log::Debug(1000,false)<<"PH_HEPEVT returning: "<<( (decay_idx) ? decay_idx : 1 )<<" from "<<idx-1<<" particles."<<endl;
         return (decay_idx) ? decay_idx : 1;
 }

Member Data Documentation

int Photospp::PH_HEPEVT_Interface::decay_idx =0

static

Index of decaying particle

Definition at line 62 of file PH_HEPEVT_Interface.h.

vector< PhotosParticle * > Photospp::PH_HEPEVT_Interface::m_particle_list

staticprivate

List of particles added to HEPEVT

Definition at line 72 of file PH_HEPEVT_Interface.h.

int Photospp::PH_HEPEVT_Interface::ME_channel =0

static

Number of channel to be used - flag for fortran routines

Definition at line 64 of file PH_HEPEVT_Interface.h.

The documentation for this class was generated from the following files:

Static Public Member Functions

Static Public Attributes

Static Private Member Functions

Static Private Attributes

Detailed Description

Member Function Documentation

Member Data Documentation