AddFlowByShifting Class Reference

#include <AddFlowByShifting.h>

Inheritance diagram for AddFlowByShifting:

Collaboration diagram for AddFlowByShifting:

Public Member Functions
	AddFlowByShifting (const std::string &name, ISvcLocator *pSvcLocator)
StatusCode	initialize ()
StatusCode	execute (const EventContext &ctx)
	Execute method.
StatusCode	finalize ()
virtual StatusCode	sysInitialize () override
	Override sysInitialize.
virtual bool	isClonable () const override
	Specify if the algorithm is clonable.
virtual StatusCode	sysExecute (const EventContext &ctx) override
	Execute an algorithm.
virtual const DataObjIDColl &	extraOutputDeps () const override
	Return the list of extra output dependencies.
virtual bool	filterPassed (const EventContext &ctx) const
	Get filter decision:
virtual void	setFilterPassed (bool state, const EventContext &ctx) const
	Set filter decision:
ServiceHandle< StoreGateSvc > &	evtStore ()
	The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc.
const ServiceHandle< StoreGateSvc > &	detStore () const
	The standard StoreGateSvc/DetectorStore Returns (kind of) a pointer to the StoreGateSvc.
virtual StatusCode	sysStart () override
	Handle START transition.
virtual std::vector< Gaudi::DataHandle * >	inputHandles () const override
	Return this algorithm's input handles.
virtual std::vector< Gaudi::DataHandle * >	outputHandles () const override
	Return this algorithm's output handles.
Gaudi::Details::PropertyBase &	declareProperty (Gaudi::Property< T, V, H > &t)
void	updateVHKA (Gaudi::Details::PropertyBase &)
MsgStream &	msg () const
bool	msgLvl (const MSG::Level lvl) const

Static Public Member Functions
static double	vn_func (double x, void *params)

Protected Member Functions
virtual bool	isReEntrant () const override final
	Legacy algorithms are not thread-safe.
void	renounceArray (SG::VarHandleKeyArray &handlesArray)
	remove all handles from I/O resolution
std::enable_if_t< std::is_void_v< std::result_of_t< decltype(&T::renounce)(T)> > &&!std::is_base_of_v< SG::VarHandleKeyArray, T > &&std::is_base_of_v< Gaudi::DataHandle, T >, void >	renounce (T &h)
void	extraDeps_update_handler (Gaudi::Details::PropertyBase &ExtraDeps)
	Add StoreName to extra input/output deps as needed.

Private Types
enum	Harmonic {   NumHar =6 , v1 =0 , v2 =1 , v3 =2 ,   v4 =3 , v5 =4 , v6 =5 }
typedef ServiceHandle< StoreGateSvc >	StoreGateSvc_t

Private Member Functions
CLHEP::HepRandomEngine *	getRandomEngine (const std::string &streamName, const EventContext &ctx) const
double	SetParentToRanPhi (HepMC::GenParticlePtr parent, CLHEP::HepRandomEngine *rndmEngine)
double	AddFlowToParent (HepMC::GenParticlePtr parent, const HijingEventParams *hijing_pars)
void	MoveDescendantsToParent (HepMC::GenParticlePtr parent, double phishift)
void	jjia_minbias_new (double b, double eta, double pt)
void	jjia_minbias_new_v2only (double b, double eta, double pt)
void	fixed_vn (double b, double eta, double pt)
void	fixed_v2 (double b, double eta, double pt)
void	jjia_minbias_old (double b, double eta, double pt)
void	ao_test (double b, double eta, double pt)
void	custom_vn (double b, double eta, double pt)
void	p_Pb_cent_eta_indep (double b, double eta, double pt)
void	OO_eta_indep (double b, double eta, double pt)
void	Set_EbE_Fluctuation_Multipliers (HepMC::GenVertexPtr mainvtx, float b, CLHEP::HepRandomEngine *rndmEngine)
Gaudi::Details::PropertyBase &	declareGaudiProperty (Gaudi::Property< T, V, H > &hndl, const SG::VarHandleKeyType &)
	specialization for handling Gaudi::Property<SG::VarHandleKey>

Private Attributes
void(AddFlowByShifting::*	m_flow_function )(double b, double eta, double pt)
TGraph *	m_graph_fluc {}
ServiceHandle< IAthRNGSvc >	m_rndmSvc {this, "RndmSvc", "AthRNGSvc"}
StringProperty	m_inkey {this, "McTruthKey" , "GEN_EVENT" }
StringProperty	m_outkey {this, "McFlowKey" , "FLOW_EVENT" }
IntegerProperty	m_ranphi_sw {this, "RandomizePhi" , 0 }
StringProperty	m_flow_function_name {this, "FlowFunctionName" , "jjia_minbias_new"}
StringProperty	m_flow_implementation {this, "FlowImplementation", "exact" }
int	m_flow_implementation_type {0}
BooleanProperty	m_flow_fluctuations {this, "FlowFluctuations" , false }
IntegerProperty	m_floweta_sw {this, "FlowEtaSwitch" , 0 }
FloatProperty	m_flow_maxeta {this, "FlowMaxEtaCut" , 10.0 }
FloatProperty	m_flow_mineta {this, "FlowMinEtaCut" , 0.f }
IntegerProperty	m_flowpt_sw {this, "FlowPtSwitch" , 0 }
FloatProperty	m_flow_maxpt {this, "FlowMaxPtCut" , 1000000.f }
FloatProperty	m_flow_minpt {this, "FlowMinPtCut" , 0.f }
IntegerProperty	m_flowb_sw {this, "FlowBSwitch" , 0 }
FloatProperty	m_custom_v1 {this, "custom_v1" , 0.f }
FloatProperty	m_custom_v2 {this, "custom_v2" , 0.f }
FloatProperty	m_custom_v3 {this, "custom_v3" , 0.f }
FloatProperty	m_custom_v4 {this, "custom_v4" , 0.f }
FloatProperty	m_custom_v5 {this, "custom_v5" , 0.f }
FloatProperty	m_custom_v6 {this, "custom_v6" , 0.f }
int	m_particles_processed {0}
float	m_psi_n [Harmonic::NumHar]
float	m_v_n [Harmonic::NumHar]
float	m_EbE_Multiplier_vn [Harmonic::NumHar]
DataObjIDColl	m_extendedExtraObjects
	Extra output dependency collection, extended by AthAlgorithmDHUpdate to add symlinks.
StoreGateSvc_t	m_evtStore
	Pointer to StoreGate (event store by default).
StoreGateSvc_t	m_detStore
	Pointer to StoreGate (detector store by default).
std::vector< SG::VarHandleKeyArray * >	m_vhka
bool	m_varHandleArraysDeclared

Detailed Description

Definition at line 39 of file AddFlowByShifting.h.

Member Typedef Documentation

StoreGateSvc_t

typedef ServiceHandle<StoreGateSvc> AthCommonDataStore< AthCommonMsg< Gaudi::Algorithm > >::StoreGateSvc_t

privateinherited

Definition at line 388 of file AthCommonDataStore.h.

Member Enumeration Documentation

Harmonic

enum AddFlowByShifting::Harmonic

private

Enumerator
NumHar
v1
v2
v3
v4
v5
v6

Definition at line 101 of file AddFlowByShifting.h.

               {
    NumHar=6,
    v1=0,
    v2=1,
    v3=2,
    v4=3,
    v5=4,
    v6=5,
  };

Constructor & Destructor Documentation

AddFlowByShifting()

AddFlowByShifting::AddFlowByShifting	(	const std::string &	name,
		ISvcLocator *	pSvcLocator )

Definition at line 56 of file AddFlowByShifting.cxx.

                                                                                    :
  AthAlgorithm(name, pSvcLocator)
{
  m_flow_function= NULL;
  for(int ihar = 0; ihar< Harmonic::NumHar; ihar++){
    m_psi_n[ihar] =0.0;
    m_v_n  [ihar] =0.0;
    m_EbE_Multiplier_vn[ihar]=1.0;
  }
}

Member Function Documentation

AddFlowToParent()

double AddFlowByShifting::AddFlowToParent ( HepMC::GenParticlePtr parent, const HijingEventParams * hijing_pars )

private

Definition at line 458 of file AddFlowByShifting.cxx.

{
  m_particles_processed++;
  CLHEP::HepLorentzVector momentum(parent->momentum().px(),
                                   parent->momentum().py(),
                                   parent->momentum().pz(),
                                   parent->momentum().e());
  double pt    = parent->momentum().perp();
  double eta   = parent->momentum().pseudoRapidity();
  double phi_0 = parent->momentum().phi();
 
  int error_=0;
  if(pt    !=pt)    {ATH_MSG_ERROR("ERROR pt  of track  is not defined");error_=1;} //true if pt==nan
  if(eta   !=eta)   {ATH_MSG_ERROR("ERROR eta of track  is not defined");error_=1;}
  if(phi_0 !=phi_0) {ATH_MSG_ERROR("ERROR phi of track  is not defined");error_=1;}
  if(error_==1){
    ATH_MSG_ERROR("Original Particle Momentum(px,py,pz,e,m)="<<parent->momentum().px()<<"  "
                  <<parent->momentum().py()<<"  "
                  <<parent->momentum().pz()<<"  "
                  <<parent->momentum().e() <<"  "
                  <<parent->momentum().m() <<"  ");
  }
 
 
 
 
  //Call the appropriate function to set the vn values
  for(int ihar = 0; ihar< Harmonic::NumHar; ihar++){m_v_n  [ihar]=0.0;} //reset the vn for this particle
  double b = hijing_pars->get_b();
  (*this.*m_flow_function)(b,eta,pt);//Set the vn for this particle
 
  //add EbE fluctuations
  if(m_flow_fluctuations){
    for(int ihar = 0; ihar< Harmonic::NumHar; ihar++){
      m_v_n[ihar] *= m_EbE_Multiplier_vn[ihar];
      if(m_v_n[ihar]>0.5){
        ATH_MSG_WARNING(" Vn Too large  "<<ihar+1<<"  "<<m_EbE_Multiplier_vn[ihar]<<"  "<<m_v_n[ihar]);m_v_n[ihar]=0.5;
      }
    }
  }
 
  double phishift=0;
 
  // Old fashioned rotation(approximate)
  if (m_flow_implementation_type==0){
    float phi=phi_0;
    phishift=  -2*( m_v_n[Harmonic::v1]*sin(1*(phi-m_psi_n[Harmonic::v1]))/1.0 + 
                    m_v_n[Harmonic::v2]*sin(2*(phi-m_psi_n[Harmonic::v2]))/2.0 +
                    m_v_n[Harmonic::v3]*sin(3*(phi-m_psi_n[Harmonic::v3]))/3.0 + 
                    m_v_n[Harmonic::v4]*sin(4*(phi-m_psi_n[Harmonic::v4]))/4.0 +
                    m_v_n[Harmonic::v5]*sin(5*(phi-m_psi_n[Harmonic::v5]))/5.0 + 
                    m_v_n[Harmonic::v6]*sin(6*(phi-m_psi_n[Harmonic::v6]))/6.0 );
 
  }
 
  // New fashioned rotation(exact)
  else if (m_flow_implementation_type==1){
    // Thread-safe according to https://www.gnu.org/software/gsl/doc/html/roots.html
    const gsl_root_fsolver_type *T ATLAS_THREAD_SAFE = gsl_root_fsolver_brent;
    gsl_root_fsolver *s = gsl_root_fsolver_alloc (T);
    double x_lo=-2*M_PI,x_hi=2*M_PI;
    float params[13];
    for(int ipar=0;ipar<13;ipar++) {params[ipar]=0;}
    gsl_function F;
    F.function = &(AddFlowByShifting::vn_func);
    F.params =&params;
    gsl_root_fsolver_set (s, &F, x_lo, x_hi);
    int iter=0;
    params[ 0]=phi_0;
    params[ 1]=m_v_n  [Harmonic::v1];
    params[ 2]=m_v_n  [Harmonic::v2];
    params[ 3]=m_v_n  [Harmonic::v3];
    params[ 4]=m_v_n  [Harmonic::v4];
    params[ 5]=m_v_n  [Harmonic::v5];
    params[ 6]=m_v_n  [Harmonic::v6];
    params[ 7]=m_psi_n[Harmonic::v1];
    params[ 8]=m_psi_n[Harmonic::v2];
    params[ 9]=m_psi_n[Harmonic::v3];
    params[10]=m_psi_n[Harmonic::v4];
    params[11]=m_psi_n[Harmonic::v5];
    params[12]=m_psi_n[Harmonic::v6];
    int status;
    double phi=phi_0;
    do
      {
        iter++;
        status = gsl_root_fsolver_iterate (s);
        phi = gsl_root_fsolver_root (s);
        x_lo = gsl_root_fsolver_x_lower (s);
        x_hi = gsl_root_fsolver_x_upper (s);
        status = gsl_root_test_interval (x_lo, x_hi,0, 0.00001);
      }
    while (status == GSL_CONTINUE && iter < 1000);
    gsl_root_fsolver_free (s);
 
    if (iter>=1000) return 0;
 
    phishift = phi-phi_0;
  }
 
  if(std::abs(phishift) > 1e-7) {
    momentum.rotateZ(phishift*Gaudi::Units::rad);
    parent->set_momentum( HepMC::FourVector(momentum.px(),momentum.py(),momentum.pz(),momentum.e()) );
  }
  ATH_MSG_DEBUG( "Parent particle:"   <<
                 " V1 = "             << m_v_n[Harmonic::v1] <<
                 " V2 = "             << m_v_n[Harmonic::v2] <<
                 " V3 = "             << m_v_n[Harmonic::v3] <<
                 " V4 = "             << m_v_n[Harmonic::v4] <<
                 " V5 = "             << m_v_n[Harmonic::v5] <<
                 " V6 = "             << m_v_n[Harmonic::v6] <<
                 " Phi shift = "      << phishift            <<
                 " Phi shifted = "    << momentum.phi() );
 
  return phishift;
}

ao_test()

void AddFlowByShifting::ao_test ( double b, double eta, double pt )

private

Definition at line 655 of file AddFlowByShifting.cxx.

{
  for(int ihar=0;ihar<Harmonic::NumHar;ihar++){
    m_v_n[ihar]=0.0;
  }
 
  pt/=1000;
  if(pt>2) pt = 2; // flat max at pt > 2
 
  m_v_n[Harmonic::v2] = 0.02 * b * pt;
}

custom_vn()

void AddFlowByShifting::custom_vn ( double b, double eta, double pt )

private

Definition at line 668 of file AddFlowByShifting.cxx.

{
  m_v_n[Harmonic::v1]=m_custom_v1;
  m_v_n[Harmonic::v2]=m_custom_v2;
  m_v_n[Harmonic::v3]=m_custom_v3;
  m_v_n[Harmonic::v4]=m_custom_v4;
  m_v_n[Harmonic::v5]=m_custom_v5;
  m_v_n[Harmonic::v6]=m_custom_v6;
}

declareGaudiProperty()

Gaudi::Details::PropertyBase & AthCommonDataStore< AthCommonMsg< Gaudi::Algorithm > >::declareGaudiProperty ( Gaudi::Property< T, V, H > & hndl, const SG::VarHandleKeyType &  )

inlineprivateinherited

specialization for handling Gaudi::Property<SG::VarHandleKey>

Definition at line 156 of file AthCommonDataStore.h.

  {
    return *AthCommonDataStore<PBASE>::declareProperty(hndl.name(),
                                                       hndl.value(), 
                                                       hndl.documentation());
 
  }

declareProperty()

Gaudi::Details::PropertyBase & AthCommonDataStore< AthCommonMsg< Gaudi::Algorithm > >::declareProperty ( Gaudi::Property< T, V, H > & t )

inlineinherited

Definition at line 145 of file AthCommonDataStore.h.

                                                                     {
    typedef typename SG::HandleClassifier<T>::type htype;
    return AthCommonDataStore<PBASE>::declareGaudiProperty(t, htype());
  }

detStore()

const ServiceHandle< StoreGateSvc > & AthCommonDataStore< AthCommonMsg< Gaudi::Algorithm > >::detStore ( ) const

inlineinherited

The standard StoreGateSvc/DetectorStore Returns (kind of) a pointer to the StoreGateSvc.

Definition at line 95 of file AthCommonDataStore.h.

{ return m_detStore; }

evtStore()

ServiceHandle< StoreGateSvc > & AthCommonDataStore< AthCommonMsg< Gaudi::Algorithm > >::evtStore ( )

inlineinherited

The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc.

Definition at line 85 of file AthCommonDataStore.h.

{ return m_evtStore; }

execute()

StatusCode AddFlowByShifting::execute ( const EventContext & ctx )

virtual

Execute method.

Implements AthAlgorithm.

Definition at line 231 of file AddFlowByShifting.cxx.

                                                             {
  ATH_MSG_INFO(">>> AddFlowByShifting from execute");
 
  CLHEP::HepRandomEngine *rndmEngine = getRandomEngine("FLOW", ctx);
  // Get hijing event parameters
  const HijingEventParams *hijing_pars;
  if( evtStore()->retrieve(hijing_pars, "Hijing_event_params").isFailure() ) {
    ATH_MSG_ERROR("Could not retrieve Hijing_event_params");
    return StatusCode::FAILURE;
  }
  ATH_MSG_INFO("Event parameters: B = " << hijing_pars->get_b()<<
               "  BPhi = " << hijing_pars->get_bphi());
 
 
  // FIXME: changing data in the event store
  HijingEventParams *hijing_pars_nc = const_cast<HijingEventParams*> (hijing_pars);
 
 
  // Read Data from Transient Store
  const McEventCollection* mcCollptr;
  if ( evtStore()->retrieve(mcCollptr, m_inkey).isFailure() ) {
    ATH_MSG_ERROR("Could not retrieve truth McEventCollection");
    return StatusCode::FAILURE;
  }
 
 
  // Loop over all events in original McEventCollection and
  // Copy to a new (modifiable) collection
  McEventCollection::const_iterator citr;
  McEventCollection* mcFlowCollptr = new  McEventCollection();
  for (citr = mcCollptr->begin(); citr!=mcCollptr->end(); ++citr) {
    mcFlowCollptr->push_back(new HepMC::GenEvent(*(*citr)));
  }
 
 
  //Geneate the event-plane angles (some of them may or may not be used later on)
  //Store the angles into the hijing event parameters
  for(int ihar=0;ihar<6;ihar++){
    m_psi_n[ihar] =(CLHEP::RandFlat::shoot(rndmEngine)-0.5)*2*M_PI / (ihar+1);   //Principal value must be within -PI/n to PI/n
    hijing_pars_nc->set_psi(ihar+1,m_psi_n[ihar]);
  }
  m_psi_n[1]=hijing_pars->get_bphi()                   ;//the psi2 plane is aligned with the impact parameter
  m_psi_n[1]=std::atan2(std::sin(2*m_psi_n[1]),std::cos(2*m_psi_n[1]))/2.0;//ensure that Psi2 is within [-PI/2,PI/2]
  hijing_pars_nc->set_psi(2,m_psi_n[1]);
  ATH_MSG_DEBUG(" Psi2 for event : "<<(*hijing_pars).get_psi(2));
 
 
  // Add flow by phi angle shifting
  McEventCollection::iterator itr;
  for (itr = mcFlowCollptr->begin(); itr!=mcFlowCollptr->end(); ++itr) {
    ATH_MSG_DEBUG("Next event in the bag ...");
 
 
 
    auto mainvtx=(*itr)->vertices().front();
    if(m_flow_fluctuations) Set_EbE_Fluctuation_Multipliers(mainvtx,hijing_pars->get_b(),rndmEngine);
    int particles_in_event = (*itr)->particles().size();
    m_particles_processed = 0;
    for ( auto parent: mainvtx->particles_out())
      {
        // Process particles from main vertex
        CLHEP::HepLorentzVector momentum(parent->momentum().px(),
                                         parent->momentum().py(),
                                         parent->momentum().pz(),
                                         parent->momentum().e());
        ATH_MSG_DEBUG("Parent particle: " << parent        <<
                      " Eta = "           << momentum.pseudoRapidity()<<
                      " Phi = "           << momentum.phi()            );
 
        //skip particle if eta is outside implementation range
        if(m_floweta_sw){
          float eta=std::abs(momentum.pseudoRapidity());
          if (eta<m_flow_mineta || eta> m_flow_maxeta) continue;
        }
 
        //skip particle if pT is outside implementation range
        if(m_flowpt_sw){
          float pT=momentum.perp();
          if (pT<m_flow_minpt || pT> m_flow_maxpt) continue;
        }
 
        // Randomize phi if explicitely requested
        if(m_ranphi_sw) {
          double phishift = SetParentToRanPhi(parent, rndmEngine);
          MoveDescendantsToParent(parent, phishift)  ;// adjust descendants to parent position
        }
 
        // Add flow to particles from main vertex
        double phishift = AddFlowToParent(parent, hijing_pars);
        MoveDescendantsToParent(std::move(parent), phishift);// adjust descendants to parent position
      }
 
    // correct for double counting
    if(m_ranphi_sw) m_particles_processed /= 2;
    // correct for incoming particles
    ATH_MSG_INFO( " Particles in event: " << particles_in_event <<
                  " Processed for flow: " << m_particles_processed+2);
    if(particles_in_event != (m_particles_processed+2)){
      ATH_MSG_WARNING( " Particles in event: " << particles_in_event <<
                       " Processed for flow: " << m_particles_processed+2);
    }
  }
 
  if(evtStore()->record(mcFlowCollptr, m_outkey).isFailure()){
    ATH_MSG_ERROR("Could not record flow McEventCollection");
    return StatusCode::FAILURE;
  }
  return StatusCode::SUCCESS;
}

extraDeps_update_handler()

void AthCommonDataStore< AthCommonMsg< Gaudi::Algorithm > >::extraDeps_update_handler ( Gaudi::Details::PropertyBase & ExtraDeps )

protectedinherited

Add StoreName to extra input/output deps as needed.

use the logic of the VarHandleKey to parse the DataObjID keys supplied via the ExtraInputs and ExtraOuputs Properties to add the StoreName if it's not explicitly given

extraOutputDeps()

const DataObjIDColl & AthCommonAlgorithm< Gaudi::Algorithm >::extraOutputDeps ( ) const

overridevirtualinherited

Return the list of extra output dependencies.

This list is extended to include symlinks implied by inheritance relations.

Definition at line 89 of file AthCommonAlgorithm.cxx.

{
  // If we didn't find any symlinks to add, just return the collection
  // from the base class.  Otherwise, return the extended collection.
  if (!m_extendedExtraObjects.empty()) {
    return m_extendedExtraObjects;
  }
  return BaseAlg::extraOutputDeps();
}

filterPassed()

virtual bool AthCommonAlgorithm< Gaudi::Algorithm >::filterPassed ( const EventContext & ctx ) const

inlinevirtualinherited

Get filter decision:

Definition at line 93 of file AthCommonAlgorithm.h.

                                                           {
    return execState( ctx ).filterPassed();
  }

finalize()

StatusCode AddFlowByShifting::finalize

(

)

Definition at line 342 of file AddFlowByShifting.cxx.

                                       {
  ATH_MSG_INFO(">>> AddFlowByShifting from finalize <<<");
  if(m_graph_fluc){
    delete m_graph_fluc ;
    m_graph_fluc=nullptr;
  }
 
  // End of finalization step
  return StatusCode::SUCCESS;
}

fixed_v2()

void AddFlowByShifting::fixed_v2 ( double b, double eta, double pt )

private

Definition at line 632 of file AddFlowByShifting.cxx.

{
  for(int ihar=0;ihar<Harmonic::NumHar;ihar++){
    m_v_n[ihar]=0.0;
  }
  m_v_n[Harmonic::v2]=0.0500;
}

fixed_vn()

void AddFlowByShifting::fixed_vn ( double b, double eta, double pt )

private

Definition at line 620 of file AddFlowByShifting.cxx.

{
  m_v_n[Harmonic::v1]=0.0000;
  m_v_n[Harmonic::v2]=0.0500;
  m_v_n[Harmonic::v3]=0.0280;
  m_v_n[Harmonic::v4]=0.0130;
  m_v_n[Harmonic::v5]=0.0045;
  m_v_n[Harmonic::v6]=0.0015;
}

getRandomEngine()

CLHEP::HepRandomEngine * AddFlowByShifting::getRandomEngine ( const std::string & streamName, const EventContext & ctx ) const

private

Definition at line 221 of file AddFlowByShifting.cxx.

{
  ATHRNG::RNGWrapper* rngWrapper = m_rndmSvc->getEngine(this, streamName);
  std::string rngName = name()+streamName;
  rngWrapper->setSeed( rngName, ctx );
  return rngWrapper->getEngine(ctx);
}

initialize()

StatusCode AddFlowByShifting::initialize

(

)

Definition at line 68 of file AddFlowByShifting.cxx.

                                        {
  ATH_MSG_INFO(">>> AddFlowByShifting from Initialize <<<");
 
  ATH_CHECK(m_rndmSvc.retrieve());
 
  ATH_MSG_INFO("**********Settings for Afterburner************");
  ATH_MSG_INFO("McTruthKey         : " << m_inkey              );
  ATH_MSG_INFO("McFlowKey          : " << m_outkey             );
 
  ATH_MSG_INFO("FlowFunctionName   : " << m_flow_function_name );
  ATH_MSG_INFO("FlowInplementation : " << m_flow_implementation);
  ATH_MSG_INFO("FlowFluctuations   : " << m_flow_fluctuations  );
 
  ATH_MSG_INFO("RandomizePhi       : " << m_ranphi_sw          );
 
  ATH_MSG_INFO("FlowEtaSwitch      : " << m_floweta_sw         );
  ATH_MSG_INFO("FlowMinEtaCut      : " << m_flow_mineta        );
  ATH_MSG_INFO("FlowMaxEtaCut      : " << m_flow_maxeta        );
 
  ATH_MSG_INFO("FlowPtSwitch       : " << m_flowpt_sw          );
  ATH_MSG_INFO("FlowMinPtCut       : " << m_flow_minpt         );
  ATH_MSG_INFO("FlowMaxPtCut       : " << m_flow_maxpt         );
 
  ATH_MSG_INFO("FlowV1             : " << m_custom_v1          );
  ATH_MSG_INFO("FlowV2             : " << m_custom_v2          );
  ATH_MSG_INFO("FlowV3             : " << m_custom_v3          );
  ATH_MSG_INFO("FlowV4             : " << m_custom_v4          );
  ATH_MSG_INFO("FlowV5             : " << m_custom_v5          );
  ATH_MSG_INFO("FlowV6             : " << m_custom_v6          );
  ATH_MSG_INFO("FlowBSwitch        : " << m_flowb_sw           );
  ATH_MSG_INFO("********************************r*************");
 
 
  // Select the flow-implementing function based of the function-choice variable
  m_flow_function=&AddFlowByShifting::jjia_minbias_new;
  if     (m_flow_function_name=="jjia_minbias_new"       ) m_flow_function=&AddFlowByShifting::jjia_minbias_new;
  else if(m_flow_function_name=="jjia_minbias_new_v2only") m_flow_function=&AddFlowByShifting::jjia_minbias_new_v2only;
  else if(m_flow_function_name=="fixed_vn"               ) m_flow_function=&AddFlowByShifting::fixed_vn;
  else if(m_flow_function_name=="fixed_v2"               ) m_flow_function=&AddFlowByShifting::fixed_v2;
  else if(m_flow_function_name=="jjia_minbias_old"       ) m_flow_function=&AddFlowByShifting::jjia_minbias_old;
  else if(m_flow_function_name=="ao_test"                ) m_flow_function=&AddFlowByShifting::ao_test;
  else if(m_flow_function_name=="custom"                 ) m_flow_function=&AddFlowByShifting::custom_vn;
  else if(m_flow_function_name=="p_Pb_cent_eta_indep"    ) m_flow_function=&AddFlowByShifting::p_Pb_cent_eta_indep;
  else if(m_flow_function_name=="OO_eta_indep"           ) m_flow_function=&AddFlowByShifting::OO_eta_indep;
  else{
    ATH_MSG_ERROR("Unimplemented option for setting 'FlowFunctionName' : " << m_flow_function_name);
    return StatusCode::FAILURE;
  }
 
 
  m_flow_implementation_type=1;
  if     (m_flow_implementation=="approximate"){
    m_flow_implementation_type=0;
    ATH_MSG_WARNING("'FlowInplementation=\"approximate\"' is obsolete, please switch to 'FlowInplementation=\"exact\"' ");
  }
  else if(m_flow_implementation=="exact"      ){
    m_flow_implementation_type=1;
  }
  else{
    ATH_MSG_ERROR("Unimplemented option for  setting 'FlowInplementation' : " << m_flow_implementation);
    return StatusCode::FAILURE;
  }
 
 
 
  //TGraph storing the v2_RP/delta Vs b_imp values to be used in implementing the EbyE fluctuations
  //the values below are b_imp-low,b_imp-high, delta/v2_RP for different centralities
  //underflow and overflow bins are added for smooth extrapolation
  if(m_flow_fluctuations){
    //Fluctuations for Pb+Pb
    if(m_flow_function_name=="jjia_minbias_new"        ||
       m_flow_function_name=="jjia_minbias_new_v2only")
    {
      //The delta/v2_RP values are taken from Fig15 of EbE vn paper (arXiv:1305.2942)
      //                 <0  ,  0-1 ,  1-2 , 2-3  ,  3-4 ,  4-5 , 5-10 , 10-15, 15-20, 20-25, 25-30,
      float b_lo[21]={  -1.00, 0.000, 1.483, 2.098, 2.569, 2.966, 3.317, 4.687, 5.739, 6.627, 7.409,
                        8.117, 8.767, 9.373, 9.943,10.479,10.991,11.479,11.947,15.00 ,100.0};
      //                30-35, 35-40, 40-45, 45-50, 50-55, 55-60, 60-65, 65-70, 70-
      float b_hi[21]={ -1.00 , 1.483, 2.098, 2.569, 2.966, 3.317, 4.687, 5.739, 6.627, 7.409, 8.117,//bimp_high
                        8.767, 9.373, 9.943,10.479,10.991,11.479,11.947,12.399,15.00 ,100.0};
      float val [21]={  5.600, 5.600, 5.600,1.175 ,0.8253,0.7209,0.5324,0.4431,0.3984,0.3844,0.3847,
                       0.3935,0.4106,0.4310,0.4574,0.4674,0.4873,0.4796,0.4856,0.5130,0.5130};
      float bimp_vals[21];
 
      for(int i=0;i<21;i++){
        bimp_vals[i]=(b_lo[i]+b_hi[i])/2.0f;
        val      [i]=1.0f/val[i];//change to v2_RP/delta
        val      [i]=1.0/val[i];//change to v2_RP/delta
      }
 
      m_graph_fluc=new TGraph(21,bimp_vals,val);
    }
 
    //Fluctuations for O+O
    if(m_flow_function_name=="OO_eta_indep")
    {
      //The v2{2} and v2{4} values are taken from the OO paper (arXiv:2509.05171)
      //The v2{4} values beyond the 50-60% bin are put in by hand
      //The v2{2} values beyond the 75-80% bin are put in by hand
      //The impact parameters are obtained from HIJING
 
      //4-particle cumulant v2 (v2{4})
      //                       <0,       0-1,       1-2,       2-3,       3-4,       4-5,      5-10,     10-15,     15-20,     20-25,     25-30,
      //                    30-35,     35-40,     40-50,     50-60,     60-80,    80-100,      >100,
      float v2_4 [18]={ 0.0444322, 0.0444322, 0.0476009, 0.0495882, 0.0513053, 0.0527020, 0.0563533, 0.0614435, 0.0650522, 0.0675199, 0.0693970,
                        0.0705456, 0.0706748, 0.0686654, 0.0571383, 0.04     , 0.02     , 0.0      };//last three bins are put in by hand
      float bimp4[18]={-1        ,  0.481843, 0.886609 , 1.13462  , 1.33981  , 1.52064  , 1.95342  , 2.51881  , 2.97665  , 3.37416  , 3.72876  ,
                        4.05377  ,  4.35564 , 4.77463  , 5.30719  , 6.08296  , 7.6415   , 15       };
 
 
      // template-fit v2  (v2{2})
      //                       <0,       0-1,       1-2,       2-3,       3-4,       4-5,      5-10,     10-15,     15-20,     20-25,     25-30,
      //                    30-35,     35-40,     40-45,     45-50,     50-55,     55-60,     60-65,     65-70,     70-75,     75-80,    80-100,>100
      float v2_2 [23]={ 0.0749381, 0.0749381, 0.0782992, 0.0801742, 0.0814646, 0.0826611, 0.0853616, 0.0889407, 0.0913229, 0.0928072, 0.0933774,
                        0.0933299, 0.0926254, 0.0910148, 0.0889322, 0.086161 , 0.0828212, 0.0793132, 0.0755036, 0.0714726, 0.0668264, 0.055    , 0.05};//last two bins are put in by hand
      float bimp2[23]={-1        ,  0.481843, 0.886609 , 1.13462  , 1.33981  , 1.52064  , 1.95342  , 2.51881  , 2.97665  , 3.37416  , 3.72876  ,
                        4.05377  ,  4.35564 , 4.63868  , 4.91057  , 5.17773  , 5.43665  , 5.68739  , 5.94258  , 6.20668  , 6.49519  , 7.6415   , 15};
 
 
      //Evaluate ratio of v2{4}/delta 
      TGraph *graph_v2_4=new TGraph(18,bimp4,v2_4);
      TGraph *graph_v2_2=new TGraph(23,bimp2,v2_2);
      float ratio[23]={0};
      for(int i=0;i<23;i++){
        double vn_RP=graph_v2_4->Eval(bimp2[i]);
        double vn_2 =graph_v2_2->Eval(bimp2[i]);
        double delta=vn_2*vn_2 - vn_RP*vn_RP;
        if(delta>0){
          ratio[i]= vn_RP/sqrt(delta/2.0);
        }
        else{
          ATH_MSG_ERROR("vn{2}<vn{4} for b_imp = " << bimp2[i]);
          return StatusCode::FAILURE;
        }
      }
 
      m_graph_fluc=new TGraph(23,bimp2,ratio);
      delete graph_v2_4;
      delete graph_v2_2;
    }
    //Fluctuations for Ne+Ne (TODO)
 
    else{
       ATH_MSG_ERROR("Flow fluctuations are not implemented for the following case: " << m_flow_function_name);
    }
  }
 
 
  // Initialization terminated
  return StatusCode::SUCCESS;
}

inputHandles()

virtual std::vector< Gaudi::DataHandle * > AthCommonDataStore< AthCommonMsg< Gaudi::Algorithm > >::inputHandles ( ) const

overridevirtualinherited

Return this algorithm's input handles.

We override this to include handle instances from key arrays if they have not yet been declared. See comments on updateVHKA.

isClonable()

virtual bool AthCommonAlgorithm< Gaudi::Algorithm >::isClonable ( ) const

inlineoverridevirtualinherited

Specify if the algorithm is clonable.

Only relevant for non-reentrant algorithms. Actual number of clones needs to be set via the "Cardinality" property.

Reimplemented in AFP_DigiTop, AlgB, AlgT, BCM_Digitization, CscDigitBuilder, CscDigitToCscRDO, G4AtlasAlg, G4RunAlg, HGTD_Digitization, HiveAlgBase, InDet::GNNSeedingTrackMaker, InDet::SCT_Clusterization, InDet::SiSPGNNTrackMaker, InDet::SiSPSeededTrackFinder, InDet::SiTrackerSpacePointFinder, ISF::SimKernelMT, ITk::StripDigitization, ITkPixelCablingAlg, ITkStripCablingAlg, LArHitEMapMaker, LArTTL1Maker, LUCID_DigiTop, LVL1::L1TopoSimulation, MergeCalibHits, MergeGenericMuonSimHitColl, MergeHijingPars, MergeMcEventCollection, MergeTrackRecordCollection, MergeTruthJets, MergeTruthParticles, MuonDigitizer, PileUpMTAlg, PixelDigitization, RoIBResultToxAOD, SCT_ByteStreamErrorsTestAlg, SCT_CablingCondAlgFromCoraCool, SCT_CablingCondAlgFromText, SCT_ConditionsParameterTestAlg, SCT_ConditionsSummaryTestAlg, SCT_ConfigurationConditionsTestAlg, SCT_Digitization, SCT_FlaggedConditionTestAlg, SCT_LinkMaskingTestAlg, SCT_MajorityConditionsTestAlg, SCT_ModuleVetoTestAlg, SCT_MonitorConditionsTestAlg, SCT_PrepDataToxAOD, SCT_RawDataToxAOD, SCT_ReadCalibChipDataTestAlg, SCT_ReadCalibDataTestAlg, SCT_RODVetoTestAlg, SCT_SensorsTestAlg, SCT_SiliconConditionsTestAlg, SCT_StripVetoTestAlg, SCT_TdaqEnabledTestAlg, SCT_TestCablingAlg, SCTEventFlagWriter, SCTRawDataProvider, SCTSiLorentzAngleTestAlg, SCTSiPropertiesTestAlg, SGInputLoader, Simulation::BeamEffectsAlg, TileHitVecToCnt, TileMuonFitter, TilePulseForTileMuonReceiver, TileRawChannelMaker, TRTDigitization, and ZDC_DigiTop.

Definition at line 68 of file AthCommonAlgorithm.h.

                                           {
    return true;
  }

isReEntrant()

virtual bool AthAlgorithm::isReEntrant ( ) const

inlinefinaloverrideprotectedvirtualinherited

Legacy algorithms are not thread-safe.

Definition at line 47 of file AthAlgorithm.h.

{ return false; }

jjia_minbias_new()

void AddFlowByShifting::jjia_minbias_new ( double b, double eta, double pt )

private

Definition at line 578 of file AddFlowByShifting.cxx.

{
  pt=pt/1000.0; //convert to GeV
 
  float a1,a2,a3,a4;
  a1=0.4397*std::exp(-(b-4.526)*(b-4.526)/72.0) + 0.636;
  a2=1.916/(b+2) +0.1;
  a3=4.79*0.0001*(b-0.621)*(b-10.172)*(b-23)+1.2;   // this is >0 for b>0
  a4=0.135*std::exp(-0.5*(b-10.855)*(b-10.855)/4.607/4.607) +0.0120;
 
  float temp1 = std::pow(pt    , a1) / (1+std::exp( (pt-3.0)/a3));
  float temp2 = std::pow(pt+0.1,-a2) / (1+std::exp(-(pt-4.5)/a3));
  float temp3 =  0.01           / (1+std::exp(-(pt-4.5)/a3));
 
  m_v_n[Harmonic::v2] = ( a4*(temp1+temp2) + temp3 )* std::exp(-0.5* eta*eta /6.27/6.27) ;
 
  float fb=0.97 +1.06*std::exp(-0.5*b*b/3.2/3.2);
  m_v_n[Harmonic::v3]=std::pow(fb*std::sqrt(m_v_n[1]),3);
 
  float gb= 1.096 +1.36 *std::exp(-0.5*b*b/3.0/3.0);
  gb=gb*sqrt(m_v_n[1]);
  m_v_n[Harmonic::v4]=pow(gb,4);
  m_v_n[Harmonic::v5]=pow(gb,5);
  m_v_n[Harmonic::v6]=pow(gb,6);
  m_v_n[Harmonic::v1]=0;
}

jjia_minbias_new_v2only()

void AddFlowByShifting::jjia_minbias_new_v2only ( double b, double eta, double pt )

private

Definition at line 607 of file AddFlowByShifting.cxx.

{
  jjia_minbias_new(b, eta, pt);
  //Set all harmonics except v2 to 0
  m_v_n[Harmonic::v1] = 0;
  m_v_n[Harmonic::v3] = 0;
  m_v_n[Harmonic::v4] = 0;
  m_v_n[Harmonic::v5] = 0;
  m_v_n[Harmonic::v6] = 0;
}

jjia_minbias_old()

void AddFlowByShifting::jjia_minbias_old ( double b, double eta, double pt )

private

Definition at line 642 of file AddFlowByShifting.cxx.

{
  m_v_n[Harmonic::v1] = 0;
  m_v_n[Harmonic::v2] = 0.03968 * b
    * (1 - 2.1/(1 + std::exp(1.357*(pt/1000))))
    * std::exp(-(eta*eta)/(2*6.37*6.37));
  m_v_n[Harmonic::v3]=0.0000;
  m_v_n[Harmonic::v4]=0.0000;
  m_v_n[Harmonic::v5]=0.0000;
  m_v_n[Harmonic::v6]=0.0000;
}

MoveDescendantsToParent()

void AddFlowByShifting::MoveDescendantsToParent ( HepMC::GenParticlePtr parent, double phishift )

private

Definition at line 381 of file AddFlowByShifting.cxx.

{
  // Move the branch of descendant vertices and particles
  // by phishift to parent particle position
  auto endvtx = parent->end_vertex();
  if ( endvtx ) {
    ATH_MSG_DEBUG("Processing branch of parent particle "<< parent);
 
    //Added October 2025
    // --- Rotate the parent’s end vertex itself ---
    if (std::abs(phishift) > 1e-7) {
      CLHEP::HepLorentzVector pos(endvtx->position().x(),
                                  endvtx->position().y(),
                                  endvtx->position().z(),
                                  endvtx->position().t());
      pos.rotateZ(phishift * Gaudi::Units::rad);
      endvtx->set_position(HepMC::FourVector(pos.x(), pos.y(), pos.z(), pos.t()));
    }
 
 
    //Added October 2025
    // --- Rotate the parent’s immediate daughters (outgoing of endvtx) ---
    for (auto child : endvtx->particles_out()) 
    {
      CLHEP::HepLorentzVector p(child->momentum().px(),
                                child->momentum().py(),
                                child->momentum().pz(),
                                child->momentum().e());
      if (std::abs(phishift) > 1e-7) {
        p.rotateZ(phishift * Gaudi::Units::rad);
        child->set_momentum(HepMC::FourVector(p.px(), p.py(), p.pz(), p.e()));
      }
      ++m_particles_processed;
    }
 
 
    // now rotate descendant vertices
    for (HepMC::GenVertexPtr descvtx:  HepMC::descendant_vertices(endvtx)) { //}
 
      ATH_MSG_DEBUG("Processing vertex " << descvtx);
 
      // rotate vertex
      if(std::abs(phishift) > 1e-7) {
        CLHEP::HepLorentzVector position(descvtx->position().x(),
                                         descvtx->position().y(),
                                         descvtx->position().z(),
                                         descvtx->position().t());
        position.rotateZ(phishift*Gaudi::Units::rad);
        descvtx->set_position(HepMC::FourVector( position.x(),position.y(),position.z(),position.t()) );
      }
 
      // now rotate their associated particles
      for (auto descpart: descvtx->particles_out())
      {
        CLHEP::HepLorentzVector momentum(descpart->momentum().px(),
                                          descpart->momentum().py(),
                                          descpart->momentum().pz(),
                                          descpart->momentum().e());
        ATH_MSG_DEBUG("Descendant particle: " << descpart <<
                      " Eta = "               << descpart->momentum().pseudoRapidity() <<
                      " Phi = "               << descpart->momentum().phi() );
        m_particles_processed++;
        // rotate particle
        if(std::abs(phishift) > 1e-7) {
          momentum.rotateZ(phishift*Gaudi::Units::rad);
          descpart->set_momentum( HepMC::FourVector(momentum.px(),momentum.py(),momentum.pz(),momentum.e()) );
          ATH_MSG_DEBUG(" Phi shift =   " << phishift<<
                        " Phi shifted = " << momentum.phi());
        }
      }
    }
  }
  return;
}

msg()

MsgStream & AthCommonMsg< Gaudi::Algorithm >::msg ( ) const

inlineinherited

Definition at line 24 of file AthCommonMsg.h.

                                {
    return this->msgStream();
  }

msgLvl()

bool AthCommonMsg< Gaudi::Algorithm >::msgLvl ( const MSG::Level lvl ) const

inlineinherited

Definition at line 30 of file AthCommonMsg.h.

                                               {
    return this->msgLevel(lvl);
  }

OO_eta_indep()

void AddFlowByShifting::OO_eta_indep ( double b, double eta, double pt )

private

Definition at line 710 of file AddFlowByShifting.cxx.

                                                                 {
  pt=pt/1000.0; //convert to GeV
 
 
  //For Testing
  /*
  m_v_n[Harmonic::v1]=0;
  m_v_n[Harmonic::v2]=0.1;
  m_v_n[Harmonic::v3]=0.05;
  m_v_n[Harmonic::v4]=0.03;
  m_v_n[Harmonic::v5]=0;
  m_v_n[Harmonic::v6]=0;
  return;
  */
  
 
  //pT differential vn for 0-5% centrality
  static TGraph gr_v2pt;
  static TGraph gr_v3pt;
  static TGraph gr_v4pt;
  //centrality dependent scale-factors for other centralities
  static TGraph gr_v2_cent_scale;
  static TGraph gr_v3_cent_scale;
  static TGraph gr_v4_cent_scale;
 
  static bool is_initialized=false;
  if(is_initialized==false){
    //----------------------------------------------------------------------------------
    //The pt Dependent v2,v3,v4 for the 0-5% centrality interval from template-fit method
    //https://atlas.web.cern.ch/Atlas/GROUPS/PHYSICS/PAPERS/HION-2025-02/fig_05.png
    const int NumPtBins=23;
    float v2_pt[NumPtBins]={
      0.0      ,//extrapolation
      0.049533 ,0.0622719,0.0776591,0.0906141,0.101341 ,0.110451 ,0.117545,
      0.123217 ,0.126643 ,0.129142 ,0.129176 ,0.128795 ,0.126009 ,0.120966,
      0.109416 ,0.0966881,0.0825195,0.0708585,0.0615164,0.0508748,
      0.015     ,0.015  //extrapolation to  1.5% v2 at high-pT
    }; 
    float v3_pt[NumPtBins]={
      0.0      ,//extrapolation
      0.0186328,0.0262968,0.0366394,0.0467804,0.0562739,0.0639566,0.0714195,
      0.0783429,0.0838629,0.0886202,0.0907969,0.0922756,0.0948961,0.0947266,
      0.0898434,0.0839913,0.0801961,0.0685197,0.0648406,0.0600871,
      0.00     ,0.00//extrapolation
    };
    //----------------------------------------------------
    //For v3, we use the average of the template and 2PC values
    float v3_pt_2PC[NumPtBins]={
      0.0      ,//extrapolation
      0.0182106,0.0255873,0.0355349,0.0451583,0.0537545,0.0611133,0.0675791,
      0.0732434,0.0778371,0.0812549,0.0823093,0.0826907,0.082954 ,0.0800712,
      0.0719686,0.062255 ,0.0498296,0.0370541,0.027638 ,0.0139785,
      0.00     ,0.00//extrapolation
    };
    for(int i=0;i<NumPtBins;i++) v3_pt[i]=(v3_pt[i]+v3_pt_2PC[i])/2.0;
    //----------------------------------------------------
 
    float v4_pt[NumPtBins]={
      0.0       ,//extrapolation
      0.00446155,0.00710366,0.0117637,0.0158203,0.01931  ,0.0234234,0.0268995,
      0.0299693 ,0.0321662 ,0.0335182,0.0358789,0.0374195,0.0359359,0.0327255,
      0.0254762 ,0.0273239 ,0.0212901,0.0152591, -0      ,-0       ,//-ve values replaced by -0
      0.0       , 0.00//extrapolation
    };
    float pt_bins[NumPtBins]={
      0.0      ,//extrapolation
      0.55     ,0.7      ,0.9      ,1.1      ,      1.3,      1.5,     1.7,
      1.9      ,2.1      ,2.3      ,2.5      ,      2.7,      2.9,    3.25,
      3.75     ,4.25     ,4.75     ,5.25     ,     5.75,      6.5,    
      10.0     ,10000    //extrapolation
    };
    gr_v2pt=TGraph(NumPtBins,pt_bins,v2_pt);
    gr_v3pt=TGraph(NumPtBins,pt_bins,v3_pt);
    gr_v4pt=TGraph(NumPtBins,pt_bins,v4_pt);
    gr_v2pt.SetBit(TGraph::kIsSortedX);//For faster return From Eval()
    gr_v3pt.SetBit(TGraph::kIsSortedX);//For faster return From Eval()
    gr_v4pt.SetBit(TGraph::kIsSortedX);//For faster return From Eval()
    //----------------------------------------------------------------------------------
  
  
    //----------------------------------------------------------------------------------
    //centrality dependent v2, v3, v4 (template-fit method)
    //https://atlas.web.cern.ch/Atlas/GROUPS/PHYSICS/PAPERS/HION-2025-02/figaux_03.png
    const int NumCentBins=23;
    float v2_cent[NumCentBins]={
      0.0749381,//underflow
      0.0749381,0.0782992,0.0801742,0.0814646,0.0826611,0.0853616,0.0889407,
      0.0913229,0.0928072,0.0933774,0.0933299,0.0926254,0.0910148,0.0889322,
      0.086161 ,0.0828212,0.0793132,0.0755036,0.0714726,0.0668264,
      0.055    , 0.05 //80-100 and overflow, put by hand ; using 5% v2 asymptotically
    };
    float v3_cent[NumCentBins]={
      0.0409694,//underflow
      0.0409694,0.0399241,0.0394262,0.0388545,0.0385726,0.0375281,0.0360179,
      0.0345174,0.0330085,0.0314149,0.0298081,0.0282141,0.02664  ,0.0253094,
      0.0239221,0.0220779,0.0213821,0.0203176,0.0175276,0.0151146,
      0.010    , 0.0 //80-100 and overflow, put by hand; using 0% v3 asymptotically 
    };
    float v4_cent[NumCentBins]={
      0.0131927 ,//underflow
      0.0131927 ,0.0129838,0.0135232 ,0.0132689 ,0.0130435 ,0.0134254 ,0.013706 ,
      0.0136456 ,0.0136656,0.013245  ,0.0130861 ,0.0125356 ,0.0118006 ,0.0110316,
      0.00987005,0.0106963,0.00933145,0.00784534,0.00393584,0.00308131,
      0.003    , 0.0 //80-100 and overflow, put by hand; using 0% v4 asymptotically
    };
    //b_imp values for differen centralities
    float b_imp[NumCentBins]={
      -1        ,  //underflow
      0.481843  ,0.886609 , 1.13462  , 1.33981  , 1.52064  , 1.95342  , 2.51881  , 
      2.97665   ,3.37416  , 3.72876  , 4.05377  ,  4.35564 , 4.63868  , 4.91057  , 
      5.17773   ,5.43665  , 5.68739  , 5.94258  , 6.20668  , 6.49519  , 
      7.6415   , 15 , //80-100 and overflow
    };
 
    //Integrated v2, v3, v4, for 0-5% centrality bin
    float vn_0_5[3]={0.078963056, 0.039746672, 0.013202649};
    //divide by integrated vn for 0-5% centrality
    for(int i=0;i<NumCentBins;i++){
      v2_cent[i]/=vn_0_5[0];
      v3_cent[i]/=vn_0_5[1]; 
      v4_cent[i]/=vn_0_5[2];
    }
    gr_v2_cent_scale=TGraph(NumCentBins, b_imp, v2_cent);
    gr_v3_cent_scale=TGraph(NumCentBins, b_imp, v3_cent);
    gr_v4_cent_scale=TGraph(NumCentBins, b_imp, v4_cent);
    gr_v2_cent_scale.SetBit(TGraph::kIsSortedX);//For faster return From Eval()
    gr_v3_cent_scale.SetBit(TGraph::kIsSortedX);//For faster return From Eval()
    gr_v4_cent_scale.SetBit(TGraph::kIsSortedX);//For faster return From Eval()
    //----------------------------------------------------------------------------------
    is_initialized=true;
  }
 
 
  m_v_n[Harmonic::v1]=0;
  m_v_n[Harmonic::v2]=gr_v2pt.Eval(pt) * gr_v2_cent_scale.Eval(b);
  m_v_n[Harmonic::v3]=gr_v3pt.Eval(pt) * gr_v3_cent_scale.Eval(b);
  m_v_n[Harmonic::v4]=gr_v4pt.Eval(pt) * gr_v4_cent_scale.Eval(b);
  m_v_n[Harmonic::v5]=0;
  m_v_n[Harmonic::v6]=0;
}

outputHandles()

virtual std::vector< Gaudi::DataHandle * > AthCommonDataStore< AthCommonMsg< Gaudi::Algorithm > >::outputHandles ( ) const

overridevirtualinherited

Return this algorithm's output handles.

We override this to include handle instances from key arrays if they have not yet been declared. See comments on updateVHKA.

p_Pb_cent_eta_indep()

void AddFlowByShifting::p_Pb_cent_eta_indep ( double b, double eta, double pt )

private

Definition at line 679 of file AddFlowByShifting.cxx.

{
  pt=pt/1000.0; //convert to GeV
 
  float an_val[4][3];
 
  an_val[0][0] = 0.1149;
  an_val[0][1] = 1.181;
  an_val[0][2] = 0.3767;
 
  an_val[1][0] = 0.0498;
  an_val[1][1] = 1.688;
  an_val[1][2] = 0.5046;
 
  an_val[2][0] = 0.02095;
  an_val[2][1] = 2.196;
  an_val[2][2] = 0.6259;
 
  an_val[3][0] = 0.00682*0.5;//added in 0.5 factor by hand
  an_val[3][1] = 4.938;
  an_val[3][2] = 1.237;
 
  m_v_n[Harmonic::v1]=0;
  m_v_n[Harmonic::v2]=an_val[0][0]*std::pow(pt,an_val[0][1])*std::exp(-an_val[0][2]*pt);
  m_v_n[Harmonic::v3]=an_val[1][0]*std::pow(pt,an_val[1][1])*std::exp(-an_val[1][2]*pt);
  m_v_n[Harmonic::v4]=an_val[2][0]*std::pow(pt,an_val[2][1])*std::exp(-an_val[2][2]*pt);
  m_v_n[Harmonic::v5]=an_val[3][0]*std::pow(pt,an_val[3][1])*std::exp(-an_val[3][2]*pt);
  m_v_n[Harmonic::v6]=0;
}

renounce()

std::enable_if_t< std::is_void_v< std::result_of_t< decltype(&T::renounce)(T)> > &&!std::is_base_of_v< SG::VarHandleKeyArray, T > &&std::is_base_of_v< Gaudi::DataHandle, T >, void > AthCommonDataStore< AthCommonMsg< Gaudi::Algorithm > >::renounce ( T & h )

inlineprotectedinherited

Definition at line 380 of file AthCommonDataStore.h.

  {
    h.renounce();
    PBASE::renounce (h);
  }

renounceArray()

void AthCommonDataStore< AthCommonMsg< Gaudi::Algorithm > >::renounceArray ( SG::VarHandleKeyArray & handlesArray )

inlineprotectedinherited

remove all handles from I/O resolution

Definition at line 364 of file AthCommonDataStore.h.

                                                          {
    handlesArray.renounce();
  }

Set_EbE_Fluctuation_Multipliers()

void AddFlowByShifting::Set_EbE_Fluctuation_Multipliers ( HepMC::GenVertexPtr mainvtx, float b, CLHEP::HepRandomEngine * rndmEngine )

private

Definition at line 851 of file AddFlowByShifting.cxx.

                                                                                                                    {
  for(int ihar=0;ihar<Harmonic::NumHar;ihar++){
    m_EbE_Multiplier_vn[ihar]=1.0;
  }
 
 
  for(int ihar=0;ihar<Harmonic::NumHar;ihar++){
    float vn_rp=0,delta=0;//BG parameterizations
 
    //in the following we assume that the vn in the event is sqrt(<vn^2>)
    //This is because the vn(pT) were tuned to 2PC/EP measurements
    //then : vn_evt=vn_rp^2 + 2*(delta^2)
    //which is used together with the "alpha" to get the "vn_rp" and "delta"
 
    //No EbE fluctuation for v1
    if     (ihar==Harmonic::v1) continue; 
    //v2
    else if(ihar==Harmonic::v2){    
      //alpha stores ratio of V2_RP over delta
      float alpha=m_graph_fluc->Eval(b);
      delta=1.0f/sqrt(2.0f+alpha*alpha); //stores delta  /v2{2}
      vn_rp=alpha*delta              ; //stores v2^{RP}/v2{2}       
    }
    //v3-v6
    else{ 
      vn_rp =0;
      delta=1.0f/sqrt(2.0f);
    }
    float X=CLHEP::RandGaussQ::shoot(rndmEngine,vn_rp,delta);
    float Y=CLHEP::RandGaussQ::shoot(rndmEngine,0.0  ,delta);
    m_EbE_Multiplier_vn[ihar]=sqrt(X*X+ Y*Y);
    ATH_MSG_INFO("EbE_Multiplier_v"<<ihar+1<<"="<<m_EbE_Multiplier_vn[ihar]);
  }
}

setFilterPassed()

virtual void AthCommonAlgorithm< Gaudi::Algorithm >::setFilterPassed ( bool state, const EventContext & ctx ) const

inlinevirtualinherited

Set filter decision:

Reimplemented in AthFilterAlgorithm.

Definition at line 99 of file AthCommonAlgorithm.h.

                                                                            {
    execState( ctx ).setFilterPassed( state );
  }

SetParentToRanPhi()

double AddFlowByShifting::SetParentToRanPhi ( HepMC::GenParticlePtr parent, CLHEP::HepRandomEngine * rndmEngine )

private

Definition at line 354 of file AddFlowByShifting.cxx.

{
  // Set particle to random phi
  // Return phi shift
  m_particles_processed++;
 
  double phi, phishift;
  CLHEP::HepLorentzVector momentum(parent->momentum().px(),
                                   parent->momentum().py(),
                                   parent->momentum().pz(),
                                   parent->momentum().e());
  phi = momentum.phi();
 
  double rannum = CLHEP::RandFlat::shoot(rndmEngine);
  double ranphi = (rannum-0.5)*2*M_PI;
  phishift = ranphi - phi;
 
  momentum.setPhi(ranphi*Gaudi::Units::rad);
  parent->set_momentum(  HepMC::FourVector(momentum.px(),momentum.py(),momentum.pz(),momentum.e()) );
 
  ATH_MSG_DEBUG("Parent phi randomized = " << momentum.phi());
 
  return phishift;
}

sysExecute()

StatusCode AthCommonAlgorithm< Gaudi::Algorithm >::sysExecute ( const EventContext & ctx )

overridevirtualinherited

Execute an algorithm.

We override this in order to work around an issue with the Algorithm base class storing the event context in a member variable that can cause crashes in MT jobs.

Reimplemented in AthAnalysisAlgorithm.

Definition at line 80 of file AthCommonAlgorithm.cxx.

{
  return BaseAlg::sysExecute (ctx);
}

sysInitialize()

StatusCode AthCommonAlgorithm< Gaudi::Algorithm >::sysInitialize ( )

overridevirtualinherited

Override sysInitialize.

Override sysInitialize from the base class.

Loop through all output handles, and if they're WriteCondHandles, automatically register them and this Algorithm with the CondSvc

Scan through all outputHandles, and if they're WriteCondHandles, register them with the CondSvc

Reimplemented from AthCommonDataStore< AthCommonMsg< Gaudi::Algorithm > >.

Reimplemented in AthAnalysisAlgorithm, AthFilterAlgorithm, AthHistogramAlgorithm, HypoBase, InputMakerBase, and PyAthena::Alg.

Definition at line 60 of file AthCommonAlgorithm.cxx.

                                                      {
  StatusCode sc=AthCommonDataStore<AthCommonMsg<BaseAlg>>::sysInitialize();
 
  if (sc.isFailure()) {
    return sc;
  }
  
  ServiceHandle<ICondSvc> cs("CondSvc",name());
  for (auto h : outputHandles()) {
    if (h->isCondition() && h->mode() == Gaudi::DataHandle::Writer) {
      // do this inside the loop so we don't create the CondSvc until needed
      if ( cs.retrieve().isFailure() ) {
        ATH_MSG_WARNING("no CondSvc found: won't autoreg WriteCondHandles");
        return StatusCode::SUCCESS;
      }      
      if (cs->regHandle(this,*h).isFailure()) {
        sc = StatusCode::FAILURE;
        ATH_MSG_ERROR("unable to register WriteCondHandle " << h->fullKey()
                      << " with CondSvc");
      }
    }
  }
  return sc;  
}

sysStart()

virtual StatusCode AthCommonDataStore< AthCommonMsg< Gaudi::Algorithm > >::sysStart ( )

overridevirtualinherited

Handle START transition.

We override this in order to make sure that conditions handle keys can cache a pointer to the conditions container.

updateVHKA()

void AthCommonDataStore< AthCommonMsg< Gaudi::Algorithm > >::updateVHKA ( Gaudi::Details::PropertyBase & )

inlineinherited

Definition at line 308 of file AthCommonDataStore.h.

                                               {
    // debug() << "updateVHKA for property " << p.name() << " " << p.toString() 
    //         << "  size: " << m_vhka.size() << endmsg;
    for (auto &a : m_vhka) {
      std::vector<SG::VarHandleKey*> keys = a->keys();
      for (auto k : keys) {
        k->setOwner(this);
      }
    }
  }

vn_func()

double AddFlowByShifting::vn_func ( double x, void * params )

static

Definition at line 40 of file AddFlowByShifting.cxx.

{
  float *par_float = (float*) params;
  double phi_0  = par_float[0];
  float *vn    = par_float+1;
  float *psi_n = vn+Harmonic::NumHar;
  double val=x   +2*(   vn[Harmonic::v1]*sin(1*(x-psi_n[Harmonic::v1]))/1.0 +
                        vn[Harmonic::v2]*sin(2*(x-psi_n[Harmonic::v2]))/2.0 +
                        vn[Harmonic::v3]*sin(3*(x-psi_n[Harmonic::v3]))/3.0 +
                        vn[Harmonic::v4]*sin(4*(x-psi_n[Harmonic::v4]))/4.0 +
                        vn[Harmonic::v5]*sin(5*(x-psi_n[Harmonic::v5]))/5.0 +
                        vn[Harmonic::v6]*sin(6*(x-psi_n[Harmonic::v6]))/6.0 );
  return val-phi_0;
}

Member Data Documentation

m_custom_v1

FloatProperty AddFlowByShifting::m_custom_v1 {this, "custom_v1" , 0.f }

private

Definition at line 92 of file AddFlowByShifting.h.

{this, "custom_v1"         , 0.f               };

m_custom_v2

FloatProperty AddFlowByShifting::m_custom_v2 {this, "custom_v2" , 0.f }

private

Definition at line 93 of file AddFlowByShifting.h.

{this, "custom_v2"         , 0.f               };

m_custom_v3

FloatProperty AddFlowByShifting::m_custom_v3 {this, "custom_v3" , 0.f }

private

Definition at line 94 of file AddFlowByShifting.h.

{this, "custom_v3"         , 0.f               };

m_custom_v4

FloatProperty AddFlowByShifting::m_custom_v4 {this, "custom_v4" , 0.f }

private

Definition at line 95 of file AddFlowByShifting.h.

{this, "custom_v4"         , 0.f               };

m_custom_v5

FloatProperty AddFlowByShifting::m_custom_v5 {this, "custom_v5" , 0.f }

private

Definition at line 96 of file AddFlowByShifting.h.

{this, "custom_v5"         , 0.f               };

m_custom_v6

FloatProperty AddFlowByShifting::m_custom_v6 {this, "custom_v6" , 0.f }

private

Definition at line 97 of file AddFlowByShifting.h.

{this, "custom_v6"         , 0.f               };

m_detStore

StoreGateSvc_t AthCommonDataStore< AthCommonMsg< Gaudi::Algorithm > >::m_detStore

privateinherited

Pointer to StoreGate (detector store by default).

Definition at line 393 of file AthCommonDataStore.h.

m_EbE_Multiplier_vn

float AddFlowByShifting::m_EbE_Multiplier_vn[Harmonic::NumHar]

private

Definition at line 112 of file AddFlowByShifting.h.

m_evtStore

StoreGateSvc_t AthCommonDataStore< AthCommonMsg< Gaudi::Algorithm > >::m_evtStore

privateinherited

Pointer to StoreGate (event store by default).

Definition at line 390 of file AthCommonDataStore.h.

m_extendedExtraObjects

DataObjIDColl AthCommonAlgorithm< Gaudi::Algorithm >::m_extendedExtraObjects

privateinherited

Extra output dependency collection, extended by AthAlgorithmDHUpdate to add symlinks.

Empty if no symlinks were found.

Definition at line 108 of file AthCommonAlgorithm.h.

m_flow_fluctuations

BooleanProperty AddFlowByShifting::m_flow_fluctuations {this, "FlowFluctuations" , false }

private

Definition at line 84 of file AddFlowByShifting.h.

{this, "FlowFluctuations"  , false             };

m_flow_function

void(AddFlowByShifting::* AddFlowByShifting::m_flow_function) (double b, double eta, double pt)

private

Definition at line 60 of file AddFlowByShifting.h.

m_flow_function_name

StringProperty AddFlowByShifting::m_flow_function_name {this, "FlowFunctionName" , "jjia_minbias_new"}

private

Definition at line 81 of file AddFlowByShifting.h.

{this, "FlowFunctionName"  , "jjia_minbias_new"};

m_flow_implementation

StringProperty AddFlowByShifting::m_flow_implementation {this, "FlowImplementation", "exact" }

private

Definition at line 82 of file AddFlowByShifting.h.

{this, "FlowImplementation", "exact"           };

m_flow_implementation_type

int AddFlowByShifting::m_flow_implementation_type {0}

private

Definition at line 83 of file AddFlowByShifting.h.

{0};

m_flow_maxeta

FloatProperty AddFlowByShifting::m_flow_maxeta {this, "FlowMaxEtaCut" , 10.0 }

private

Definition at line 86 of file AddFlowByShifting.h.

{this, "FlowMaxEtaCut"     , 10.0              };

m_flow_maxpt

FloatProperty AddFlowByShifting::m_flow_maxpt {this, "FlowMaxPtCut" , 1000000.f }

private

Definition at line 89 of file AddFlowByShifting.h.

{this, "FlowMaxPtCut"      , 1000000.f         };

m_flow_mineta

FloatProperty AddFlowByShifting::m_flow_mineta {this, "FlowMinEtaCut" , 0.f }

private

Definition at line 87 of file AddFlowByShifting.h.

{this, "FlowMinEtaCut"     , 0.f               };

m_flow_minpt

FloatProperty AddFlowByShifting::m_flow_minpt {this, "FlowMinPtCut" , 0.f }

private

Definition at line 90 of file AddFlowByShifting.h.

{this, "FlowMinPtCut"      , 0.f               };

m_flowb_sw

IntegerProperty AddFlowByShifting::m_flowb_sw {this, "FlowBSwitch" , 0 }

private

Definition at line 91 of file AddFlowByShifting.h.

{this, "FlowBSwitch"       , 0                 };//currently not used

m_floweta_sw

IntegerProperty AddFlowByShifting::m_floweta_sw {this, "FlowEtaSwitch" , 0 }

private

Definition at line 85 of file AddFlowByShifting.h.

{this, "FlowEtaSwitch"     , 0                 };

m_flowpt_sw

IntegerProperty AddFlowByShifting::m_flowpt_sw {this, "FlowPtSwitch" , 0 }

private

Definition at line 88 of file AddFlowByShifting.h.

{this, "FlowPtSwitch"      , 0                 };

m_graph_fluc

TGraph* AddFlowByShifting::m_graph_fluc {}

private

Definition at line 71 of file AddFlowByShifting.h.

{};//TGraph storing the v2_RP/delta Vs b_imp

m_inkey

StringProperty AddFlowByShifting::m_inkey {this, "McTruthKey" , "GEN_EVENT" }

private

Definition at line 78 of file AddFlowByShifting.h.

{this, "McTruthKey"        , "GEN_EVENT"       }; //FIXME use Read/WriteHandles

m_outkey

StringProperty AddFlowByShifting::m_outkey {this, "McFlowKey" , "FLOW_EVENT" }

private

Definition at line 79 of file AddFlowByShifting.h.

{this, "McFlowKey"         , "FLOW_EVENT"      }; //FIXME use Read/WriteHandles

m_particles_processed

int AddFlowByShifting::m_particles_processed {0}

private

Definition at line 99 of file AddFlowByShifting.h.

{0};

m_psi_n

float AddFlowByShifting::m_psi_n[Harmonic::NumHar]

private

Definition at line 110 of file AddFlowByShifting.h.

m_ranphi_sw

IntegerProperty AddFlowByShifting::m_ranphi_sw {this, "RandomizePhi" , 0 }

private

Definition at line 80 of file AddFlowByShifting.h.

{this, "RandomizePhi"      , 0                 };

m_rndmSvc

ServiceHandle<IAthRNGSvc> AddFlowByShifting::m_rndmSvc {this, "RndmSvc", "AthRNGSvc"}

private

Definition at line 75 of file AddFlowByShifting.h.

{this, "RndmSvc", "AthRNGSvc"};

m_v_n

float AddFlowByShifting::m_v_n[Harmonic::NumHar]

private

Definition at line 111 of file AddFlowByShifting.h.

m_varHandleArraysDeclared

bool AthCommonDataStore< AthCommonMsg< Gaudi::Algorithm > >::m_varHandleArraysDeclared

privateinherited

Definition at line 399 of file AthCommonDataStore.h.

m_vhka

std::vector<SG::VarHandleKeyArray*> AthCommonDataStore< AthCommonMsg< Gaudi::Algorithm > >::m_vhka

privateinherited

Definition at line 398 of file AthCommonDataStore.h.

---

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

• AddFlowByShifting.h
• AddFlowByShifting.cxx