21 #include "CLHEP/Random/RandomEngine.h"
22 #include "CLHEP/Random/RandFlat.h"
23 #include "CLHEP/Random/RandGaussQ.h"
24 #include "CLHEP/Vector/LorentzVector.h"
30 #include <gsl/gsl_errno.h>
31 #include <gsl/gsl_math.h>
32 #include <gsl/gsl_roots.h>
34 #include "GaudiKernel/PhysicalConstants.h"
42 float *par_float = (
float*)
params;
43 float phi_0 = par_float[0];
44 float *vn = par_float+1;
46 double val=
x +2*( vn[0]*
sin(1*(
x-psi_n[0]))/1.0 + vn[1]*
sin(2*(
x-psi_n[1]))/2.0 +
47 vn[2]*
sin(3*(
x-psi_n[2]))/3.0 + vn[3]*
sin(4*(
x-psi_n[3]))/4.0 +
48 vn[4]*
sin(5*(
x-psi_n[4]))/5.0 + vn[5]*
sin(6*(
x-psi_n[5]))/6.0 );
55 float *par_float = (
float*)
params;
56 float *vn = par_float+1;
58 double val=1 +2*( vn[0]*
cos(1*(
x-psi_n[0]))/1.0 + vn[1]*
cos(2*(
x-psi_n[1]))/2.0 +
59 vn[2]*
cos(3*(
x-psi_n[2]))/3.0 + vn[3]*
cos(4*(
x-psi_n[3]))/4.0 +
60 vn[4]*
cos(5*(
x-psi_n[4]))/5.0 + vn[5]*
cos(6*(
x-psi_n[5]))/6.0 );
69 for(
int ihar = 0; ihar< 6; ihar++){
78 ATH_MSG_INFO(
">>> AddFlowByShifting from Initialize <<<");
82 ATH_MSG_INFO(
"**********Settings for Afterburner************");
107 ATH_MSG_INFO(
"********************************r*************");
132 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,
133 8.117, 8.767, 9.373, 9.943,10.479,10.991,11.479,11.947,15.00 ,100.0};
135 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,
136 8.767, 9.373, 9.943,10.479,10.991,11.479,11.947,12.399,15.00 ,100.0};
137 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,
138 0.3935,0.4106,0.4310,0.4574,0.4674,0.4873,0.4796,0.4856,0.5130,0.5130};
140 for(
int i=0;
i<21;
i++) bimp_vals[
i]=(b_lo[
i]+b_hi[
i])/2.0;
145 return StatusCode::SUCCESS;
150 const EventContext& ctx)
const
154 rngWrapper->
setSeed( rngName, ctx );
162 const EventContext& ctx = Gaudi::Hive::currentContext();
166 if(
evtStore()->
retrieve(hijing_pars,
"Hijing_event_params").isFailure() ) {
168 return StatusCode::FAILURE;
171 " BPhi = " << hijing_pars->
get_bphi());
182 return StatusCode::FAILURE;
190 for (citr = mcCollptr->
begin(); citr!=mcCollptr->
end(); ++citr) {
191 mcFlowCollptr->
push_back(
new HepMC::GenEvent(*(*citr)));
197 for(
int ihar=0;ihar<6;ihar++){
198 m_psi_n[ihar] =(CLHEP::RandFlat::shoot(rndmEngine)-0.5)*2*
M_PI / (ihar+1);
204 ATH_MSG_DEBUG(
" Psi2 for event : "<<(*hijing_pars).get_psi(2));
209 for (itr = mcFlowCollptr->
begin(); itr!=mcFlowCollptr->
end(); ++itr) {
215 auto mainvtx=(*itr)->vertices().front();
217 int particles_in_event = (*itr)->particles().size();
219 for (
auto parent: mainvtx->particles_out())
221 auto mainvtx=*((*itr)->vertices_begin());
223 int particles_in_event = (*itr)->particles_size();
225 for (
auto parent: *mainvtx)
234 " Eta = " <<
momentum.pseudoRapidity()<<
263 ATH_MSG_INFO(
" Particles in event: " << particles_in_event <<
269 return StatusCode::FAILURE;
271 return StatusCode::SUCCESS;
276 ATH_MSG_INFO(
">>> AddFlowByShifting from finalize <<<");
278 return StatusCode::SUCCESS;
283 CLHEP::HepRandomEngine *rndmEngine)
289 double phi, phishift;
296 double rannum = CLHEP::RandFlat::shoot(rndmEngine);
297 double ranphi = (rannum-0.5)*2*
M_PI;
298 phishift = ranphi -
phi;
314 auto endvtx =
parent->end_vertex();
322 for ( HepMC::GenVertex::vertex_iterator
323 descvtxit = endvtx->vertices_begin(HepMC::descendants);
324 descvtxit != endvtx->vertices_end(HepMC::descendants);
326 auto descvtx = (*descvtxit);
332 if(std::abs(phishift) > 1
e-7) {
333 CLHEP::HepLorentzVector position(descvtx->position().x(),
334 descvtx->position().y(),
335 descvtx->position().z(),
336 descvtx->position().t());
338 descvtx->set_position(HepMC::FourVector( position.x(),position.y(),position.z(),position.t()) );
343 for (
auto descpart: descvtx->particles_out())
345 for (
auto descpart: *descvtx)
349 " Eta = " << descpart->momentum().pseudoRapidity() <<
350 " Phi = " << descpart->momentum().phi() );
365 double pt =
parent->momentum().perp();
366 double eta =
parent->momentum().pseudoRapidity();
367 double phi_0 =
parent->momentum().phi();
372 if(phi_0 !=phi_0) {
ATH_MSG_ERROR(
"ERROR phi of track is not defined");error_=1;}
375 <<
parent->momentum().py()<<
" "
376 <<
parent->momentum().pz()<<
" "
377 <<
parent->momentum().e() <<
" "
378 <<
parent->momentum().m() <<
" ");
385 for(
int ihar = 0; ihar< 6; ihar++){
m_v_n [ihar]=0.0;}
386 double b = hijing_pars->
get_b();
389 for(
int ihar = 0; ihar< 6; ihar++){
410 gsl_root_fsolver *
s = gsl_root_fsolver_alloc (
T);
413 for(
int ipar=0;ipar<13;ipar++) {
params[ipar]=0;}
417 gsl_root_fsolver_set (
s, &
F, x_lo, x_hi);
427 status = gsl_root_fsolver_iterate (
s);
428 phi = gsl_root_fsolver_root (
s);
429 x_lo = gsl_root_fsolver_x_lower (
s);
430 x_hi = gsl_root_fsolver_x_upper (
s);
431 status = gsl_root_test_interval (x_lo, x_hi,0, 0.00001);
433 while (
status == GSL_CONTINUE &&
iter < 1000);
434 gsl_root_fsolver_free (
s);
436 if (
iter>=1000)
return 0;
438 phishift =
phi-phi_0;
441 if(std::abs(phishift) > 1
e-7) {
446 " V2 = " <<
m_v_n[1] <<
447 " V3 = " <<
m_v_n[2] <<
448 " V4 = " <<
m_v_n[3] <<
449 " V5 = " <<
m_v_n[4] <<
450 " V6 = " <<
m_v_n[5] <<
451 " Phi shift = " << phishift <<
452 " Phi shifted = " <<
momentum.phi() );
465 a1=0.4397*
std::exp(-(
b-4.526)*(
b-4.526)/72.0) + 0.636;
467 a3=4.79*0.0001*(
b-0.621)*(
b-10.172)*(
b-23)+1.2;
468 a4=0.135*
std::exp(-0.5*(
b-10.855)*(
b-10.855)/4.607/4.607) +0.0120;
472 float temp3 = 0.01 / (1+
std::exp(-(
pt-4.5)/a3));
476 float fb=0.97 +1.06*
std::exp(-0.5*
b*
b/3.2/3.2);
479 float gb= 1.096 +1.36 *
std::exp(-0.5*
b*
b/3.0/3.0);
480 gb=gb*sqrt(
m_v_n[1]);
494 a1=0.4397*
std::exp(-(
b-4.526)*(
b-4.526)/72.0) + 0.636;
496 a3=4.79*0.0001*(
b-0.621)*(
b-10.172)*(
b-23)+1.2;
497 a4=0.135*
std::exp(-0.5*(
b-10.855)*(
b-10.855)/4.607/4.607) +0.0120;
501 float temp3 = 0.01 / (1+
std::exp(-(
pt-4.5)/a3));
568 an_val[0][0] = 0.1149;
569 an_val[0][1] = 1.181;
570 an_val[0][2] = 0.3767;
572 an_val[1][0] = 0.0498;
573 an_val[1][1] = 1.688;
574 an_val[1][2] = 0.5046;
576 an_val[2][0] = 0.02095;
577 an_val[2][1] = 2.196;
578 an_val[2][2] = 0.6259;
580 an_val[3][0] = 0.00682*0.5;
581 an_val[3][1] = 4.938;
582 an_val[3][2] = 1.237;
594 int Total_Multiplicity=0;
598 for(
auto parent: *mainvtx) {
599 float eta=
parent->momentum().pseudoRapidity();
600 float pT =
parent->momentum().perp();
602 for(
int ihar = 0; ihar< 6; ihar++){
m_v_n [ihar]=0.0;}
604 for(
int ihar = 0; ihar< 6; ihar++){EbE_Vn[ihar] +=
m_v_n [ihar];}
605 Total_Multiplicity++;
609 for(
int ihar = 0; ihar< 6; ihar++){
m_v_n [ihar]=0.0;}
610 if(Total_Multiplicity<=0)
return;
613 for(
int ihar=0;ihar<6;ihar++){
614 EbE_Vn[ihar]/=Total_Multiplicity;
615 float vn_rp=0,delta=0;
621 if(ihar==0)
continue;
629 delta=EbE_Vn[ihar]/sqrt(2.0);
632 float X=CLHEP::RandGaussQ::shoot(rndmEngine,vn_rp,delta);
633 float Y=CLHEP::RandGaussQ::shoot(rndmEngine,0.0 ,delta);