#include <PhotosInterface.h>

Classes
struct	Scaling

Public Member Functions
HepMC::GenEvent *	apply (HepMC::GenEvent *)

void	avoidTauLeptonicDecays ()

void	configureOnlyFor (int)

void	init ()

bool	isTauLeptonicDecay (HepMC::GenVertex *)

	PhotosInterface ()

	PhotosInterface (const edm::ParameterSet &)

const std::vector< std::string > &	specialSettings ()

	~PhotosInterface ()

Private Member Functions
void	applyToBranch (HepMC::GenEvent *, int)

void	applyToVertex (HepMC::GenEvent *, int)

void	attachParticles (HepMC::GenEvent , HepMC::GenVertex , int)

Private Attributes
bool	fAvoidTauLeptonicDecays

std::vector< int >	fBarcodes

bool	fIsInitialized

int	fOnlyPDG

std::vector< int >	fSecVtxStore

std::vector< std::string >	fSpecialSettings

Detailed Description

Definition at line 20 of file PhotosInterface.h.

Constructor & Destructor Documentation

PhotosInterface::PhotosInterface ( )

Definition at line 43 of file PhotosInterface.cc.

References fAvoidTauLeptonicDecays, fIsInitialized, and fSpecialSettings.

    : fOnlyPDG(-1)
 {
    fSpecialSettings.push_back("QED-brem-off:all");
    fAvoidTauLeptonicDecays = false;
    fIsInitialized = false; 
 }

PhotosInterface::PhotosInterface ( const edm::ParameterSet & )

Definition at line 51 of file PhotosInterface.cc.

References fIsInitialized, and fSpecialSettings.

    : fOnlyPDG(-1)
 {
    fSpecialSettings.push_back("QED-brem-off:all");
    fIsInitialized = false;
 }

gen::PhotosInterface::~PhotosInterface ( )

inline

Definition at line 27 of file PhotosInterface.h.

27 {}

Member Function Documentation

HepMC::GenEvent * PhotosInterface::apply ( HepMC::GenEvent * evt )

Definition at line 83 of file PhotosInterface.cc.

References funct::abs(), applyToBranch(), applyToVertex(), fAvoidTauLeptonicDecays, fIsInitialized, fOnlyPDG, fSecVtxStore, isTauLeptonicDecay(), and phoqed_.

Referenced by heavyIonTools.ConfigureHeavyIons::__call__(), editorTools.UserCodeTool::__call__(), HiCoreTools.RestrictInputToAOD::__call__(), coreTools.RunOnData::__call__(), trackTools.MakeAODTrackCandidates::__call__(), metTools.AddMETCollection::__call__(), heavyIonTools.ProductionDefaults::__call__(), editorTools.ChangeSource::__call__(), HiCoreTools.RemoveMCMatching::__call__(), jetTools.AddJetCollection::__call__(), cmsswVersionTools.PickRelValInputFiles::__call__(), coreTools.RemoveMCMatching::__call__(), trackTools.MakePATTrackCandidates::__call__(), trigTools.SwitchOnTrigger::__call__(), heavyIonTools.SelectionDefaults::__call__(), HiCoreTools.RemoveAllPATObjectsBut::__call__(), heavyIonTools.DisbaleMonteCarloDeps::__call__(), HiCoreTools.RemoveSpecificPATObjects::__call__(), trigTools.SwitchOnTriggerStandAlone::__call__(), trackTools.MakeTrackCandidates::__call__(), tauTools.AddTauCollection::__call__(), trigTools.SwitchOnTriggerMatching::__call__(), HiCoreTools.RemoveCleaning::__call__(), HiCoreTools.AddCleaning::__call__(), trigTools.SwitchOnTriggerMatchingStandAlone::__call__(), jetTools.SwitchJetCollection::__call__(), trigTools.SwitchOnTriggerMatchEmbedding::__call__(), jetTools.AddJetID::__call__(), jetTools.SetTagInfos::__call__(), metUncertaintyTools.RunMEtUncertainties::__call__(), GenMuonRadiationAlgorithm::compFSR(), and gen::ExternalDecayDriver::decay().

 {
      
    if ( !fIsInitialized ) return evt; // conv.read_next_event();
       
    // loop over HepMC::GenEvent, find vertices
       
    // for ( int ip=0; ip<evt->particles_size(); ip++ )
    for ( int ip=0; ip<4000; ip++ ) // 4000 is the max size of the array
    {
       phoqed_.qedrad[ip]=true;
    }
       
    //
    // now do actual job
    //
    
    for ( int iv=1; iv<=evt->vertices_size(); iv++ )
    {
       
       bool legalVtx = false;
       
       fSecVtxStore.clear();
       
       HepMC::GenVertex* vtx = evt->barcode_to_vertex( -iv ) ;
       
       if ( vtx->particles_in_size() != 1 ) continue; // more complex than we need
       if ( vtx->particles_out_size() <= 1 ) continue; // no outcoming particles
       
       if ( (*(vtx->particles_in_const_begin()))->pdg_id() == 111 ) continue; // pi0 decay vtx - no point to try
       
       if ( fOnlyPDG != 1 && (*(vtx->particles_in_const_begin()))->pdg_id() != fOnlyPDG )
       {
          continue;
       }
       else
       {
          // requested for specific PDG ID only, typically tau (15)
      //
      // first check if a brem vertex, where outcoming are the same pdg id and a photon
      //
      bool same = false;
      for ( HepMC::GenVertex::particle_iterator pitr=vtx->particles_begin(HepMC::children);
                pitr != vtx->particles_end(HepMC::children); ++pitr)
          {
         if ( (*pitr)->pdg_id() == fOnlyPDG )
         {
            same = true;
            break;
         }
      }
      if ( same ) continue;
      
      // OK, we get here if incoming fOnlyPDG and something else outcoming
      // call it for the whole branch starting at vtx barcode iv, and go on
      // NOTE: theoretically, it has a danger of double counting in vertices
      // down the decay branch originating from fOnlyPDG, but in practice
      // it's unlikely that down the branchg there'll be more fOnlyPDG's
      
      // cross-check printout
      // vtx->print();
      
      // overprotection...
      //
      if ( fOnlyPDG == 15 && fAvoidTauLeptonicDecays && isTauLeptonicDecay( vtx ) ) continue; 
      
      applyToBranch( evt, -iv );
      continue;
       }
       
       // configured for all types of particles
       //
       for ( HepMC::GenVertex::particle_iterator pitr=vtx->particles_begin(HepMC::children);
             pitr != vtx->particles_end(HepMC::children); ++pitr) 
       {
 
      // quark or gluon out of this vertex - no good !
      if ( abs((*pitr)->pdg_id()) >=1 &&  abs((*pitr)->pdg_id()) <=8 ) break;
      if ( abs((*pitr)->pdg_id()) == 21 ) break;
 
          if ( (*pitr)->status() == 1 || (*pitr)->end_vertex() )
      {
         // OK, legal already !
         legalVtx = true;
         break;
      }
       }
       
       if ( !legalVtx ) continue;
       
       applyToVertex( evt, -iv );
 
    } // end of master loop
    
    // restore event number in HEPEVT (safety measure, somehow needed by Hw6)
    HepMC::HEPEVT_Wrapper::set_event_number( evt->event_number() );
 
    return evt;
       
 }

void PhotosInterface::applyToBranch	(	HepMC::GenEvent *	evt,
		int	vtxbcode
	)

private

Definition at line 276 of file PhotosInterface.cc.

References applyToVertex(), and fSecVtxStore.

Referenced by apply().

 {
   
    
    fSecVtxStore.clear();
       
    // 1st level vertex
    //
    applyToVertex( evt, vtxbcode );
    
    // now look down the branch for more vertices, if any  
    //
    // Note: fSecVtxStore gets filled up in applyToVertex, if necessary 
    //
    unsigned int vcounter = 0;  
    
    while ( vcounter < fSecVtxStore.size() )
    { 
       applyToVertex( evt, fSecVtxStore[vcounter] );
       vcounter++;
    }
    
    fSecVtxStore.clear(); 
 
    return;
 
 }

void PhotosInterface::applyToVertex	(	HepMC::GenEvent *	evt,
		int	vtxbcode
	)

private

Definition at line 184 of file PhotosInterface.cc.

References attachParticles(), fAvoidTauLeptonicDecays, fBarcodes, fSecVtxStore, getHLTprescales::index, isTauLeptonicDecay(), and photos_().

Referenced by apply(), and applyToBranch().

 {
 
    HepMC::GenVertex* vtx = evt->barcode_to_vertex( vtxbcode );
    
    if ( fAvoidTauLeptonicDecays && isTauLeptonicDecay( vtx ) ) return;
 
    // cross-check printout
    //
    // vtx->print();
             
    // first, flush out HEPEVT & tmp barcode storage
    //
    HepMC::HEPEVT_Wrapper::zero_everything();
    fBarcodes.clear();
       
    // add incoming particle
    //
    int index = 1;      
    HepMC::HEPEVT_Wrapper::set_id( index, (*(vtx->particles_in_const_begin()))->pdg_id() );
    HepMC::FourVector vec4;
    vec4 = (*(vtx->particles_in_const_begin()))->momentum();
    HepMC::HEPEVT_Wrapper::set_momentum( index, vec4.x(), vec4.y(), vec4.z(), vec4.e() );
    HepMC::HEPEVT_Wrapper::set_mass( index, (*(vtx->particles_in_const_begin()))->generated_mass() );
    HepMC::HEPEVT_Wrapper::set_position( index, vtx->position().x(), vtx->position().y(),
                                                   vtx->position().z(), vtx->position().t() );
    HepMC::HEPEVT_Wrapper::set_status( index, (*(vtx->particles_in_const_begin()))->status() );
    HepMC::HEPEVT_Wrapper::set_parents( index, 0, 0 );
    fBarcodes.push_back( (*(vtx->particles_in_const_begin()))->barcode() );
                              
    // add outcoming particles (decay products)
    //
    int lastDau = 1;
    for ( HepMC::GenVertex::particle_iterator pitr=vtx->particles_begin(HepMC::children);
             pitr != vtx->particles_end(HepMC::children); ++pitr) 
    {
 
      if ( (*pitr)->status() == 1 || (*pitr)->end_vertex() )
      {
         index++;
         vec4 = (*pitr)->momentum();
         HepMC::HEPEVT_Wrapper::set_id( index, (*pitr)->pdg_id() );
             HepMC::HEPEVT_Wrapper::set_momentum( index, vec4.x(), vec4.y(), vec4.z(), vec4.e() );
         HepMC::HEPEVT_Wrapper::set_mass( index, (*pitr)->generated_mass() );
         vec4 = (*pitr)->production_vertex()->position();
             HepMC::HEPEVT_Wrapper::set_position( index, vec4.x(), vec4.y(), vec4.z(), vec4.t() );
         HepMC::HEPEVT_Wrapper::set_status( index, (*pitr)->status() );
         HepMC::HEPEVT_Wrapper::set_parents( index, 1, 1 );
         fBarcodes.push_back( (*pitr)->barcode() );
         lastDau++;
      }
          if ( (*pitr)->end_vertex() )
          {
             fSecVtxStore.push_back( (*pitr)->end_vertex()->barcode() );
          }      
    }
             
    // store, further to set NHEP in HEPEVT
    //
    int nentries = index;
       
    // reset master pointer to mother
    index = 1;
    HepMC::HEPEVT_Wrapper::set_children ( index, 2, lastDau ); // FIXME: need to check 
                                                               // if last daughter>=2 !!!
       
    // finally, set number of entries (NHEP) in HEPEVT
    //
    HepMC::HEPEVT_Wrapper::set_number_entries( nentries );
 
    // cross-check printout HEPEVT
    //  HepMC::HEPEVT_Wrapper::print_hepevt();
      
    // OK, 1-level vertex is formed - now, call PHOTOS
    //
    photos_( index ) ;
       
    // another cross-check printout HEPEVT - after photos
    // HepMC::HEPEVT_Wrapper::print_hepevt();
 
 
    // now check if something has been generated 
    // and make all adjustments to underlying vtx/parts
    //
    attachParticles( evt, vtx, nentries );
 
    // ugh, done with this vertex !
 
    return;
 
 }

void PhotosInterface::attachParticles	(	HepMC::GenEvent *	evt,
		HepMC::GenVertex *	vtx,
		int	nentries
	)

private

Definition at line 304 of file PhotosInterface.cc.

References alignCSCRings::e, fBarcodes, configurableAnalysis::GenParticle, m, and ntuplemaker::status.

Referenced by applyToVertex().

 {
 
    if ( HepMC::HEPEVT_Wrapper::number_entries() > nentries ) 
    {
          // yes, need all corrections and adjustments -
      // figure out how many photons and what particles in 
      // the decay branch have changes;
      // also, follow up each one and correct accordingly;
      // at the same time, add photon(s) to the GenVertex
      //  
      
      // vtx->print();
      
      int largestBarcode = -1;
      int Nbcodes = fBarcodes.size();
      
      for ( int ip=1; ip<Nbcodes; ip++ )
      {
 
         int bcode = fBarcodes[ip];
         HepMC::GenParticle* prt = evt->barcode_to_particle( bcode );
         if ( bcode > largestBarcode ) largestBarcode = bcode;
         double px = HepMC::HEPEVT_Wrapper::px(ip+1);
         double py = HepMC::HEPEVT_Wrapper::py(ip+1);
         double pz = HepMC::HEPEVT_Wrapper::pz(ip+1);
         double e  = HepMC::HEPEVT_Wrapper::e(ip+1);
         double m  = HepMC::HEPEVT_Wrapper::m(ip+1);   
                 
         if ( prt->end_vertex() )
         {
            
            HepMC::GenVertex* endVtx = prt->end_vertex();
            
                std::vector<int> secVtxStorage;
                secVtxStorage.clear();
            
            secVtxStorage.push_back( endVtx->barcode() );
         
            HepMC::FourVector mom4 = prt->momentum();
         
            // now rescale all descendants
            double bet1[3], bet2[3], gam1, gam2, pb;
            double mass = mom4.m();
            bet1[0] = -(mom4.px()/mass);
            bet1[1] = -(mom4.py()/mass);
            bet1[2] = -(mom4.pz()/mass);
            bet2[0] = px/m;
            bet2[1] = py/m;
            bet2[2] = pz/m;
            gam1 = mom4.e()/mass;
            gam2 = e/m;
         
            unsigned int vcounter = 0;
                 
                while ( vcounter < secVtxStorage.size() )
            {
               
           HepMC::GenVertex* theVtx = evt->barcode_to_vertex( secVtxStorage[vcounter] );
                       
           for ( HepMC::GenVertex::particle_iterator pitr=theVtx->particles_begin(HepMC::children);
                         pitr != theVtx->particles_end(HepMC::children); ++pitr) 
                   {
            
              if ( (*pitr)->end_vertex() )
              {
                 secVtxStorage.push_back( (*pitr)->end_vertex()->barcode() );
              }
              
              if ( theVtx->particles_out_size() == 1 && (*pitr)->pdg_id() == prt->pdg_id() ) 
              {
                 // carbon copy
             (*pitr)->set_momentum( HepMC::FourVector(px,py,pz,e) );
             continue;
              }
                        
                  HepMC::FourVector dmom4 = (*pitr)->momentum();
            
                  // Boost vector to parent rest frame...
                  pb = bet1[0]*dmom4.px() + bet1[1]*dmom4.py() + bet1[2]*dmom4.pz();
                  double dpx = dmom4.px() + bet1[0] * (dmom4.e() + pb/(gam1+1.) );
                  double dpy = dmom4.py() + bet1[1] * (dmom4.e() + pb/(gam1+1.) );
                  double dpz = dmom4.pz() + bet1[2] * (dmom4.e() + pb/(gam1+1.) );
                  double de  = gam1*dmom4.e() + pb;
                  // ...and boost back to modified parent frame
                  pb = bet2[0]*dpx + bet2[1]*dpy + bet2[2]*dpz;
                  dpx += bet2[0] * ( de + pb/(gam2+1.) );
                  dpy += bet2[1] * ( de + pb/(gam2+1.) );
                  dpz += bet2[2] * ( de + pb/(gam2+1.) );
                  de *= gam2;
                  de += pb;
            
                  (*pitr)->set_momentum( HepMC::FourVector(dpx,dpy,dpz,de) );
              
               }               
           vcounter++;       
                }
            
            secVtxStorage.clear();
         }
         
         prt->set_momentum( HepMC::FourVector(px,py,pz,e) );
         
      } // ok, all affected particles update, but the photon(s) still not inserted
          
      int newlyGen =  HepMC::HEPEVT_Wrapper::number_entries() - nentries;
      
      if ( largestBarcode < evt->particles_size() )
      {
         // need to adjust barcodes down from the affected vertex/particles
         // such that we "free up" barcodes for newly generated photons
         // in the middle of the event record
         //
         for ( int ipp=evt->particles_size(); ipp>largestBarcode; ipp-- )
         {
            (evt->barcode_to_particle(ipp))->suggest_barcode( ipp+newlyGen );
         }
      }
      
      // now attach new generated photons to the vertex
      //
      for ( int ipnw=1; ipnw<=newlyGen; ipnw++ )
      {
          int nbcode = largestBarcode+ipnw;
          int pdg_id = HepMC::HEPEVT_Wrapper::id( nentries+ipnw );
          int status = HepMC::HEPEVT_Wrapper::status( nentries+ipnw );
          double px  =  HepMC::HEPEVT_Wrapper::px( nentries+ipnw );
          double py  =  HepMC::HEPEVT_Wrapper::py( nentries+ipnw );
          double pz  =  HepMC::HEPEVT_Wrapper::pz( nentries+ipnw );
          double e   =  HepMC::HEPEVT_Wrapper::e(  nentries+ipnw );
          double m   =  HepMC::HEPEVT_Wrapper::m(  nentries+ipnw );
       
          HepMC::GenParticle* NewPart = new HepMC::GenParticle( HepMC::FourVector(px,py,pz,e),
                                                                 pdg_id, status);
          NewPart->set_generated_mass( m );
          NewPart->suggest_barcode( nbcode );
          vtx->add_particle_out( NewPart ) ;
      }  
    
       //vtx->print();
       //std::cout << " leaving attachParticles() " << std::endl;
    
    } // end of global if-statement 
 
    return;
 }

void gen::PhotosInterface::avoidTauLeptonicDecays ( )

inline

Definition at line 33 of file PhotosInterface.h.

References fAvoidTauLeptonicDecays.

Referenced by gen::ExternalDecayDriver::ExternalDecayDriver().

33 { fAvoidTauLeptonicDecays=true; return; }

gen::PhotosInterface::fAvoidTauLeptonicDecays

bool fAvoidTauLeptonicDecays

Definition: PhotosInterface.h:47

void PhotosInterface::configureOnlyFor ( int ipdg )

Definition at line 58 of file PhotosInterface.cc.

References fOnlyPDG, and fSpecialSettings.

Referenced by gen::ExternalDecayDriver::ExternalDecayDriver().

 {
 
    fOnlyPDG = ipdg;
 //   std::ostringstream command;
 //   command << "QED-brem-off:" << fOnlyPDG ;
    fSpecialSettings.clear();
 //   fSpecialSettings.push_back( command.str() );
    
    return;
 
 }

void PhotosInterface::init ( void )

Definition at line 71 of file PhotosInterface.cc.

References fIsInitialized, and phoini_().

Referenced by GenMuonRadiationAlgorithm::compFSR(), and gen::ExternalDecayDriver::init().

 {
    
    if ( fIsInitialized ) return; // do init only once
    
    phoini_();
    
    fIsInitialized = true; 
 
    return;
 }

bool PhotosInterface::isTauLeptonicDecay ( HepMC::GenVertex * vtx )

Definition at line 451 of file PhotosInterface.cc.

References funct::abs().

Referenced by apply(), and applyToVertex().

 {
 
    if ( abs((*(vtx->particles_in_const_begin()))->pdg_id()) != 15 ) return false;
    
    for ( HepMC::GenVertex::particle_iterator pitr=vtx->particles_begin(HepMC::children);
          pitr != vtx->particles_end(HepMC::children); ++pitr) 
    {
       if ( abs((*pitr)->pdg_id()) == 11 || abs((*pitr)->pdg_id()) == 13 )
       {
          return true;
       }      
    } 
    
    return false;  
 
 }