ATLAS Offline Software
AddFlowByShifting.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 // File: Generators/FlowAfterburner/AddFlowByShifting.cxx
6 // Description:
7 // This code is used to introduce particle flow
8 // to particles from generated events
9 //
10 // AuthorList:
11 // Andrzej Olszewski: Initial Code February 2006
12 // 11.10.2006: Add predefined flow function by name
13 
15 
16 #include <set>
17 #include <cmath>
18 
19 // For the Athena-based random numbers
21 #include "CLHEP/Random/RandomEngine.h"
22 #include "CLHEP/Random/RandFlat.h"
23 #include "CLHEP/Random/RandGaussQ.h"
24 #include "CLHEP/Vector/LorentzVector.h"
25 //
28 #include "AtlasHepMC/Relatives.h" //descendant_vertices
29 // gnus scientific library
30 #include <gsl/gsl_errno.h>
31 #include <gsl/gsl_math.h>
32 #include <gsl/gsl_roots.h>
33 
34 #include "GaudiKernel/PhysicalConstants.h"
35 
37 
38 #include "TGraph.h"
39 
40 double AddFlowByShifting::vn_func(double x, void *params)
41 {
42  float *par_float = (float*) params;
43  float phi_0 = par_float[0];
44  float *vn = par_float+1;
45  float *psi_n = vn+6;
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 );
49  return val-phi_0;
50 }
51 
52 
54 {
55  float *par_float = (float*) params;
56  float *vn = par_float+1;
57  float *psi_n = vn+6;
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 );
61  return val;
62 }
63 
64 
65 AddFlowByShifting::AddFlowByShifting(const std::string& name, ISvcLocator* pSvcLocator) :
66  AthAlgorithm(name, pSvcLocator)
67 {
68  m_flow_function= NULL;
69  for(int ihar = 0; ihar< 6; ihar++){
70  m_psi_n[ihar] =0.0;
71  m_v_n [ihar] =0.0;
72  m_EbE_Multiplier_vn[ihar]=1.0;
73  }
74 }
75 
76 
78  ATH_MSG_INFO(">>> AddFlowByShifting from Initialize <<<");
79 
80  ATH_CHECK(m_rndmSvc.retrieve());
81 
82  ATH_MSG_INFO("**********Settings for Afterburner************");
83  ATH_MSG_INFO("McTruthKey : " << m_inkey );
84  ATH_MSG_INFO("McFlowKey : " << m_outkey );
85 
86  ATH_MSG_INFO("FlowFunctionName : " << m_flow_function_name );
87  ATH_MSG_INFO("FlowInplementation : " << m_flow_implementation);
88  ATH_MSG_INFO("FlowFluctuations : " << m_flow_fluctuations );
89 
90  ATH_MSG_INFO("RandomizePhi : " << m_ranphi_sw );
91 
92  ATH_MSG_INFO("FlowEtaSwitch : " << m_floweta_sw );
93  ATH_MSG_INFO("FlowMinEtaCut : " << m_flow_mineta );
94  ATH_MSG_INFO("FlowMaxEtaCut : " << m_flow_maxeta );
95 
96  ATH_MSG_INFO("FlowPtSwitch : " << m_flowpt_sw );
97  ATH_MSG_INFO("FlowMinPtCut : " << m_flow_minpt );
98  ATH_MSG_INFO("FlowMaxPtCut : " << m_flow_maxpt );
99 
100  ATH_MSG_INFO("FlowV1 : " << m_custom_v1 );
101  ATH_MSG_INFO("FlowV2 : " << m_custom_v2 );
102  ATH_MSG_INFO("FlowV3 : " << m_custom_v3 );
103  ATH_MSG_INFO("FlowV4 : " << m_custom_v4 );
104  ATH_MSG_INFO("FlowV5 : " << m_custom_v5 );
105  ATH_MSG_INFO("FlowV6 : " << m_custom_v6 );
106  ATH_MSG_INFO("FlowBSwitch : " << m_flowb_sw );
107  ATH_MSG_INFO("********************************r*************");
108 
109 
110  // Select the flow-implementing function based of the function-choice variable
120 
124 
125 
126 
127  //TGraph storing the v2_RP/delta Vs b_imp values to be used in implementing the EbyE fluctuations
128  //the values below are b_imp-low,b_imp-high, delta/v2_RP for different centralities
129  //underflow and overflow bins are added for smooth extrapolation
130  //The delta/v2_RP values are taken from Fig15 of EbE vn paper (arXiv:1305.2942)
131  // <0 , 0-1 , 1-2 , 2-3 , 3-4 , 4-5 , 5-10 , 10-15, 15-20, 20-25, 25-30,
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};
134  // 30-35, 35-40, 40-45, 45-50, 50-55, 55-60, 60-65, 65-70, 70-
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,//bimp_high
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};
139  float bimp_vals[21];
140  for(int i=0;i<21;i++) bimp_vals[i]=(b_lo[i]+b_hi[i])/2.0;
141  if(m_flow_fluctuations) m_graph_fluc=new TGraph(21,bimp_vals,val);
142 
143 
144  // Initialization terminated
145  return StatusCode::SUCCESS;
146 }
147 
148 
149 CLHEP::HepRandomEngine* AddFlowByShifting::getRandomEngine(const std::string& streamName,
150  const EventContext& ctx) const
151 {
152  ATHRNG::RNGWrapper* rngWrapper = m_rndmSvc->getEngine(this, streamName);
153  std::string rngName = name()+streamName;
154  rngWrapper->setSeed( rngName, ctx );
155  return rngWrapper->getEngine(ctx);
156 }
157 
158 
160  ATH_MSG_INFO(">>> AddFlowByShifting from execute");
161 
162  const EventContext& ctx = Gaudi::Hive::currentContext();
163  CLHEP::HepRandomEngine *rndmEngine = getRandomEngine("FLOW", ctx);
164  // Get hijing event parameters
165  const HijingEventParams *hijing_pars;
166  if( evtStore()->retrieve(hijing_pars, "Hijing_event_params").isFailure() ) {
167  ATH_MSG_ERROR("Could not retrieve Hijing_event_params");
168  return StatusCode::FAILURE;
169  }
170  ATH_MSG_INFO("Event parameters: B = " << hijing_pars->get_b()<<
171  " BPhi = " << hijing_pars->get_bphi());
172 
173 
174  // FIXME: changing data in the event store
175  HijingEventParams *hijing_pars_nc = const_cast<HijingEventParams*> (hijing_pars);
176 
177 
178  // Read Data from Transient Store
179  const McEventCollection* mcCollptr;
180  if ( evtStore()->retrieve(mcCollptr, m_inkey).isFailure() ) {
181  ATH_MSG_ERROR("Could not retrieve truth McEventCollection");
182  return StatusCode::FAILURE;
183  }
184 
185 
186  // Loop over all events in original McEventCollection and
187  // Copy to a new (modifiable) collection
189  McEventCollection* mcFlowCollptr = new McEventCollection();
190  for (citr = mcCollptr->begin(); citr!=mcCollptr->end(); ++citr) {
191  mcFlowCollptr->push_back(new HepMC::GenEvent(*(*citr)));
192  }
193 
194 
195  //Geneate the event-plane angles (some of them may or may not be used later on)
196  //Store the angles into the hijing event parameters
197  for(int ihar=0;ihar<6;ihar++){
198  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
199  hijing_pars_nc->set_psi(ihar+1,m_psi_n[ihar]);
200  }
201  m_psi_n[1]=hijing_pars->get_bphi() ;//the psi2 plane is aligned with the impact parameter
202  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]
203  hijing_pars_nc->set_psi(2,m_psi_n[1]);
204  ATH_MSG_DEBUG(" Psi2 for event : "<<(*hijing_pars).get_psi(2));
205 
206 
207  // Add flow by phi angle shifting
209  for (itr = mcFlowCollptr->begin(); itr!=mcFlowCollptr->end(); ++itr) {
210  ATH_MSG_DEBUG("Next event in the bag ...");
211 
212 
213 
214 #ifdef HEPMC3
215  auto mainvtx=(*itr)->vertices().front();
216  if(m_flow_fluctuations) Set_EbE_Fluctuation_Multipliers(mainvtx,hijing_pars->get_b(),rndmEngine);
217  int particles_in_event = (*itr)->particles().size();
219  for ( auto parent: mainvtx->particles_out())
220 #else
221  auto mainvtx=*((*itr)->vertices_begin());
222  if(m_flow_fluctuations) Set_EbE_Fluctuation_Multipliers(mainvtx,hijing_pars->get_b(),rndmEngine);
223  int particles_in_event = (*itr)->particles_size();
225  for ( auto parent: *mainvtx)
226 #endif
227  {
228  // Process particles from main vertex
229  CLHEP::HepLorentzVector momentum(parent->momentum().px(),
230  parent->momentum().py(),
231  parent->momentum().pz(),
232  parent->momentum().e());
233  ATH_MSG_DEBUG("Parent particle: " << parent <<
234  " Eta = " << momentum.pseudoRapidity()<<
235  " Phi = " << momentum.phi() );
236 
237  //skip particle if eta is outside implementation range
238  if(m_floweta_sw){
239  float eta=std::abs(momentum.pseudoRapidity());
240  if (eta<m_flow_mineta || eta> m_flow_maxeta) continue;
241  }
242 
243  //skip particle if pT is outside implementation range
244  if(m_flowpt_sw){
245  float pT=momentum.perp();
246  if (pT<m_flow_minpt || pT> m_flow_maxpt) continue;
247  }
248 
249  // Randomize phi if explicitely requested
250  if(m_ranphi_sw) {
251  double phishift = SetParentToRanPhi(parent, rndmEngine);
252  MoveDescendantsToParent(parent, phishift) ;// adjust decsandants to parent position
253  }
254 
255  // Add flow to particles from main vertex
256  double phishift = AddFlowToParent(parent, hijing_pars);
257  MoveDescendantsToParent(std::move(parent), phishift);// adjust decsandants to parent position
258  }
259 
260  // correct for double counting
262  // correct for incoming particles
263  ATH_MSG_INFO( " Particles in event: " << particles_in_event <<
264  " Processed for flow: " << m_particles_processed+2);
265  }
266 
267  if(evtStore()->record(mcFlowCollptr, m_outkey).isFailure()){
268  ATH_MSG_ERROR("Could not record flow McEventCollection");
269  return StatusCode::FAILURE;
270  }
271  return StatusCode::SUCCESS;
272 }
273 
274 
276  ATH_MSG_INFO(">>> AddFlowByShifting from finalize <<<");
277  // End of finalization step
278  return StatusCode::SUCCESS;
279 }
280 
281 
283  CLHEP::HepRandomEngine *rndmEngine)
284 {
285  // Set particle to random phi
286  // Return phi shift
288 
289  double phi, phishift;
290  CLHEP::HepLorentzVector momentum(parent->momentum().px(),
291  parent->momentum().py(),
292  parent->momentum().pz(),
293  parent->momentum().e());
294  phi = momentum.phi();
295 
296  double rannum = CLHEP::RandFlat::shoot(rndmEngine);
297  double ranphi = (rannum-0.5)*2*M_PI;
298  phishift = ranphi - phi;
299 
300  momentum.setPhi(ranphi*Gaudi::Units::rad);
301  parent->set_momentum( HepMC::FourVector(momentum.px(),momentum.py(),momentum.pz(),momentum.e()) );
302 
303  ATH_MSG_INFO("Parent phi randomized = " << momentum.phi());
304 
305  return phishift;
306 }
307 
308 
310 (HepMC::GenParticlePtr parent, double phishift)
311 {
312  // Move the branch of descendant vertices and particles
313  // by phishift to parent particle position
314  auto endvtx = parent->end_vertex();
315  if ( endvtx ) {
316  ATH_MSG_DEBUG("Processing branch of parent particle "<< parent);
317 
318  // now rotate descendant vertices
319 #ifdef HEPMC3
320  for (HepMC::GenVertexPtr descvtx: HepMC::descendant_vertices(endvtx)) {
321 #else
322  for ( HepMC::GenVertex::vertex_iterator
323  descvtxit = endvtx->vertices_begin(HepMC::descendants);
324  descvtxit != endvtx->vertices_end(HepMC::descendants);
325  ++descvtxit) {
326  auto descvtx = (*descvtxit);
327 #endif
328 
329  ATH_MSG_DEBUG("Processing vertex " << descvtx);
330 
331  // rotate vertex
332  if(std::abs(phishift) > 1e-7) {
333  CLHEP::HepLorentzVector position(descvtx->position().x(),
334  descvtx->position().y(),
335  descvtx->position().z(),
336  descvtx->position().t());
337  position.rotateZ(phishift*Gaudi::Units::rad);
338  descvtx->set_position(HepMC::FourVector( position.x(),position.y(),position.z(),position.t()) );
339  }
340 
341  // now rotate their associated particles
342 #ifdef HEPMC3
343  for (auto descpart: descvtx->particles_out())
344 #else
345  for (auto descpart: *descvtx)
346 #endif
347  {
348  ATH_MSG_DEBUG("Descendant particle: " << descpart <<
349  " Eta = " << descpart->momentum().pseudoRapidity() <<
350  " Phi = " << descpart->momentum().phi() );
351  }
352  }
353  }
354  return;
355 }
356 
357 
359 {
361  CLHEP::HepLorentzVector momentum(parent->momentum().px(),
362  parent->momentum().py(),
363  parent->momentum().pz(),
364  parent->momentum().e());
365  double pt = parent->momentum().perp();
366  double eta = parent->momentum().pseudoRapidity();
367  double phi_0 = parent->momentum().phi();
368 
369  int error_=0;
370  if(pt !=pt) {ATH_MSG_ERROR("ERROR pt of track is not defined");error_=1;} //true if pt==nan
371  if(eta !=eta) {ATH_MSG_ERROR("ERROR eta of track is not defined");error_=1;}
372  if(phi_0 !=phi_0) {ATH_MSG_ERROR("ERROR phi of track is not defined");error_=1;}
373  if(error_==1){
374  ATH_MSG_ERROR("Original Particle Momentum(px,py,pz,e,m)="<<parent->momentum().px()<<" "
375  <<parent->momentum().py()<<" "
376  <<parent->momentum().pz()<<" "
377  <<parent->momentum().e() <<" "
378  <<parent->momentum().m() <<" ");
379  }
380 
381 
382 
383 
384  //Call the appropriate function to set the vn values
385  for(int ihar = 0; ihar< 6; ihar++){m_v_n [ihar]=0.0;} //reset the vn for this particle
386  double b = hijing_pars->get_b();
387  (*this.*m_flow_function)(b,eta,pt);//get the vn for this particle
388  if(m_flow_fluctuations){//add EbE fluctuations
389  for(int ihar = 0; ihar< 6; ihar++){
390  m_v_n[ihar] *= m_EbE_Multiplier_vn[ihar];
391  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;}
392  }
393  }
394 
395  double phishift=0;
396 
397  // Old fashioned rotation(approximate)
399  float phi=phi_0;
400  phishift= -2*( m_v_n[0]*sin(1*(phi-m_psi_n[0]))/1.0 + m_v_n[1]*sin(2*(phi-m_psi_n[1]))/2.0 +
401  m_v_n[2]*sin(3*(phi-m_psi_n[2]))/3.0 + m_v_n[3]*sin(4*(phi-m_psi_n[3]))/4.0 +
402  m_v_n[4]*sin(5*(phi-m_psi_n[4]))/5.0 + m_v_n[5]*sin(6*(phi-m_psi_n[5]))/6.0 );
403 
404  }
405 
406  // New fashioned rotation(exact)
407  else if (m_flow_implementation_type==1){
408  // Thread-safe according to https://www.gnu.org/software/gsl/doc/html/roots.html
409  const gsl_root_fsolver_type *T ATLAS_THREAD_SAFE = gsl_root_fsolver_brent;
410  gsl_root_fsolver *s = gsl_root_fsolver_alloc (T);
411  double x_lo=-2*M_PI,x_hi=2*M_PI;
412  float params[13];
413  for(int ipar=0;ipar<13;ipar++) {params[ipar]=0;}
414  gsl_function F;
415  F.function = &(AddFlowByShifting::vn_func);
416  F.params =&params;
417  gsl_root_fsolver_set (s, &F, x_lo, x_hi);
418  int iter=0;
419  params[0]=phi_0;
420  params[1]=m_v_n [0]; params[2]=m_v_n [1]; params[3]=m_v_n [2]; params[4 ]=m_v_n [3]; params[5 ]=m_v_n [4]; params[6 ]=m_v_n [5];
421  params[7]=m_psi_n[0]; params[8]=m_psi_n[1]; params[9]=m_psi_n[2]; params[10]=m_psi_n[3]; params[11]=m_psi_n[4]; params[12]=m_psi_n[5];
422  int status;
423  double phi=phi_0;
424  do
425  {
426  iter++;
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);
432  }
433  while (status == GSL_CONTINUE && iter < 1000);
434  gsl_root_fsolver_free (s);
435 
436  if (iter>=1000) return 0;
437 
438  phishift = phi-phi_0;
439  }
440 
441  if(std::abs(phishift) > 1e-7) {
442  momentum.rotateZ(phishift*Gaudi::Units::rad);
443  parent->set_momentum( HepMC::FourVector(momentum.px(),momentum.py(),momentum.pz(),momentum.e()) );
444  }
445  ATH_MSG_DEBUG( "Parent particle: V1 = " << m_v_n[0] <<
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() );
453 
454  return phishift;
455 }
456 
457 
458 
459 // New parameterization for vn
460 void AddFlowByShifting::jjia_minbias_new(double b, double eta, double pt)
461 {
462  pt=pt/1000.0; //convert to GeV
463 
464  float a1,a2,a3,a4;
465  a1=0.4397*std::exp(-(b-4.526)*(b-4.526)/72.0) + 0.636;
466  a2=1.916/(b+2) +0.1;
467  a3=4.79*0.0001*(b-0.621)*(b-10.172)*(b-23)+1.2; // this is >0 for b>0
468  a4=0.135*std::exp(-0.5*(b-10.855)*(b-10.855)/4.607/4.607) +0.0120;
469 
470  float temp1 = std::pow(pt , a1) / (1+std::exp( (pt-3.0)/a3));
471  float temp2 = std::pow(pt+0.1,-a2) / (1+std::exp(-(pt-4.5)/a3));
472  float temp3 = 0.01 / (1+std::exp(-(pt-4.5)/a3));
473 
474  m_v_n[1] = ( a4*(temp1+temp2) + temp3 )* std::exp(-0.5* eta*eta /6.27/6.27) ;
475 
476  float fb=0.97 +1.06*std::exp(-0.5*b*b/3.2/3.2);
477  m_v_n[2]=std::pow(fb*std::sqrt(m_v_n[1]),3);
478 
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]);
481  m_v_n[3]=pow(gb,4);
482  m_v_n[4]=pow(gb,5);
483  m_v_n[5]=pow(gb,6);
484  m_v_n[0]=0;
485 }
486 
487 
488 // New parameterization for vn
490 {
491  pt=pt/1000.0; //convert to GeV
492 
493  float a1,a2,a3,a4;
494  a1=0.4397*std::exp(-(b-4.526)*(b-4.526)/72.0) + 0.636;
495  a2=1.916/(b+2) +0.1;
496  a3=4.79*0.0001*(b-0.621)*(b-10.172)*(b-23)+1.2; // this is >0 for b>0
497  a4=0.135*std::exp(-0.5*(b-10.855)*(b-10.855)/4.607/4.607) +0.0120;
498 
499  float temp1 = std::pow(pt , a1) / (1+std::exp( (pt-3.0)/a3));
500  float temp2 = std::pow(pt+0.1,-a2) / (1+std::exp(-(pt-4.5)/a3));
501  float temp3 = 0.01 / (1+std::exp(-(pt-4.5)/a3));
502 
503  m_v_n[1] = ( a4*(temp1+temp2) + temp3 )* std::exp(-0.5* eta*eta /6.27/6.27) ;
504 
505  m_v_n[0]=0;
506  m_v_n[2]=0;
507  m_v_n[3]=0;
508  m_v_n[4]=0;
509  m_v_n[5]=0;
510 }
511 
512 
513 // Fixed vn
514 void AddFlowByShifting::fixed_vn(double /*b*/, double /*eta*/, double /*pt*/)
515 {
516  m_v_n[0]=0.0000; m_v_n[1]=0.0500;
517  m_v_n[2]=0.0280; m_v_n[3]=0.0130;
518  m_v_n[4]=0.0045; m_v_n[5]=0.0015;
519 }
520 
521 
522 // Fixed 5% v2 (other vn=0)
523 void AddFlowByShifting::fixed_v2(double /*b*/, double /*eta*/, double /*pt*/)
524 {
525  m_v_n[0]=0.0000; m_v_n[1]=0.0500;
526  m_v_n[2]=0.0000; m_v_n[3]=0.0000;
527  m_v_n[4]=0.0000; m_v_n[5]=0.0000;
528 }
529 
530 // Old parameterization for v2
531 void AddFlowByShifting::jjia_minbias_old(double b, double eta, double pt)
532 {
533  m_v_n[0] = 0;
534  m_v_n[1] = 0.03968 * b
535  * (1 - 2.1/(1 + std::exp(1.357*(pt/1000))))
536  * std::exp(-(eta*eta)/(2*6.37*6.37));
537  m_v_n[2]=0.0000; m_v_n[3]=0.0000;
538  m_v_n[4]=0.0000; m_v_n[5]=0.0000;
539 }
540 
541 
542 void AddFlowByShifting::ao_test (double b, double /*eta*/, double pt)
543 {
544  pt/=1000;
545  if(pt>2) pt = 2; // flat max at pt > 2
546  m_v_n[1] = 0.02 * b * pt;
547  m_v_n[0] = 0;
548 }
549 
550 
551 void AddFlowByShifting::custom_vn (double /*b*/, double /*eta*/, double /*pt*/)
552 {
553  m_v_n[0]=m_custom_v1;
554  m_v_n[1]=m_custom_v2;
555  m_v_n[2]=m_custom_v3;
556  m_v_n[3]=m_custom_v4;
557  m_v_n[4]=m_custom_v5;
558  m_v_n[5]=m_custom_v6;
559 }
560 
561 // p_Pb vn
562 void AddFlowByShifting::p_Pb_cent_eta_indep(double /*b*/, double /*eta*/, double pt)
563 {
564  pt=pt/1000.0; //convert to GeV
565 
566  float an_val[4][3];
567 
568  an_val[0][0] = 0.1149;
569  an_val[0][1] = 1.181;
570  an_val[0][2] = 0.3767;
571 
572  an_val[1][0] = 0.0498;
573  an_val[1][1] = 1.688;
574  an_val[1][2] = 0.5046;
575 
576  an_val[2][0] = 0.02095;
577  an_val[2][1] = 2.196;
578  an_val[2][2] = 0.6259;
579 
580  an_val[3][0] = 0.00682*0.5;//added in 0.5 factor by hand
581  an_val[3][1] = 4.938;
582  an_val[3][2] = 1.237;
583 
584  m_v_n[0]=0;
585  m_v_n[1]=an_val[0][0]*std::pow(pt,an_val[0][1])*std::exp(-an_val[0][2]*pt);
586  m_v_n[2]=an_val[1][0]*std::pow(pt,an_val[1][1])*std::exp(-an_val[1][2]*pt);
587  m_v_n[3]=an_val[2][0]*std::pow(pt,an_val[2][1])*std::exp(-an_val[2][2]*pt);
588  m_v_n[4]=an_val[3][0]*std::pow(pt,an_val[3][1])*std::exp(-an_val[3][2]*pt);
589  m_v_n[5]=0;
590 }
591 
592 
593 void AddFlowByShifting::Set_EbE_Fluctuation_Multipliers(HepMC::GenVertexPtr mainvtx,float b, CLHEP::HepRandomEngine *rndmEngine){
594  int Total_Multiplicity=0;
595  double EbE_Vn[6];
596  for(int ihar=0;ihar<6;ihar++){m_EbE_Multiplier_vn[ihar]=1.0;EbE_Vn[ihar]=0.0;}
597 
598  for(auto parent: *mainvtx) {
599  float eta= parent->momentum().pseudoRapidity();
600  float pT = parent->momentum().perp();
601 
602  for(int ihar = 0; ihar< 6; ihar++){m_v_n [ihar]=0.0;}
603  (*this.*m_flow_function)(b,eta,pT);
604  for(int ihar = 0; ihar< 6; ihar++){EbE_Vn[ihar] += m_v_n [ihar];}
605  Total_Multiplicity++;
606  }
607 
608 
609  for(int ihar = 0; ihar< 6; ihar++){m_v_n [ihar]=0.0;}//keep vn as zero before we return from function
610  if(Total_Multiplicity<=0) return;
611 
612 
613  for(int ihar=0;ihar<6;ihar++){
614  EbE_Vn[ihar]/=Total_Multiplicity;
615  float vn_rp=0,delta=0;//BG parameterizations
616 
617  //in the following we assume that the vn in the event is sqrt(<vn^2>)
618  //This is because the vn(pT) were tuned to 2PC/EP measurements
619  //then : vn_evt=vn_rp^2 + 2*(delta^2)
620  //which is used together with the "alpha" to get the "vn_rp" and "delta"
621  if(ihar==0) continue; //No EbE fluctuation for v1
622  else if(ihar==1) { //v2
623  float alpha=1.0/m_graph_fluc->Eval(b);// ratio of V2_RP over delta (from Fig 15 of EbE vn paper)
624  delta=EbE_Vn[ihar]/sqrt(2.0+alpha*alpha);
625  vn_rp=alpha*delta;
626  }
627  else if(ihar>=2) { //v3-v6
628  vn_rp =0;
629  delta=EbE_Vn[ihar]/sqrt(2.0);
630  }
631  if(EbE_Vn[ihar]==0) {ATH_MSG_WARNING("Zero EbeV"<<ihar+1); continue;}
632  float X=CLHEP::RandGaussQ::shoot(rndmEngine,vn_rp,delta);
633  float Y=CLHEP::RandGaussQ::shoot(rndmEngine,0.0 ,delta);
634  m_EbE_Multiplier_vn[ihar]=sqrt(X*X+ Y*Y)/EbE_Vn[ihar];
635  ATH_MSG_INFO("EbE_Multiplier_v"<<ihar+1<<"="<<m_EbE_Multiplier_vn[ihar]);
636  }
637 }
AddFlowByShifting::m_flow_function_name
StringProperty m_flow_function_name
Definition: AddFlowByShifting.h:83
HepMC::GenVertexPtr
HepMC::GenVertex * GenVertexPtr
Definition: GenVertex.h:59
python.PyKernel.retrieve
def retrieve(aClass, aKey=None)
Definition: PyKernel.py:110
AddFlowByShifting::m_flow_function
void(AddFlowByShifting::* m_flow_function)(double b, double eta, double pt)
Definition: AddFlowByShifting.h:61
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
ATHRNG::RNGWrapper::setSeed
void setSeed(const std::string &algName, const EventContext &ctx)
Set the random seed using a string (e.g.
Definition: RNGWrapper.h:169
CalculateHighPtTerm.pT
pT
Definition: ICHEP2016/CalculateHighPtTerm.py:57
AddFlowByShifting::finalize
StatusCode finalize()
Definition: AddFlowByShifting.cxx:275
createLinkingScheme.iter
iter
Definition: createLinkingScheme.py:62
AddFlowByShifting::m_flow_implementation
StringProperty m_flow_implementation
Definition: AddFlowByShifting.h:84
AddFlowByShifting::m_flow_implementation_type
int m_flow_implementation_type
Definition: AddFlowByShifting.h:85
DataModel_detail::const_iterator
Const iterator class for DataVector/DataList.
Definition: DVLIterator.h:82
AddFlowByShifting::Set_EbE_Fluctuation_Multipliers
void Set_EbE_Fluctuation_Multipliers(HepMC::GenVertexPtr mainvtx, float b, CLHEP::HepRandomEngine *rndmEngine)
Definition: AddFlowByShifting.cxx:593
add-xsec-uncert-quadrature-N.alpha
alpha
Definition: add-xsec-uncert-quadrature-N.py:110
phi
Scalar phi() const
phi method
Definition: AmgMatrixBasePlugin.h:67
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
HijingEventParams.h
AddFlowByShifting::m_graph_fluc
TGraph * m_graph_fluc
Definition: AddFlowByShifting.h:71
AddFlowByShifting::getRandomEngine
CLHEP::HepRandomEngine * getRandomEngine(const std::string &streamName, const EventContext &ctx) const
Definition: AddFlowByShifting.cxx:149
eta
Scalar eta() const
pseudorapidity method
Definition: AmgMatrixBasePlugin.h:83
AddFlowByShifting::m_rndmSvc
ServiceHandle< IAthRNGSvc > m_rndmSvc
Definition: AddFlowByShifting.h:75
HepMC::GenParticlePtr
GenParticle * GenParticlePtr
Definition: GenParticle.h:37
AddFlowByShifting::AddFlowToParent
double AddFlowToParent(HepMC::GenParticlePtr parent, const HijingEventParams *hijing_pars)
Definition: AddFlowByShifting.cxx:358
McEventCollection
McEventCollection
Definition: GeneratorObjectsTPCnv.cxx:60
test_pyathena.pt
pt
Definition: test_pyathena.py:11
M_PI
#define M_PI
Definition: ActiveFraction.h:11
AddFlowByShifting::m_custom_v2
FloatProperty m_custom_v2
Definition: AddFlowByShifting.h:98
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
AddFlowByShifting::vn_func_derivative
static double vn_func_derivative(double x, void *params)
Definition: AddFlowByShifting.cxx:53
drawFromPickle.exp
exp
Definition: drawFromPickle.py:36
x
#define x
AddFlowByShifting::MoveDescendantsToParent
void MoveDescendantsToParent(HepMC::GenParticlePtr parent, double phishift)
Definition: AddFlowByShifting.cxx:310
AddFlowByShifting::AddFlowByShifting
AddFlowByShifting(const std::string &name, ISvcLocator *pSvcLocator)
Definition: AddFlowByShifting.cxx:65
AddFlowByShifting::m_custom_v4
FloatProperty m_custom_v4
Definition: AddFlowByShifting.h:100
Monitored::X
@ X
Definition: HistogramFillerUtils.h:24
AddFlowByShifting::m_flowb_sw
IntegerProperty m_flowb_sw
Definition: AddFlowByShifting.h:96
AthCommonDataStore< AthCommonMsg< Algorithm > >::evtStore
ServiceHandle< StoreGateSvc > & evtStore()
The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc.
Definition: AthCommonDataStore.h:85
AddFlowByShifting::m_EbE_Multiplier_vn
float m_EbE_Multiplier_vn[6]
Definition: AddFlowByShifting.h:108
AddFlowByShifting::m_psi_n
float m_psi_n[6]
Definition: AddFlowByShifting.h:107
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
ParticleGun_EoverP_Config.momentum
momentum
Definition: ParticleGun_EoverP_Config.py:63
McEventCollection.h
DataModel_detail::iterator
(Non-const) Iterator class for DataVector/DataList.
Definition: DVLIterator.h:184
lumiFormat.i
int i
Definition: lumiFormat.py:85
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
AddFlowByShifting::m_flow_fluctuations
BooleanProperty m_flow_fluctuations
Definition: AddFlowByShifting.h:86
AddFlowByShifting::m_custom_v6
FloatProperty m_custom_v6
Definition: AddFlowByShifting.h:102
HijingEventParams
Definition: HijingEventParams.h:23
AddFlowByShifting::m_v_n
float m_v_n[6]
Definition: AddFlowByShifting.h:107
AddFlowByShifting::jjia_minbias_old
void jjia_minbias_old(double b, double eta, double pt)
Definition: AddFlowByShifting.cxx:531
AddFlowByShifting::m_flowpt_sw
IntegerProperty m_flowpt_sw
Definition: AddFlowByShifting.h:92
AddFlowByShifting::custom_vn
void custom_vn(double b, double eta, double pt)
Definition: AddFlowByShifting.cxx:551
HijingEventParams::set_psi
void set_psi(int ihar, float psi)
Definition: HijingEventParams.h:107
test_pyathena.parent
parent
Definition: test_pyathena.py:15
AddFlowByShifting::m_inkey
StringProperty m_inkey
Definition: AddFlowByShifting.h:78
AddFlowByShifting::m_flow_mineta
FloatProperty m_flow_mineta
Definition: AddFlowByShifting.h:90
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
McEventCollection
This defines the McEventCollection, which is really just an ObjectVector of McEvent objects.
Definition: McEventCollection.h:32
python.SystemOfUnits.rad
float rad
Definition: SystemOfUnits.py:126
AddFlowByShifting::execute
StatusCode execute()
Definition: AddFlowByShifting.cxx:159
HijingEventParams::get_bphi
float get_bphi() const
Definition: HijingEventParams.h:92
AddFlowByShifting::m_ranphi_sw
IntegerProperty m_ranphi_sw
Definition: AddFlowByShifting.h:81
AthAlgorithm
Definition: AthAlgorithm.h:47
AddFlowByShifting.h
Monitored::Y
@ Y
Definition: HistogramFillerUtils.h:24
WriteHiveWithMetaData.streamName
string streamName
Definition: WriteHiveWithMetaData.py:23
AddFlowByShifting::m_flow_maxeta
FloatProperty m_flow_maxeta
Definition: AddFlowByShifting.h:89
AddFlowByShifting::jjia_minbias_new_v2only
void jjia_minbias_new_v2only(double b, double eta, double pt)
Definition: AddFlowByShifting.cxx:489
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
ATHRNG::RNGWrapper
A wrapper class for event-slot-local random engines.
Definition: RNGWrapper.h:56
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:76
AddFlowByShifting::m_custom_v5
FloatProperty m_custom_v5
Definition: AddFlowByShifting.h:101
AddFlowByShifting::jjia_minbias_new
void jjia_minbias_new(double b, double eta, double pt)
Definition: AddFlowByShifting.cxx:460
AddFlowByShifting::initialize
StatusCode initialize()
Definition: AddFlowByShifting.cxx:77
AddFlowByShifting::m_outkey
StringProperty m_outkey
Definition: AddFlowByShifting.h:79
DataVector::push_back
value_type push_back(value_type pElem)
Add an element to the end of the collection.
ATHRNG::RNGWrapper::getEngine
CLHEP::HepRandomEngine * getEngine(const EventContext &ctx) const
Retrieve the random engine corresponding to the provided EventContext.
Definition: RNGWrapper.h:134
RNGWrapper.h
HijingEventParams::get_b
float get_b() const
Definition: HijingEventParams.h:91
DataVector::end
const_iterator end() const noexcept
Return a const_iterator pointing past the end of the collection.
AddFlowByShifting::fixed_vn
void fixed_vn(double b, double eta, double pt)
Definition: AddFlowByShifting.cxx:514
AddFlowByShifting::vn_func
static double vn_func(double x, void *params)
Definition: AddFlowByShifting.cxx:40
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
Pythia8_RapidityOrderMPI.val
val
Definition: Pythia8_RapidityOrderMPI.py:14
F
#define F(x, y, z)
Definition: MD5.cxx:112
Relatives.h
AddFlowByShifting::m_particles_processed
int m_particles_processed
Definition: AddFlowByShifting.h:104
python.SystemOfUnits.s
float s
Definition: SystemOfUnits.py:147
AddFlowByShifting::m_custom_v1
FloatProperty m_custom_v1
Definition: AddFlowByShifting.h:97
merge.status
status
Definition: merge.py:16
ATLAS_THREAD_SAFE
#define ATLAS_THREAD_SAFE
Definition: checker_macros.h:211
AddFlowByShifting::m_custom_v3
FloatProperty m_custom_v3
Definition: AddFlowByShifting.h:99
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
PowhegControl_ttFCNC_NLO.params
params
Definition: PowhegControl_ttFCNC_NLO.py:226
AddFlowByShifting::p_Pb_cent_eta_indep
void p_Pb_cent_eta_indep(double b, double eta, double pt)
Definition: AddFlowByShifting.cxx:562
AddFlowByShifting::m_floweta_sw
IntegerProperty m_floweta_sw
Definition: AddFlowByShifting.h:88
checker_macros.h
Define macros for attributes used to control the static checker.
pow
constexpr int pow(int base, int exp) noexcept
Definition: ap_fixedTest.cxx:15
AddFlowByShifting::m_flow_minpt
FloatProperty m_flow_minpt
Definition: AddFlowByShifting.h:94
TSU::T
unsigned long long T
Definition: L1TopoDataTypes.h:35
AddFlowByShifting::ao_test
void ao_test(double b, double eta, double pt)
Definition: AddFlowByShifting.cxx:542
DataVector::begin
const_iterator begin() const noexcept
Return a const_iterator pointing at the beginning of the collection.
AddFlowByShifting::m_flow_maxpt
FloatProperty m_flow_maxpt
Definition: AddFlowByShifting.h:93
AddFlowByShifting::fixed_v2
void fixed_v2(double b, double eta, double pt)
Definition: AddFlowByShifting.cxx:523
AddFlowByShifting::SetParentToRanPhi
double SetParentToRanPhi(HepMC::GenParticlePtr parent, CLHEP::HepRandomEngine *rndmEngine)
Definition: AddFlowByShifting.cxx:282