StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StFastGlauberMcMaker.cxx
1 /******************************************************************************
2  * $Id: StFastGlauberMcMaker.cxx,v 1.6 2014/10/21 01:03:38 hmasui Exp $
3  * $Log: StFastGlauberMcMaker.cxx,v $
4  * Revision 1.6 2014/10/21 01:03:38 hmasui
5  * Add 14.5 GeV cross section
6  *
7  * Revision 1.5 2013/02/06 18:58:20 hmasui
8  * Added 2.76 TeV cross section
9  *
10  * Revision 1.4 2013/02/05 22:49:50 hmasui
11  * Added Pb and Cu nuclei
12  *
13  * Revision 1.3 2012/04/25 05:22:42 hmasui
14  * Added deformation parameters, higher order participant eccentricity with r^2 weight
15  *
16  * Revision 1.2 2010/11/20 19:03:11 hmasui
17  * Mode mode flag to StCentralityMaker. Use STAR logger rather than STD iostream
18  *
19  ******************************************************************************/
20 
21 #include <assert.h>
22 
23 #include "TF1.h"
24 #include "TF3.h"
25 #include "TH1.h"
26 #include "TGraph.h"
27 #include "TMath.h"
28 #include "TVector3.h"
29 
30 #include "StMessMgr.h"
31 
32 #include "StMessMgr.h"
33 #include "Nucleon.h"
34 #include "StCentralityMaker/StCentralityMaker.h"
35 #include "StCentralityMaker/StNegativeBinomial.h"
36 #include "StGlauberTree/StGlauberTree.h"
37 #include "StGlauberUtilities/StGlauberUtilities.h"
38 #include "StFastGlauberMcMaker.h"
39 
40 using std::vector ;
41 
42 ClassImp(StFastGlauberMcMaker)
43 
44  const UInt_t StFastGlauberMcMaker::mVersion = 2 ;
45 
46  //____________________________________________________________________________________________________
47  // Default constructor
49  : mEnergy(200), mOutputFileName("fastglaubermc.root"),
50  mInelasticNNCrossSection(4.2) // fm^2 --> 42 mb
51 {
53  const UInt_t A = 197 ;
54  const Double_t R = 6.38 ;
55  const Double_t d = 0.535 ;
56 
57  mMode = 0 ;
58  mIsDeformed[0] = kFALSE ;
59  mIsDeformed[1] = kFALSE ;
60 
61  Init(A, R, d, 0.0, 0.0, A, R, d, 0.0, 0.0);
62 }
63 
64 //____________________________________________________________________________________________________
66  const TString outputFileName,
67  const TString system,
68  const Double_t energy,
69  const TString type,
70  const Bool_t isDeformed
71  )
72 : mEnergy(energy), mOutputFileName(outputFileName),
73  mInelasticNNCrossSection(GetInelasticNNCrossSection(mEnergy, type))
74 {
75  // Set deformation. Spherical + Spherical or Deformed + Deformed collision
76  // Use other constructor if Spherical + Deformed collision is needed
77  mIsDeformed[0] = isDeformed ;
78  mIsDeformed[1] = isDeformed ;
79 
80  // Initialize nuclei
81  if ( system.CompareTo("auau", TString::kIgnoreCase) == 0 ){
82  LOG_INFO << "StFastGlauberMcMaker Initialize AuAu collisions" << endm;
83  InitAuAu(type);
84  }
85  else if ( system.CompareTo("smsm", TString::kIgnoreCase) == 0 ){
86  LOG_INFO << "StFastGlauberMcMaker Initialize SmSm collisions" << endm;
87  InitSmSm(type);
88  }
89  else if ( system.CompareTo("uu", TString::kIgnoreCase) == 0 ){
90  LOG_INFO << "StFastGlauberMcMaker Initialize UU collisions" << endm;
91  InitUU(type);
92  }
93  else if ( system.CompareTo("pbpb", TString::kIgnoreCase) == 0 ){
94  LOG_INFO << "StFastGlauberMcMaker Initialize PbPb collisions" << endm;
95  InitPbPb(type);
96  }
97  else if ( system.CompareTo("cucu", TString::kIgnoreCase) == 0 ){
98  LOG_INFO << "StFastGlauberMcMaker Initialize CuCu collisions" << endm;
99  InitCuCu(type);
100  }
101  else if ( system.CompareTo("zrzr_case1", TString::kIgnoreCase) == 0 ){
102  LOG_INFO << "StFastGlauberMcMaker Initialize ZrZr_Case1 collisions" << endm;
103  InitZrZr(type,1);
104  }
105  else if ( system.CompareTo("ruru_case1", TString::kIgnoreCase) == 0 ){
106  LOG_INFO << "StFastGlauberMcMaker Initialize RuRu_Case1 collisions" << endm;
107  InitRuRu(type,1);
108  }
109  else if ( system.CompareTo("zrzr_case2", TString::kIgnoreCase) == 0 ){
110  LOG_INFO << "StFastGlauberMcMaker Initialize ZrZr_Case2 collisions" << endm;
111  InitZrZr(type,2);
112  }
113  else if ( system.CompareTo("ruru_case2", TString::kIgnoreCase) == 0 ){
114  LOG_INFO << "StFastGlauberMcMaker Initialize RuRu_Case2 collisions" << endm;
115  InitRuRu(type,2);
116  }
117  else if ( system.CompareTo("zrzr_case3", TString::kIgnoreCase) == 0 ){
118  LOG_INFO << "StFastGlauberMcMaker Initialize ZrZr_Case3 collisions" << endm;
119  InitZrZr(type,3);
120  }
121  else if ( system.CompareTo("ruru_case3", TString::kIgnoreCase) == 0 ){
122  LOG_INFO << "StFastGlauberMcMaker Initialize RuRu_Case3 collisions" << endm;
123  InitRuRu(type,3);
124  }
125  else{
126  Error("StFastGlauberMcMaker", "Unknown system %s. abort", system.Data());
127  assert(0);
128  }
129 }
130 
131 //____________________________________________________________________________________________________
133  const TString outputFileName, // Output filename
134  const UInt_t massNumber, // Mass number of nucleus
135  const Double_t radius, // Radius of nucleus
136  const Double_t skinDepth, // Skin depth of nucleus
137  const Double_t beta2, // 2nd order deformation parameter
138  const Double_t beta4, // 4th order deformation parameter
139  const Double_t inelasticCrossSection, // Inelastic NN cross section
140  const Double_t energy // sqrt(sNN)
141  )
142 : mEnergy(energy), mOutputFileName(outputFileName),
143  mInelasticNNCrossSection(inelasticCrossSection)
144 {
145  mMode = 0 ;
146 
147  // Set deformation
148  mIsDeformed[0] = kFALSE ;
149  mIsDeformed[1] = kFALSE ;
150  if(beta2 != 0.0 || beta4 != 0.0 ){
151  mIsDeformed[0] = kTRUE ;
152  mIsDeformed[1] = kTRUE ;
153  }
154 
156  Init(massNumber, radius, skinDepth, beta2, beta4, massNumber, radius, skinDepth, beta2, beta4) ;
157 }
158 
159 //____________________________________________________________________________________________________
161  const TString outputFileName, // Output filename
162  const UInt_t massNumberA, // Mass number of nucleus for nucleus A
163  const Double_t radiusA, // Radius of nucleus for nucleus A
164  const Double_t skinDepthA, // Skin depth of nucleus for nucleus A
165  const Double_t beta2A, // 2nd order deformation parameter for nucleus A
166  const Double_t beta4A, // 4th order deformation parameter for nucleus A
167  const UInt_t massNumberB, // Mass number of nucleus for nucleus B
168  const Double_t radiusB, // Radius of nucleus for nucleus B
169  const Double_t skinDepthB, // Skin depth of nucleus for nucleus B
170  const Double_t beta2B, // 2nd order deformation parameter for nucleus B
171  const Double_t beta4B, // 4th order deformation parameter for nucleus B
172  const Double_t inelasticCrossSection, // Inelastic NN cross section
173  const Double_t energy // sqrt(sNN)
174  )
175 : mEnergy(energy), mOutputFileName(outputFileName),
176  mInelasticNNCrossSection(inelasticCrossSection)
177 {
178  mMode = 0 ;
179 
180  // Set deformation
181  mIsDeformed[0] = kFALSE ;
182  mIsDeformed[1] = kFALSE ;
183  if( beta2A != 0.0 || beta4A != 0.0 ) mIsDeformed[0] = kTRUE ;
184  if( beta2B != 0.0 || beta4B != 0.0 ) mIsDeformed[1] = kTRUE ;
185 
187  Init(massNumberA, radiusA, skinDepthA, beta2A, beta4A, massNumberB, radiusB, skinDepthB, beta2B, beta4B) ;
188 }
189 
190 //____________________________________________________________________________________________________
191 // Default destructor
193 {
194 }
195 
196 //____________________________________________________________________________________________________
197 void StFastGlauberMcMaker::SetRepulsionDistance(const Double_t repulsionDistance)
198 {
199  mRepulsionDistance = repulsionDistance ;
200  LOG_INFO << "StFastGlauberMcMaker::SetRepulsionDistance Set repulsion distance between two nucleons, d = "
201  << mRepulsionDistance << " fm"
202  << endm;
203 }
204 
205 //____________________________________________________________________________________________________
206 Int_t StFastGlauberMcMaker::Clear()
207 {
208  mGlauberTree->Clear() ;
209 
210  for(UInt_t in=0;in<2;in++){
211  for(UInt_t i=0;i<mNucleons[in].size();i++){
212  mNucleons[in][i]->Reset() ;
213  }
214  }
215 
216  return 1 ;
217 }
218 
219 //____________________________________________________________________________________________________
220 Int_t StFastGlauberMcMaker::InitTree()
221 {
222  // Initialize Glauber tree in write mode
223  mGlauberTree = new StGlauberTree(1);
224 
225  // Clear data members in tree
226  Clear() ;
227 
228  // Initialize tree
229  mGlauberTree->Open(mOutputFileName) ;
230 
231  // Histograms
232  for(Int_t in=0;in<4;in++){
233  TString nucleus("Projectile");
234  if( in % 2 == 1 ) nucleus = "Target";
235 
236  TString title(Form("Woods-saxon density profile (%s)", nucleus.Data()));
237  if(in>=2) title = Form("Woods-saxon density profile after repulsion (%s)", nucleus.Data());
238 
239  mhWoodsSaxon[in] = new TH1D(Form("hWoodsSaxon_%d", in), title, 400, 0, 20);
240  mhWoodsSaxon[in]->SetXTitle("r (fm)");
241  }
242 
243  // Sort output histograms/tree
244  mGlauberTree->Sort();
245 
246  return 1;
247 }
248 
249 //____________________________________________________________________________________________________
250 Int_t StFastGlauberMcMaker::Init(
251  const UInt_t massNumberA, // Mass number of nucleus for nucleus A
252  const Double_t radiusA, // Radius of nucleus for nucleus A
253  const Double_t skinDepthA, // Skin depth of nucleus for nucleus A
254  const Double_t beta2A, // 2nd order deformation parameter for nucleus A
255  const Double_t beta4A, // 4th order deformation parameter for nucleus A
256  const UInt_t massNumberB, // Mass number of nucleus for nucleus B
257  const Double_t radiusB, // Radius of nucleus for nucleus B
258  const Double_t skinDepthB, // Skin depth of nucleus for nucleus B
259  const Double_t beta2B, // 2nd order deformation parameter for nucleus B
260  const Double_t beta4B, // 4th order deformation parameter for nucleus B
261  const TString type
262  ){
263  LOG_INFO << "StFastGlauberMcMaker::Init Inelastic N-N cross section = "
264  << mInelasticNNCrossSection * 10.0 << " mb"
265  << " for sqrt(sNN) = " << mEnergy << " (GeV)"
266  << endm;
267 
268  mDebug = 0 ;
269 
270  mNeventsThrow = 0 ; // Number of all events
271  mNeventsAccept = 0 ; // Number of accepted events (Ncoll>0)
272 
273  // Centrlaity maker
274  const Char_t* system(Form("%s%s_%dGeV", GetName(massNumberA), GetName(massNumberB),
275  static_cast<Int_t>(mEnergy))) ;
276  LOG_INFO << "StFastGlauberMcMaker::init Initialize system " << system << endm;
277  mCentralityMaker = new StCentralityMaker(system);
278 
280  mRepulsionDistance = 0.0 ;
281 
282  mHardCoreSmearing = kFALSE ;
283  mGaussianSmearing = kFALSE ;
284 
285  mCollisionProfile = mkHardCoreProfile ;
286  if ( type.CompareTo("gauss", TString::kIgnoreCase) == 0 ){
287  //DoGaussianCollision() ;
288  }
289 
291  const Double_t dmaxh = TMath::Sqrt(mInelasticNNCrossSection/TMath::Pi());
292  mfHardCore = new TF3("fHardCore", GlauberUtilities::StepFunction, -dmaxh, dmaxh, -dmaxh, dmaxh, -dmaxh, dmaxh, 1);
293  mfHardCore->SetParameter(0, mInelasticNNCrossSection);
294 
295  if(mDebug){
296  LOG_INFO << "StFastGlauberMcMaker::Init Initialize hard-core smearing function (will be used only if hard-core smearing is ON)"
297  << endm;
298  }
299 
301  const Double_t sigma = 0.79/TMath::Sqrt(3.0) ;
302  const Double_t dmaxg = sigma*5.0 ;
303  mfGaussian = new TF3("fGaussian", GlauberUtilities::Gaussian, -dmaxg, dmaxg, -dmaxg, dmaxg, -dmaxg, dmaxg, 1);
304  mfGaussian->SetParameter(0, sigma); // width
305 
306  if(mDebug){
307  LOG_INFO << "StFastGlauberMcMaker::Init Initialize gaussian smearing function (will be used only if gaussian smearing is ON)"
308  << endm;
309  }
310 
311  for(UInt_t in=0; in<2; in++){
312  if(mIsDeformed[in]){
313  // Make sure deformation is finite
314  const Bool_t isNoDeformation = (in==0) ? (beta2A==0.0 && beta4A==0.0) : (beta2B==0.0 && beta4B==0.0) ;
315  if( isNoDeformation ){
316  Error("StFastGlauberMcMaker::Init", "No deformation: beta2=beta4=0. abort");
317  assert(0);
318  }
319 
321  mfWoodsSaxon[in] = 0 ; // NULL pointer for spherical woods-saxon
322 
323  mfWoodsSaxon2D[in] = new TF2(Form("fWoodsSaxon2D_%d", in), GlauberUtilities::WoodsSaxon2D, 0, 20, -1.0, 1.0, 4);
324  mfWoodsSaxon2D[in]->SetParName(0, "Radius");
325  mfWoodsSaxon2D[in]->SetParName(1, "Skin depth");
326  mfWoodsSaxon2D[in]->SetParName(2, "#beta_{2}");
327  mfWoodsSaxon2D[in]->SetParName(3, "#beta_{4}");
328  // mfWoodsSaxon2D[in]->SetParName(4, "#beta_{6}");
329 
330  // Default number of sampling points are too small (probably for y-axis) in TF2
331  // this causes the step-like structures if you get (r,theta) from TF2::GetRandom2()
332  // --> Increase Npx, Npy by a factor of 2
333  mfWoodsSaxon2D[in]->SetNpx(200);
334  mfWoodsSaxon2D[in]->SetNpy(200);
335 
336  if( in == 0 ){
337  mfWoodsSaxon2D[in]->SetParameter(0, radiusA) ;
338  mfWoodsSaxon2D[in]->SetParameter(1, skinDepthA) ;
339  mfWoodsSaxon2D[in]->SetParameter(2, beta2A) ;
340  mfWoodsSaxon2D[in]->SetParameter(3, beta4A) ;
341  }
342  else{
343  mfWoodsSaxon2D[in]->SetParameter(0, radiusB) ;
344  mfWoodsSaxon2D[in]->SetParameter(1, skinDepthB) ;
345  mfWoodsSaxon2D[in]->SetParameter(2, beta2B) ;
346  mfWoodsSaxon2D[in]->SetParameter(3, beta4B) ;
347  }
348 
349  LOG_INFO << "StFastGlauberMcMaker::Init Initialize Deformed Woods-saxon density profile" << endm;
350  for(Int_t ipar=0; ipar<mfWoodsSaxon2D[in]->GetNpar(); ipar++){
351  LOG_INFO << "StFastGlauberMcMaker::Init par[" << ipar << "] = " << mfWoodsSaxon2D[in]->GetParameter(ipar)
352  << ", parName = " << mfWoodsSaxon2D[in]->GetParName(ipar)
353  << endm;
354  }
355  }
356  else{
358  // Radius and skin depth
359  mfWoodsSaxon2D[in] = 0 ; // NULL pointer for deformed woods-saxon
360 
361  mfWoodsSaxon[in] = new TF1(Form("fWoodsSaxon_%d", in), GlauberUtilities::WoodsSaxon, 0, 20, 2);
362  mfWoodsSaxon[in]->SetParName(0, "Radius");
363  mfWoodsSaxon[in]->SetParName(1, "Skin depth");
364  if( in == 0 ){
365  mfWoodsSaxon[in]->SetParameter(0, radiusA) ;
366  mfWoodsSaxon[in]->SetParameter(1, skinDepthA) ;
367  }
368  else{
369  mfWoodsSaxon[in]->SetParameter(0, radiusB) ;
370  mfWoodsSaxon[in]->SetParameter(1, skinDepthB) ;
371  }
372 
373  LOG_INFO << "StFastGlauberMcMaker::Init Initialize Spherical Woods-saxon density profile" << endm;
374  for(Int_t ipar=0; ipar<mfWoodsSaxon[in]->GetNpar(); ipar++){
375  LOG_INFO << "StFastGlauberMcMaker::Init par[" << ipar << "] = " << mfWoodsSaxon[in]->GetParameter(ipar)
376  << ", parName = " << mfWoodsSaxon[in]->GetParName(ipar)
377  << endm;
378  }
379  }
380 
382  mNucleons[in].clear() ;
383  const UInt_t A = (in==0) ? massNumberA : massNumberB ;
384  LOG_INFO << "StFastGlauberMcMaker::Init Initialize " << A << " nucleons ..." << endm;
385 
386  for(UInt_t i=0;i<A;i++){
387  mNucleons[in].push_back( new Nucleon() ) ;
388  }
389  }
390 
391  // Initialize output tree
392  InitTree() ;
393 
394  return 1 ;
395 }
396 
397 //____________________________________________________________________________________________________
398 Int_t StFastGlauberMcMaker::InitAuAu(const TString type)
399 {
401  const UInt_t A = 197 ;
402  const Double_t R = 6.38 ;
403  const Double_t d = 0.535 ;
404 
405  // Parameterization with hard-core smearing
406  // from PRC79, 064904 (2009) T. Hirano et. al.
407  // const Double_t R = 6.42 ;
408  // const Double_t d = 0.544 ;
409  Double_t RError = 0.0 ; // Absolute error on R
410  Double_t dError = 0.0 ; // Absolute error on d
411 
412  // Deformation parameters
413  Double_t beta2 = 0.0 ;
414  Double_t beta4 = 0.0 ;
415  // Double_t beta6 = 0.0 ;
416 
417  // Deformation
418  if( mIsDeformed[0] || mIsDeformed[1] ){
419  LOG_INFO << "StFastGlauberMcMaker::InitAuAu Set deformation for Au nuclei" << endm ;
420  beta2 = -0.131 ; // Deformation parameter beta2
421  beta4 = -0.031 ; // Deformation parameter beta4
422  // beta6 = 0.007 ; // Deformation parameter beta6
423  }
424 
425  // Error on R and d is taken from Atomic Data and Nuclear Table 36, 495 (1987)
426  // R = 6.38 +/- 0.06, d = 0.535 +/- 0.027
427  // --> take 2 sigma
428  if ( type.CompareTo("large", TString::kIgnoreCase) == 0 ){
429  LOG_INFO << "StFastGlauberMcMaker::InitAuAu Set large R, small d" << endm;
430  RError = 0.12 ;
431  dError = -0.054 ;
432  // RError = 0.32 ; //5%
433  // dError = -0.0267 ;//5%
434  }
435  else if ( type.CompareTo("small", TString::kIgnoreCase) == 0 ){
436  LOG_INFO << "StFastGlauberMcMaker::InitAuAu Set small R, large d" << endm;
437  RError = -0.12 ;
438  dError = 0.054 ;
439  // RError = -0.32 ; //5%
440  // dError = 0.0267 ; //5%
441  }
442  else if ( type.CompareTo("gauss", TString::kIgnoreCase) == 0 ){
443  LOG_INFO << "StFastGlauberMcMaker::InitAuAu Gaussian collision profile" << endm;
445  }
446  else if ( type.CompareTo("smallNpp", TString::kIgnoreCase) == 0 ){
447  mMode = 1 ;
448  LOG_INFO << "StFastGlauberMcMaker::InitAuAu Use small npp, large x for multiplicity (mMode = "
449  << mMode << ")"
450  << endm;
451  }
452  else if ( type.CompareTo("largeNpp", TString::kIgnoreCase) == 0 ){
453  mMode = 2 ;
454  LOG_INFO << "StFastGlauberMcMaker::InitAuAu Use large npp, small x for multiplicity (mMode = "
455  << mMode << ")"
456  << endm;
457  }
458 
459  return Init(A, R+RError, d+dError, beta2, beta4, A, R+RError, d+dError, beta2, beta4, type);
460 }
461 
462 //____________________________________________________________________________________________________
463 Int_t StFastGlauberMcMaker::InitSmSm(const TString type)
464 {
466  // - Sm144 (spherical, A=144, N=82, Z=62)
467  // - Sm154 (deformed, A=154, N=92, Z=62)
468  //
469  // Sm144
470  // R = 5.71 = 1.12*(144)^{1/3} - 0.86*(144)^{-1/3}
471  // d = 0.543 = 0.522 * 5.9387 / 5.71
472  const UInt_t A = (mIsDeformed) ? 154 : 144 ;
473  // const Double_t R = (mIsDeformed) ? 5.9387 : 5.71 ;
474  // const Double_t d = (mIsDeformed) ? 0.522 : 0.543 ;
475  const Double_t R = (mIsDeformed) ? 5.9387 : 5.9387 ;
476  const Double_t d = (mIsDeformed) ? 0.522 : 0.522 ;
477  Double_t RError = 0.0 ; // Absolute error on R
478  Double_t dError = 0.0 ; // Absolute error on d
479 
480  // Deformation parameters
481  Double_t beta2 = 0.0 ;
482  Double_t beta4 = 0.0 ;
483  // Double_t beta6 = 0.0 ;
484 
485  // Deformation
486  if( mIsDeformed ){
487  LOG_INFO << "StFastGlauberMcMaker::InitSmSm Set deformation for Sm nuclei" << endm ;
488  beta2 = 0.270 ; // Deformation parameter beta2
489  beta4 = 0.113 ; // Deformation parameter beta4
490  // beta6 = -0.005 ; // Deformation parameter beta6
491  }
492 
493  // Error on R and d is taken from Atomic Data and Nuclear Table 36, 495 (1987)
494  // R = 5.9387 (no error, assign +/- 0.0050), d = 0.522 +/- 0.015 for Sm154
495  // No R and d parameters for Sm144
496  // --> take 2 sigma
497  if ( type.CompareTo("large", TString::kIgnoreCase) == 0 ){
498  LOG_INFO << "StFastGlauberMcMaker::InitSmSm Set large R, small d" << endm;
499  RError = 0.01 ;
500  dError = -0.030 ;
501  }
502  else if ( type.CompareTo("small", TString::kIgnoreCase) == 0 ){
503  LOG_INFO << "StFastGlauberMcMaker::InitSmSm Set small R, large d" << endm;
504  RError = -0.01 ;
505  dError = 0.030 ;
506  }
507  else if ( type.CompareTo("gauss", TString::kIgnoreCase) == 0 ){
508  LOG_INFO << "StFastGlauberMcMaker::InitSmSm Gaussian collision profile" << endm;
510  }
511  else if ( type.CompareTo("smallNpp", TString::kIgnoreCase) == 0 ){
512  mMode = 1 ;
513  LOG_INFO << "StFastGlauberMcMaker::InitSmSm Use small npp, large x for multiplicity (mMode = "
514  << mMode << ")"
515  << endm;
516  }
517  else if ( type.CompareTo("largeNpp", TString::kIgnoreCase) == 0 ){
518  mMode = 2 ;
519  LOG_INFO << "StFastGlauberMcMaker::InitSmSm Use large npp, small x for multiplicity (mMode = "
520  << mMode << ")"
521  << endm;
522  }
523 
524  return Init(A, R+RError, d+dError, beta2, beta4, A, R+RError, d+dError, beta2, beta4, type);
525 }
526 
527 //____________________________________________________________________________________________________
528 Int_t StFastGlauberMcMaker::InitUU(const TString type)
529 {
531  const UInt_t A = 238 ;
532  const Double_t R = 6.81 ;
533  const Double_t d = 0.55 ;
534  Double_t RError = 0.0 ; // Absolute error on R
535  Double_t dError = 0.0 ; // Absolute error on d
536 
537  // Deformation parameters
538  Double_t beta2 = 0.0 ;
539  Double_t beta4 = 0.0 ;
540  // Double_t beta6 = 0.0 ;
541 
542  // Deformation
543  if( mIsDeformed ){
544  LOG_INFO << "StFastGlauberMcMaker::InitUU Set deformation for U nuclei" << endm ;
545  beta2 = 0.28 ; // Deformation parameter beta2
546  beta4 = 0.093 ; // Deformation parameter beta4
547  // beta6 = 0.007 ; // Deformation parameter beta6
548  }
549 
550  // Error on R and d is taken from Atomic Data and Nuclear Table 36, 495 (1987)
551  // R = 6.8054, d = 0.605 +/- 0.016 (PRC13, 1083, 1976)
552  // R = 6.874, d = 0.556 (unpublished)
553  // difference gives
554  // \Delta R = 0.0686 fm, \Delta d = 0.049 fm
555  // --> take 2 sigma
556  if ( type.CompareTo("large", TString::kIgnoreCase) == 0 ){
557  LOG_INFO << "StFastGlauberMcMaker::InitUU Set large R, small d" << endm;
558  RError = 0.137 ;
559  dError = -0.098 ;
560  }
561  else if ( type.CompareTo("small", TString::kIgnoreCase) == 0 ){
562  LOG_INFO << "StFastGlauberMcMaker::InitUU Set small R, large d" << endm;
563  RError = -0.137 ;
564  dError = 0.098 ;
565  }
566  else if ( type.CompareTo("gauss", TString::kIgnoreCase) == 0 ){
567  LOG_INFO << "StFastGlauberMcMaker::InitUU Gaussian collision profile" << endm;
569  }
570  else if ( type.CompareTo("smallNpp", TString::kIgnoreCase) == 0 ){
571  mMode = 1 ;
572  LOG_INFO << "StFastGlauberMcMaker::InitUU Use small npp, large x for multiplicity (mMode = "
573  << mMode << ")"
574  << endm;
575  }
576  else if ( type.CompareTo("largeNpp", TString::kIgnoreCase) == 0 ){
577  mMode = 2 ;
578  LOG_INFO << "StFastGlauberMcMaker::InitUU Use large npp, small x for multiplicity (mMode = "
579  << mMode << ")"
580  << endm;
581  }
582 
583  return Init(A, R+RError, d+dError, beta2, beta4, A, R+RError, d+dError, beta2, beta4, type);
584 }
585 
586 //____________________________________________________________________________________________________
587 Int_t StFastGlauberMcMaker::InitPbPb(const TString type)
588 {
590  const UInt_t A = 208 ;
591  const Double_t R = 6.62 ;
592  const Double_t d = 0.546 ;
593  Double_t RError = 0.0 ; // Absolute error on R
594  Double_t dError = 0.0 ; // Absolute error on d
595 
596  // Deformation parameters
597  Double_t beta2 = 0.0 ;
598  Double_t beta4 = 0.0 ;
599  // Double_t beta6 = 0.0 ;
600 
601  // Deformation
602  if( mIsDeformed ){
603  LOG_INFO << "StFastGlauberMcMaker::InitPbPb Set deformation for Pb nuclei" << endl;
604  beta2 = 0.0;
605  beta4 = 0.0;
606  }
607 
608  if ( type.CompareTo("large", TString::kIgnoreCase) == 0 ){
609  LOG_INFO << "StFastGlauberMcMaker::InitPbPb Set large R, small d" << endm;
610  RError = 0.120 ;
611  dError = -0.020 ;
612  }
613  else if ( type.CompareTo("small", TString::kIgnoreCase) == 0 ){
614  LOG_INFO << "StFastGlauberMcMaker::InitPbPb Set small R, large d" << endm;
615  RError = -0.120 ;
616  dError = 0.020 ;
617  }
618  else if ( type.CompareTo("gauss", TString::kIgnoreCase) == 0 ){
619  LOG_INFO << "StFastGlauberMcMaker::InitPbPb Gaussian collision profile" << endm;
621  }
622  else if ( type.CompareTo("smallNpp", TString::kIgnoreCase) == 0 ){
623  mMode = 1 ;
624  LOG_INFO << "StFastGlauberMcMaker::InitPbPb Use small npp, large x for multiplicity (mMode = "
625  << mMode << ")"
626  << endm;
627  }
628  else if ( type.CompareTo("largeNpp", TString::kIgnoreCase) == 0 ){
629  mMode = 2 ;
630  LOG_INFO << "StFastGlauberMcMaker::InitPbPb Use large npp, small x for multiplicity (mMode = "
631  << mMode << ")"
632  << endm;
633  }
634 
635  return Init(A, R+RError, d+dError, beta2, beta4, A, R+RError, d+dError, beta2, beta4, type);
636 }
637 
638 //____________________________________________________________________________________________________
639 Int_t StFastGlauberMcMaker::InitCuCu(const TString type)
640 {
642  const UInt_t A = 63 ;
643  const Double_t R = 4.19 ;
644  const Double_t d = 0.60 ;
645  Double_t RError = 0.0 ; // Absolute error on R
646  Double_t dError = 0.0 ; // Absolute error on d
647 
648  // Deformation parameters
649  Double_t beta2 = 0.0 ;
650  Double_t beta4 = 0.0 ;
651  // Double_t beta6 = 0.0 ;
652 
653  // Deformation
654  if( mIsDeformed ){
655  LOG_INFO << "StFastGlauberMcMaker::InitCuCu Set deformation for Cu nuclei" << endm ;
656  beta2 = 0.162 ; // Deformation parameter beta2
657  beta4 = -0.006 ; // Deformation parameter beta4
658  }
659 
660  // Error on R and d is taken from Atomic Data and Nuclear Table 36, 495 (1987)
661  if ( type.CompareTo("large", TString::kIgnoreCase) == 0 ){
662  LOG_INFO << "StFastGlauberMcMaker::InitCuCu Set large R, small d" << endm;
663  RError = 0.041 ;
664  dError = -0.016 ;
665  }
666  else if ( type.CompareTo("small", TString::kIgnoreCase) == 0 ){
667  LOG_INFO << "StFastGlauberMcMaker::InitCuCu Set small R, large d" << endm;
668  RError = -0.041 ;
669  dError = 0.016 ;
670  }
671  else if ( type.CompareTo("gauss", TString::kIgnoreCase) == 0 ){
672  LOG_INFO << "StFastGlauberMcMaker::InitCuCu Gaussian collision profile" << endm;
674  }
675  else if ( type.CompareTo("smallNpp", TString::kIgnoreCase) == 0 ){
676  mMode = 1 ;
677  LOG_INFO << "StFastGlauberMcMaker::InitCuCu Use small npp, large x for multiplicity (mMode = "
678  << mMode << ")"
679  << endm;
680  }
681  else if ( type.CompareTo("largeNpp", TString::kIgnoreCase) == 0 ){
682  mMode = 2 ;
683  LOG_INFO << "StFastGlauberMcMaker::InitCuCu Use large npp, small x for multiplicity (mMode = "
684  << mMode << ")"
685  << endm;
686  }
687 
688  return Init(A, R+RError, d+dError, beta2, beta4, A, R+RError, d+dError, beta2, beta4, type);
689 }
690 
691 //____________________________________________________________________________________________________
692 Int_t StFastGlauberMcMaker::InitZrZr(const TString type, int Case)
693 {
695  const UInt_t A = 96 ;
696  Double_t R;
697  Double_t d;
698  if(Case==1){R = 5.02 ; d = 0.46 ;}
699  else if(Case==2){R = 5.02 ; d = 0.46 ;}
700  else if(Case==3){R = 4.965 ; d = 0.556 ;}
701  Double_t RError = 0.0 ; // Absolute error on R
702  Double_t dError = 0.0 ; // Absolute error on d
703 
704  // Deformation parameters
705  Double_t beta2 = 0.0 ;
706  Double_t beta4 = 0.0 ;
707  // Double_t beta6 = 0.0 ;
708 
709  // Deformation parameters taken from https://journals.aps.org/prc/pdf/10.1103/PhysRevC.99.044908
710  // B Schenke, C Shen, P Tribedy. Phys. Rev. C 99, 044908 (2019)
711  // Deformation
712  if( mIsDeformed ){
713  LOG_INFO << "StFastGlauberMcMaker::InitZrZr Set deformation for Zr nuclei" << endm ;
714  if(Case==1){beta2 = 0.08 ; beta4 = 0.0 ; }
715  else if(Case==2){beta2 = 0.217 ; beta4 = 0.0 ; }
716  else if(Case==3){beta2 = 0.0 ; beta4 = 0.0 ; }
717  }
718 
719  if ( type.CompareTo("large", TString::kIgnoreCase) == 0 ){
720  LOG_INFO << "StFastGlauberMcMaker::InitZrZr Set large R, small d" << endm;
721  RError = 0.0 ;
722  dError = 0.0 ;
723  }
724  else if ( type.CompareTo("small", TString::kIgnoreCase) == 0 ){
725  LOG_INFO << "StFastGlauberMcMaker::InitZrZr Set small R, large d" << endm;
726  RError = 0.0 ;
727  dError = 0.0 ;
728  }
729  else if ( type.CompareTo("gauss", TString::kIgnoreCase) == 0 ){
730  LOG_INFO << "StFastGlauberMcMaker::InitZrZr Gaussian collision profile" << endm;
732  }
733  else if ( type.CompareTo("smallNpp", TString::kIgnoreCase) == 0 ){
734  mMode = 1 ;
735  LOG_INFO << "StFastGlauberMcMaker::InitZrZr Use small npp, large x for multiplicity (mMode = "
736  << mMode << ")"
737  << endm;
738  }
739  else if ( type.CompareTo("largeNpp", TString::kIgnoreCase) == 0 ){
740  mMode = 2 ;
741  LOG_INFO << "StFastGlauberMcMaker::InitZrZr Use large npp, small x for multiplicity (mMode = "
742  << mMode << ")"
743  << endm;
744  }
745 
746  return Init(A, R+RError, d+dError, beta2, beta4, A, R+RError, d+dError, beta2, beta4, type);
747 }
748 
749 //____________________________________________________________________________________________________
750 Int_t StFastGlauberMcMaker::InitRuRu(const TString type, int Case)
751 {
753  const UInt_t A = 96 ;
754  Double_t R;
755  Double_t d;
756  if(Case==1){R = 5.085 ; d = 0.46 ;}
757  else if(Case==2){R = 5.085 ; d = 0.46 ;}
758  else if(Case==3){R = 5.067 ; d = 0.500 ;}
759  Double_t RError = 0.0 ; // Absolute error on R
760  Double_t dError = 0.0 ; // Absolute error on d
761 
762  // Deformation parameters
763  Double_t beta2 = 0.0 ;
764  Double_t beta4 = 0.0 ;
765  // Double_t beta6 = 0.0 ;
766 
767  // Deformation parameters taken from https://journals.aps.org/prc/pdf/10.1103/PhysRevC.99.044908
768  // B Schenke, C Shen, P Tribedy. Phys. Rev. C 99, 044908 (2019)
769  // Deformation
770  if( mIsDeformed ){
771  LOG_INFO << "StFastGlauberMcMaker::InitRuRu Set deformation for Ru nuclei" << endm ;
772  if(Case==1){beta2 = 0.158 ; beta4 = 0.0 ; }
773  else if(Case==2){beta2 = 0.053 ; beta4 = 0.0 ; }
774  else if(Case==3){beta2 = 0.0 ; beta4 = 0.0 ; }
775  }
776 
777  if ( type.CompareTo("large", TString::kIgnoreCase) == 0 ){
778  LOG_INFO << "StFastGlauberMcMaker::InitRuRu Set large R, small d" << endm;
779  RError = 0.0 ;
780  dError = 0.0 ;
781  }
782  else if ( type.CompareTo("small", TString::kIgnoreCase) == 0 ){
783  LOG_INFO << "StFastGlauberMcMaker::InitRuRu Set small R, large d" << endm;
784  RError = 0.0 ;
785  dError = 0.0 ;
786  }
787  else if ( type.CompareTo("gauss", TString::kIgnoreCase) == 0 ){
788  LOG_INFO << "StFastGlauberMcMaker::InitRuRu Gaussian collision profile" << endm;
790  }
791  else if ( type.CompareTo("smallNpp", TString::kIgnoreCase) == 0 ){
792  mMode = 1 ;
793  LOG_INFO << "StFastGlauberMcMaker::InitRuRu Use small npp, large x for multiplicity (mMode = "
794  << mMode << ")"
795  << endm;
796  }
797  else if ( type.CompareTo("largeNpp", TString::kIgnoreCase) == 0 ){
798  mMode = 2 ;
799  LOG_INFO << "StFastGlauberMcMaker::InitRuRu Use large npp, small x for multiplicity (mMode = "
800  << mMode << ")"
801  << endm;
802  }
803 
804  return Init(A, R+RError, d+dError, beta2, beta4, A, R+RError, d+dError, beta2, beta4, type);
805 }
806 
807 
808 
809 
810 //____________________________________________________________________________________________________
811 Bool_t StFastGlauberMcMaker::IsCollision(Nucleon* nucleon0, Nucleon* nucleon1) const
812 {
813  const Double_t dx = nucleon0->GetX() - nucleon1->GetX() ;
814  const Double_t dy = nucleon0->GetY() - nucleon1->GetY() ;
815  const Double_t dt = TMath::Sqrt(dx*dx + dy*dy) ;
816  const Double_t cutoff = TMath::Sqrt(mInelasticNNCrossSection/TMath::Pi());
817 
818  switch ( mCollisionProfile ){
820  case mkHardCoreProfile:
821  {
822  return dt <= cutoff ;
823  }
824 
826  case mkGaussianProfile:
827  {
828  const Double_t sigma = cutoff ;
829  const Double_t arg = 0.5*dt/sigma ;
830  return ( StGlauberUtilities::instance()->GetUniform() <= TMath::Exp(-arg*arg) );
831  }
832 
833  default:
834  {
835  Error("StFastGlauberMcMaker::IsCollision", "No collision profile specified. see below");
836  LOG_INFO << " profile description function" << endm;
837  LOG_INFO << " Hard-core profile cut off distance = sqrt(sigma/pi) DoHardCoreCollision()" << endm;
838  LOG_INFO << " Gaussian profile sigma of gaussian = sqrt(sigma/pi) DoGaussianCollision()" << endm;
839  LOG_INFO << endm;
840  LOG_INFO << " d: transverse distrance between two nucleons (fm)" << endm;
841  LOG_INFO << " sigma: Inelastic N-N cross section (fm^2)" << endm;
842  assert(0);
843  }
844  }
845 
846  return kFALSE ;
847 }
848 
849 //____________________________________________________________________________________________________
851 {
853  Clear() ;
854 
855  StGlauberUtilities* glauberUtilities = StGlauberUtilities::instance() ;
856 
858  const Double_t impactParameter = glauberUtilities->GetImpactParameter() ;
859  mGlauberTree->SetB( impactParameter );
860 
861  //----------------------------------------------------------------------------------------------------
865  //----------------------------------------------------------------------------------------------------
866 
867  for(UInt_t in=0;in<2;in++){
868  // Rotation angle (theta, phi) for deformed nuclei
869  const Double_t Theta = (mIsDeformed[in]) ? glauberUtilities->GetTheta() : 0.0 ;
870  const Double_t Phi = (mIsDeformed[in]) ? glauberUtilities->GetPhi() : 0.0 ;
871  mGlauberTree->SetTheta(in, Theta);
872  mGlauberTree->SetPhi(in, Phi);
873 
875  const Double_t b = (in==0) ? -impactParameter/2.0 : impactParameter/2.0 ;
876 
877  UInt_t nNucleons = 0 ;
878  while( nNucleons < mNucleons[in].size() ){
879  Double_t r = 0.0 ;
880  Double_t theta = 0.0 ;
881  Double_t phi = 0.0 ;
882 
883  // Get (r,theta,phi) for nucleons
884  GetRThetaPhi(in, r, theta, phi) ;
885 
886  mhWoodsSaxon[in]->Fill(r, 1.0/(r*r)) ;
887 
888  if( mRepulsionDistance == 0.0 ){
890  mNucleons[in][nNucleons]->Set(r, theta, phi, b, Theta, Phi, kTRUE) ;
891  nNucleons++;
892  }
893  else{
894  if(nNucleons==0){
896  mNucleons[in][nNucleons]->Set(r, theta, phi, b, Theta, Phi, kTRUE) ;
897  nNucleons++;
898  }
899  else{
903  Bool_t isOk = kFALSE ;
904  do {
905  const Double_t x = r*TMath::Sin(theta)*TMath::Cos(phi) ;
906  const Double_t y = r*TMath::Sin(theta)*TMath::Sin(phi) ;
907  const Double_t z = r*TMath::Cos(theta);
908 
909  Bool_t isOverlap = kFALSE ;
910  for(UInt_t i=0; i<nNucleons; i++){
911  const Double_t x0 = mNucleons[in][i]->GetX() ;
912  const Double_t y0 = mNucleons[in][i]->GetY() ;
913  const Double_t z0 = mNucleons[in][i]->GetZ() ;
914  const Double_t dx = x - x0 ;
915  const Double_t dy = y - y0 ;
916  const Double_t dz = z - z0 ;
917 
918  // If any nucleons are found to be overlapped to the current one, stop the loop to save time
919  if( TMath::Sqrt(dx*dx + dy*dy + dz*dz) <= mRepulsionDistance ){
920  isOverlap = kTRUE ;
921  break ;
922  }
923  }// Looping over other nucleons
924 
925  if(isOverlap){
926  // There are overlap nucleons, try again
927  isOk = kFALSE ;
928 
929  // Get (r,theta,phi) for nucleons
930  GetRThetaPhi(in, r, theta, phi) ;
931  }
932  else{
933  // if no nucleons are overlapped with the current one
934  isOk = kTRUE ;
935  }
936 
937  } while( !isOk ) ; // Find any overlaped nucleons ?
938 
939  // Add positions and increment total
940  mNucleons[in][nNucleons]->Set(r, theta, phi, b, Theta, Phi, kTRUE) ;
941  nNucleons++;
942  }// Check overlap from 2nd nucleons
943  }// mRepulsionDistance > 0
944  }// Looping over all nucleons
945  }// Nucleus loop
946 
947  //----------------------------------------------------------------------------------------------------
949  //----------------------------------------------------------------------------------------------------
950  UInt_t Ncoll = 0 ;
951  for(UInt_t in=0; in<mNucleons[0].size(); in++){
952  Nucleon* nucleon0 = mNucleons[0][in] ;
953  for(UInt_t jn=0; jn<mNucleons[1].size(); jn++){
954  Nucleon* nucleon1 = mNucleons[1][jn] ;
955 
956  if(IsCollision(nucleon0, nucleon1)){
957  nucleon0->IncrementNcoll() ;
958  nucleon1->IncrementNcoll() ;
959 
960  Ncoll++;
961  }
962  }
963  }
964 
965  //----------------------------------------------------------------------------------------------------
967  //----------------------------------------------------------------------------------------------------
968  if( Ncoll == 0 ) return 0 ;
969 
970  //----------------------------------------------------------------------------------------------------
973  //----------------------------------------------------------------------------------------------------
974  Double_t nSum[4] ;
975  Double_t sumx[4] ;
976  Double_t sumy[4] ;
977  Double_t sumx2[4] ;
978  Double_t sumy2[4] ;
979  Double_t sumxy[4] ;
980 
981  for(Int_t i=0;i<4;i++){
982  nSum[i] = 0.0 ;
983  sumx[i] = 0.0 ;
984  sumy[i] = 0.0 ;
985  sumx2[i] = 0.0 ;
986  sumy2[i] = 0.0 ;
987  sumxy[i] = 0.0 ;
988  }
989 
990  UInt_t Npart = 0 ;
991  for(UInt_t in=0;in<2;in++){
992  for(UInt_t i=0;i<mNucleons[in].size();i++){
993  Nucleon* nucleon = mNucleons[in][i] ;
994  const UInt_t ncoll = nucleon->GetNcoll() ;
995  if( ncoll > 0 ){
996  // Participant weight
997  sumx[0] += nucleon->GetXYZ("x");
998  sumy[0] += nucleon->GetXYZ("y");
999  sumx2[0] += nucleon->GetXYZ("xx");
1000  sumy2[0] += nucleon->GetXYZ("yy");
1001  sumxy[0] += nucleon->GetXYZ("xy");
1002  nSum[0]++ ; // Should be identical to Npart
1003 
1004  // Ncoll weight
1005  sumx[1] += ncoll * nucleon->GetXYZ("x");
1006  sumy[1] += ncoll * nucleon->GetXYZ("y");
1007  sumx2[1] += ncoll * nucleon->GetXYZ("xx");
1008  sumy2[1] += ncoll * nucleon->GetXYZ("yy");
1009  sumxy[1] += ncoll * nucleon->GetXYZ("xy");
1010  nSum[1] += ncoll ;
1011 
1012  // Multiplicity weight
1013  const Double_t mult = mCentralityMaker->GetNegativeBinomial()->GetTwoComponentMultiplicity(1.0, ncoll) ;
1014  nucleon->SetMultiplicity(mult);
1015  sumx[2] += mult * nucleon->GetXYZ("x");
1016  sumy[2] += mult * nucleon->GetXYZ("y");
1017  sumx2[2] += mult * nucleon->GetXYZ("xx");
1018  sumy2[2] += mult * nucleon->GetXYZ("yy");
1019  sumxy[2] += mult * nucleon->GetXYZ("xy");
1020  nSum[2] += mult ;
1021 
1022  nucleon->IncrementNpart() ;
1023  Npart++;
1024  }
1025  else{
1026  // Spectator
1027  sumx[3] += nucleon->GetXYZ("x");
1028  sumy[3] += nucleon->GetXYZ("y");
1029  sumx2[3] += nucleon->GetXYZ("xx");
1030  sumy2[3] += nucleon->GetXYZ("yy");
1031  sumxy[3] += nucleon->GetXYZ("xy");
1032 
1033  nSum[3]++ ; // Should be identical to MassNumber - Npart
1034  }
1035  }// Nucleon loop
1036  }// Nucleus loop
1037 
1038  // Set multiplicity
1039  mGlauberTree->SetNpart(Npart) ;
1040  mGlauberTree->SetNcoll(Ncoll) ;
1041 
1042  // Get StNegativeBinomial with different (npp, x)
1043  const StNegativeBinomial* nbinomial = mCentralityMaker->GetNegativeBinomial() ;
1044  mGlauberTree->SetMultiplicity( nbinomial->GetMultiplicity(Npart, Ncoll) );
1045 
1046  // Get average and calculate eccentricity
1047  for(UInt_t i=0;i<4;i++){
1048  if(nSum[i]==0.0) continue ;
1049 
1050  sumx[i] /= nSum[i] ;
1051  sumy[i] /= nSum[i] ;
1052  sumx2[i] /= nSum[i] ;
1053  sumy2[i] /= nSum[i] ;
1054  sumxy[i] /= nSum[i] ;
1055 
1056  mGlauberTree->SetSumX(i, sumx[i]);
1057  mGlauberTree->SetSumY(i, sumy[i]);
1058  mGlauberTree->SetSumX2(i, sumx2[i]);
1059  mGlauberTree->SetSumY2(i, sumy2[i]);
1060  mGlauberTree->SetSumXY(i, sumxy[i]);
1061 
1062  // Eccentricity
1063  const Double_t sigmaX2 = sumx2[i] - sumx[i]*sumx[i] ;
1064  const Double_t sigmaY2 = sumy2[i] - sumy[i]*sumy[i] ;
1065  const Double_t sumSigma = sigmaX2 + sigmaY2 ;
1066  Double_t eccRP2 = -9999. ;
1067  Double_t eccPP2 = -9999. ;
1068  Double_t eccPP3 = -9999. ;
1069  Double_t eccPP4 = -9999. ;
1070  Double_t PP2 = -9999. ;
1071  Double_t PP3 = -9999. ;
1072  Double_t PP4 = -9999. ;
1073 
1074  if( sumSigma == 0.0 ){
1075  // Check denominator
1076  eccRP2 = -9999. ;
1077  eccPP2 = -9999. ;
1078  eccPP3 = -9999. ;
1079  eccPP4 = -9999. ;
1080  PP2 = -9999. ;
1081  PP3 = -9999. ;
1082  PP4 = -9999. ;
1083  }
1084  else{
1085  //----------------------------------------------------------------------------------------------------
1086  // Reaction plane eccentricity
1087  //----------------------------------------------------------------------------------------------------
1088  eccRP2 = (sigmaY2-sigmaX2)/sumSigma ;
1089 
1090  //----------------------------------------------------------------------------------------------------
1091  // The n-th order participant plane eccentricity
1092  //----------------------------------------------------------------------------------------------------
1093  for(Int_t io=0; io<3; io++){
1094  const Double_t order = io + 2.0 ;
1095  TGraph* qpart = GetParticipantEccentricity(order, sumx[i], sumy[i], nSum[i], i);
1096  if( io == 0 ) { PP2 = qpart->GetY()[2] ; eccPP2 = qpart->GetY()[3] ; }
1097  if( io == 1 ) { PP3 = qpart->GetY()[2] ; eccPP3 = qpart->GetY()[3] ; }
1098  if( io == 2 ) { PP4 = qpart->GetY()[2] ; eccPP4 = qpart->GetY()[3] ; }
1099  delete qpart ;
1100  }
1101  }
1102 
1103  mGlauberTree->SetEccRP2(i, eccRP2);
1104  mGlauberTree->SetEccPP2(i, eccPP2);
1105  mGlauberTree->SetEccPP3(i, eccPP3);
1106  mGlauberTree->SetEccPP4(i, eccPP4);
1107  mGlauberTree->SetPP2(i, PP2);
1108  mGlauberTree->SetPP3(i, PP3);
1109  mGlauberTree->SetPP4(i, PP4);
1110  }
1111 
1112  return 1;
1113 }
1114 
1115 //____________________________________________________________________________________________________
1116 Int_t StFastGlauberMcMaker::Run(const UInt_t nevents)
1117 {
1118  LOG_INFO << "StFastGlauberMcMaker::Run Run " << nevents << " events ..." << endm;
1119 
1120  while( mNeventsAccept < nevents ){
1121  if ( Make() == 1 ){
1122  // Debug
1123  if ( mNeventsAccept % 100 == 0 ){
1124  LOG_INFO << Form("StFastGlauberMcMaker::Run (accept/throw, Npart, Ncoll, b) = (%5d/%5d, %4d, %4d, %1.2f fm)",
1125  mNeventsAccept, mNeventsThrow, mGlauberTree->GetNpart(), mGlauberTree->GetNcoll(),
1126  mGlauberTree->GetB())
1127  << endm;
1128  }
1129 
1130  if(mDebug){
1131  LOG_INFO << Form("StFastGlauberMcMaker::Run (accept/throw, Npart, Ncoll, b) = (%5d/%5d, %4d, %4d, %1.2f fm)",
1132  mNeventsAccept, mNeventsThrow, mGlauberTree->GetNpart(), mGlauberTree->GetNcoll(),
1133  mGlauberTree->GetB())
1134  << endm;
1135  }
1136 
1137  // Fill tree
1138  mGlauberTree->Fill() ;
1139 
1140  mNeventsAccept++;
1141  }
1142 
1143  mNeventsThrow++;
1144  }
1145 
1146  return 1;
1147 }
1148 
1149 //____________________________________________________________________________________________________
1151 {
1152  // Store header info.
1153 
1154  // Name of nucleus
1155  const UInt_t A[2] = {mNucleons[0].size(), mNucleons[0].size()};
1156  mGlauberTree->SetNameNucleusA(GetName(A[0]));
1157  mGlauberTree->SetNameNucleusB(GetName(A[1]));
1158 
1159  // Mass numbers
1160  mGlauberTree->SetMassNumberA(A[0]);
1161  mGlauberTree->SetMassNumberB(A[1]);
1162 
1163  // Radius
1164  mGlauberTree->SetRadiusA((mIsDeformed[0]) ? mfWoodsSaxon2D[0]->GetParameter(0) : mfWoodsSaxon[0]->GetParameter(0));
1165  mGlauberTree->SetRadiusB((mIsDeformed[1]) ? mfWoodsSaxon2D[1]->GetParameter(0) : mfWoodsSaxon[1]->GetParameter(0));
1166 
1167  // Skin depth
1168  mGlauberTree->SetSkinDepthA((mIsDeformed[0]) ? mfWoodsSaxon2D[0]->GetParameter(1) : mfWoodsSaxon[0]->GetParameter(1));
1169  mGlauberTree->SetSkinDepthB((mIsDeformed[1]) ? mfWoodsSaxon2D[1]->GetParameter(1) : mfWoodsSaxon[1]->GetParameter(1));
1170 
1171  // Deformation parameters
1172  if(mIsDeformed[0]){
1173  mGlauberTree->SetBeta2A(mfWoodsSaxon2D[0]->GetParameter(2));
1174  mGlauberTree->SetBeta4A(mfWoodsSaxon2D[0]->GetParameter(3));
1175  }
1176  if(mIsDeformed[1]){
1177  mGlauberTree->SetBeta2A(mfWoodsSaxon2D[1]->GetParameter(2));
1178  mGlauberTree->SetBeta4A(mfWoodsSaxon2D[1]->GetParameter(3));
1179  }
1180 
1181  // sigmaNN, sqrt(sNN)
1182  mGlauberTree->SetSigmaNN( mInelasticNNCrossSection * 10.0 ) ; // mb
1183  mGlauberTree->SetSqrtSNN( mEnergy ) ;
1184 
1185  // Repulsion distance
1186  mGlauberTree->SetRepulsionD(mRepulsionDistance);
1187 
1188  // Total cross section (mb) skip if no events have been processed
1189  Double_t totalX = 0.0 ;
1190  Double_t totalXError = 0.0 ;
1191  const Double_t bmax = StGlauberUtilities::instance()->GetMaximumImpactParameter() ;
1192 
1193  if( mNeventsAccept != 0 ){
1194  const Double_t scaleFactor = TMath::Pi() * bmax * bmax * 10.0 ;
1195  totalX = (Double_t)mNeventsAccept / (Double_t)mNeventsThrow * scaleFactor ;
1196  totalXError = totalX * TMath::Sqrt(1.0/mNeventsAccept + 1.0/mNeventsThrow) ;
1197  }
1198  mGlauberTree->SetTotalXsec(totalX) ;
1199  mGlauberTree->SetTotalXsecError(totalXError) ;
1200 
1201  // Smearing options (0 or 1)
1202  // UInt_t smearHardCore = (mHardCoreSmearing) ? 1 : 0 ;
1203  // UInt_t smearGaussian = (mGaussianSmearing) ? 1 : 0 ;
1204  mGlauberTree->SetSmearHardCore(mHardCoreSmearing);
1205  mGlauberTree->SetSmearGaussian(mGaussianSmearing);
1206 
1207  // Collision profile options (0 or 1)
1208  UInt_t collisionHardCore = (mCollisionProfile==mkHardCoreProfile) ? 1 : 0 ;
1209  UInt_t collisionGaussian = (mCollisionProfile==mkGaussianProfile) ? 1 : 0 ;
1210  mGlauberTree->SetCollisionHardCore(collisionHardCore);
1211  mGlauberTree->SetCollisionGaussian(collisionGaussian);
1212 
1213  // Maximum impact parameter
1214  mGlauberTree->SetBMax(bmax);
1215 
1216  // Number of events
1217  mGlauberTree->SetNeventsAccept(mNeventsAccept);
1218  mGlauberTree->SetNeventsThrow(mNeventsThrow);
1219 
1220  // NBD parameters
1221  const StNegativeBinomial* nb = mCentralityMaker->GetNegativeBinomial() ;
1222  mGlauberTree->SetNpp( nb->GetNpp() );
1223  mGlauberTree->SetK( nb->GetK() );
1224  mGlauberTree->SetX( nb->GetX() );
1225  mGlauberTree->SetEfficiency( nb->GetEfficiency() );
1226  mGlauberTree->SetIsConstEfficiency( nb->IsConstEfficiency() );
1227 
1228  // version
1229  mGlauberTree->SetVersion(mVersion);
1230  mGlauberTree->FillHeader() ;
1231 
1232  LOG_INFO << "StFastGlauberMcMaker::Finish Close ROOT file" << endm;
1233  mGlauberTree->Close();
1234 
1235  return 1;
1236 }
1237 
1238 //____________________________________________________________________________________________________
1240 {
1241  mHardCoreSmearing = kTRUE ;
1242  LOG_INFO << "StFastGlauberMcMaker::DoHardCoreSmearing Hard-core smearing in (x,y,z) ";
1243  LOG_INFO << "by using the step function with inelastic nn cross section" << endm;
1244 
1246  mGaussianSmearing = kFALSE ;
1247  if(mDebug){
1248  LOG_INFO << "StFastGlauberMcMaker::DoHardCoreSmearing Gaussian smearing turned off" << endm;
1249  }
1250 }
1251 
1252 //____________________________________________________________________________________________________
1254 {
1255  mGaussianSmearing = kTRUE ;
1256  LOG_INFO << "StFastGlauberMcMaker::DoGaussianSmearing gaussian smearing in (x,y,z) ";
1257  LOG_INFO << "by using the gaussian function with fixed sigma" << endm;
1258 
1260  mHardCoreSmearing = kFALSE ;
1261  if(mDebug){
1262  LOG_INFO << "StFastGlauberMcMaker::DoGaussianSmearing hard-core smearing turned off" << endm;
1263  }
1264 }
1265 
1266 //____________________________________________________________________________________________________
1268 {
1270  mCollisionProfile = mkHardCoreProfile ;
1271  LOG_INFO << "StFastGlauberMcMaker::DoHardCoreCollision hard-core n-n collision "
1272  << " cut off distance is sqrt(sigma/pi) = " << TMath::Sqrt(mInelasticNNCrossSection/TMath::Pi())
1273  << endm;
1274 }
1275 
1276 //____________________________________________________________________________________________________
1278 {
1280  mCollisionProfile = mkGaussianProfile ;
1281  LOG_INFO << "StFastGlauberMcMaker::DoGaussianCollision gaussian profile n-n collision "
1282  << " sigma of gaussian is sqrt(sigma/pi) = " << TMath::Sqrt(mInelasticNNCrossSection/TMath::Pi())
1283  << endm;
1284 }
1285 
1286 //----------------------------------------------------------------------------------------------------
1287 // Utilities
1288 //----------------------------------------------------------------------------------------------------
1289 
1290 //____________________________________________________________________________________________________
1291 const Char_t* StFastGlauberMcMaker::GetName(const UInt_t massNumber) const
1292 {
1293  switch ( massNumber ){
1294  case 63: return "Cu";
1295  case 144: return "Sm";
1296  case 154: return "Sm";
1297  case 197: return "Au";
1298  case 208: return "Pb";
1299  case 238: return "U";
1300  //ZS The following should be changed to include Ru
1301  case 96: return "Zr";
1302  default:
1303  Error("StFastGlauberMcMaker::GetName", "can't find mass number = %d. Return white space", massNumber);
1304  return "";
1305  }
1306 
1307  return "";
1308 }
1309 
1310 //____________________________________________________________________________________________________
1311 Double_t StFastGlauberMcMaker::GetInelasticNNCrossSection(const Double_t energy, const TString type) const
1312 {
1314  const UInt_t energyInt = static_cast<UInt_t>(energy) ;
1315 
1317  Double_t error = 0.0 ; // absolute error on cross section
1318  if( type.CompareTo("smallXsec", TString::kIgnoreCase) == 0 ) error = -0.1 ;
1319  else if( type.CompareTo("largeXsec", TString::kIgnoreCase) == 0 ) error = 0.1 ;
1320 
1321  Double_t sigmaNN = 0.0 ;
1322  switch ( energyInt ) {
1323  case 2760: sigmaNN = 6.4 + error ; break ; // 64 mb
1324  case 200: sigmaNN = 4.2 + error ; break ; // 42 mb
1325  case 62: sigmaNN = 3.6 + error ; break ; // 36 mb
1326  case 39: sigmaNN = 3.4 + error ; break ; // 34 mb
1327  case 27: sigmaNN = 3.3 + error ; break ; // 33 mb
1328  case 19: sigmaNN = 3.2 + error ; break ; // 32 mb
1329  case 14: sigmaNN = 3.15 + error; break ; // 31.5 mb (for 14.5 GeV)
1330  case 11: sigmaNN = 3.12 + error; break ; // 31.2 mb (for 11.5 GeV)
1331  case 7: sigmaNN = 3.08 + error; break ; // 30.8 mb (for 7.7 GeV)
1332  default:
1333  Error("StFastGlauberMcMaker::GetInelasticNNCrossSection", "No energy = %1.1f GeV is available. abort", energy);
1334  assert(0);
1335  }
1336 
1337  if(mDebug){
1338  LOG_INFO << "StFastGlauberMcMaker::GetInelasticNNCrossSection ";
1339  LOG_INFO << "Inelastic NN cross section = " << sigmaNN * 10 << " (mb) for sqrt(sNN) = " << energy
1340  << endm;
1341  }
1342 
1343  return sigmaNN ;
1344 }
1345 
1346 //____________________________________________________________________________________________________
1347 void StFastGlauberMcMaker::GetRThetaPhi(const UInt_t inucleus, Double_t& r, Double_t& theta, Double_t& phi) const
1348 {
1349  if( mIsDeformed[inucleus] ){
1350  // Deformed nuclei
1351  Double_t cosTheta = -9999. ;
1352 
1353  // Random (r, cos(theta))
1354  mfWoodsSaxon2D[inucleus]->GetRandom2(r, cosTheta) ;
1355  theta = TMath::ACos(cosTheta) ; // convert into theta
1356  }
1357  else{
1358  // Spherical nuclei
1359  r = mfWoodsSaxon[inucleus]->GetRandom() ;
1360  theta = StGlauberUtilities::instance()->GetTheta() ; // cos(theta) profile in 0 < theta < pi
1361  }
1362  phi = StGlauberUtilities::instance()->GetPhi() ; // flat phi in -pi < phi < pi
1363 
1364  // Smearing if the switch is ON
1365  Smearing(r, theta, phi) ;
1366 }
1367 
1368 //____________________________________________________________________________________________________
1369 void StFastGlauberMcMaker::Smearing(Double_t& r, Double_t& theta, Double_t& phi) const
1370 {
1371  // Return original if both smearing options are OFF
1372  if(!mHardCoreSmearing && !mGaussianSmearing) return;
1373 
1374  // Hard-core smearing
1375  if(mHardCoreSmearing){
1376  Double_t dx, dy, dz ;
1377  mfHardCore->GetRandom3(dx, dy, dz) ;
1378 
1379  if(mDebug){
1380  LOG_INFO << Form("StFastGlauberMcMaker::Smearing Do hard-core smearing : (dx,dy,dz)=(%1.3f,%1.3f,%1.3f) fm",
1381  dx,dy,dz) << endm;
1382  }
1383 
1384  const Double_t x = r * TMath::Sin(theta) * TMath::Cos(phi) + dx ;
1385  const Double_t y = r * TMath::Sin(theta) * TMath::Sin(phi) + dy ;
1386  const Double_t z = r * TMath::Cos(theta) + dz ;
1387 
1388  r = TMath::Sqrt(x*x + y*y + z*z) ;
1389  theta = TMath::ACos(z/r) ;
1390  phi = TMath::ATan2(y, x) ;
1391  return;
1392  }
1393 
1394  // Gaussian smearing
1395  if(mGaussianSmearing){
1396  Double_t dx, dy, dz ;
1397  mfGaussian->GetRandom3(dx, dy, dz) ;
1398 
1399  if(mDebug){
1400  LOG_INFO << Form("StFastGlauberMcMaker::Smearing Do gaussian smearing : (dx,dy,dz)=(%1.3f,%1.3f,%1.3f) fm",
1401  dx,dy,dz) << endm;
1402  }
1403 
1404  const Double_t x = r * TMath::Sin(theta) * TMath::Cos(phi) + dx ;
1405  const Double_t y = r * TMath::Sin(theta) * TMath::Sin(phi) + dy ;
1406  const Double_t z = r * TMath::Cos(theta) + dz ;
1407 
1408  r = TMath::Sqrt(x*x + y*y + z*z) ;
1409  theta = TMath::ACos(z/r) ;
1410  phi = TMath::ATan2(y, x) ;
1411  return;
1412  }
1413 }
1414 
1415 //____________________________________________________________________________________________________
1416 TGraph* StFastGlauberMcMaker::GetParticipantEccentricity(const Double_t order, const Double_t sumx, const Double_t sumy,
1417  const Double_t sumw, const UInt_t weightId) const
1418 {
1419  Double_t qx = 0.0 ;
1420  Double_t qy = 0.0 ;
1421  Double_t qw = 0.0 ;
1422 
1423  for(UInt_t in=0;in<2;in++){
1424  for(UInt_t j=0;j<mNucleons[in].size();j++){
1425  Nucleon* nucleon = mNucleons[in][j] ;
1426  const UInt_t ncoll = nucleon->GetNcoll() ;
1427 
1428  const Bool_t isOk = (weightId==3 && ncoll==0) // Spectator
1429  || (weightId!=3 && ncoll > 0) ; // Participants
1430 
1431  if(!isOk) continue ;
1432 
1433  const Double_t x = nucleon->GetX() - sumx ;
1434  const Double_t y = nucleon->GetY() - sumy ;
1435  const Double_t z = nucleon->GetZ() ;
1436 
1437  const TVector3 vpart(x, y, z);
1438  const Double_t rt = vpart.Perp() ;
1439  const Double_t phi = vpart.Phi() ;
1440  Double_t weight = 1.0 ;
1441  if( weightId == 1 ) weight = ncoll ;
1442  if( weightId == 2 ) weight = nucleon->GetMultiplicity() ;
1443  if( weightId == 3 ) weight = 1.0 ;
1444 
1445  qx += -weight * rt*rt * TMath::Cos(order*phi) ; // -sign added to keep the same definition of participant plane
1446  qy += weight * rt*rt * TMath::Sin(order*phi) ;
1447  qw += weight * rt*rt ;
1448  }
1449  }
1450 
1451  qx /= sumw ;
1452  qy /= sumw ;
1453  qw /= sumw ;
1454 
1455  TGraph* g = new TGraph() ; // needs to be deleted later
1456  g->SetPoint(0, 0, qx);
1457  g->SetPoint(1, 1, qy);
1458  g->SetPoint(2, 2, TMath::ATan2(qy, qx));
1459  g->SetPoint(3, 3, TMath::Sqrt(qx*qx + qy*qy)/qw) ;
1460 
1461  return g ;
1462 }
1463 
1464 //____________________________________________________________________________________________________
1465 void StFastGlauberMcMaker::DebugOn()
1466 {
1467  mDebug = 1 ;
1468  LOG_INFO << "StFastGlauberMcMaker::DebugOn Set debug mode, will see a bunch of debugging messages" << endm;
1469 }
1470 
1471 //____________________________________________________________________________________________________
1472 void StFastGlauberMcMaker::Print(const TString option) const
1473 {
1474  if( option.CompareTo("type", TString::kIgnoreCase) == 0 ){
1475  LOG_INFO << endm;
1476  LOG_INFO << "#----------------------------------------------------------------------------------------------------" << endm;
1477  LOG_INFO << "StFastGlauberMcMaker::Print Current available types are:" << endm;
1478  LOG_INFO << "type description" << endm;
1479  LOG_INFO << " default" << endm;
1480  LOG_INFO << " large Large R(+2%), small d(-10%)" << endm;
1481  LOG_INFO << " small Small R(-2%), large d(+10%)" << endm;
1482  LOG_INFO << " largeXsec Large inelastic NN cross section (+1mb)" << endm;
1483  LOG_INFO << " smallXsec Small inelastic NN cross section (-1mb)" << endm;
1484  LOG_INFO << " gauss Use gaussian collision profile" << endm;
1485  LOG_INFO << " largeNpp Use large Npp, small x" << endm;
1486  LOG_INFO << " smallNpp Use small Npp, large x" << endm;
1487  LOG_INFO << "----------------------------------------------------------------------------------------------------" << endm;
1488  LOG_INFO << "NOTE: Types below are not relevant for generating trees." << endm;
1489  LOG_INFO << "NOTE: Only important for making histograms by StGlauberAnalysisMaker." << endm;
1490  LOG_INFO << "NOTE: Users should use the default outputs and modify type to see " << endm;
1491  LOG_INFO << "NOTE: the effect of different total cross section, for example" << endm;
1492  LOG_INFO << endm;
1493  LOG_INFO << " largeTotal +5% total cross section" << endm;
1494  LOG_INFO << " smallTotal -5% total cross section" << endm;
1495  LOG_INFO << " lowrw +2(-2) sigma p0 (p1) parameter for re-weighting correction" << endm;
1496  LOG_INFO << " highrw -2(+2) sigma p0 (p1) parameter for re-weighting correction" << endm;
1497  LOG_INFO << "#----------------------------------------------------------------------------------------------------" << endm;
1498  LOG_INFO << endm;
1499  }
1500 }
1501 
1502 //____________________________________________________________________________________________________
1504 {
1505  LOG_INFO << endm << endm;
1506  LOG_INFO << "StFastGlauberMcMaker::Version Current StFastGlauberMcMaker version is " << mVersion << endm;
1507  Print("type");
1508  LOG_INFO << endm << endm;
1509 
1510  return mVersion ;
1511 }
1512 
Int_t Clear()
Default destructor.
Int_t Finish()
Run Make() by nevents.
Definition: FJcore.h:367
void DoGaussianCollision()
Hard-core collision (default)
StFastGlauberMcMaker()
Current version.
void Sort()
Open ROOT file, initialize tree.
Double_t GetTheta() const
Get theta (polar angle)
Double_t GetMaximumImpactParameter() const
Get maximum impact parameter.
Int_t FillHeader()
Fill event-wise tree.
void Print(const TString option="") const
Gaussion profile collision.
void IncrementNcoll()
Increment Npart(x,y)
Definition: Nucleon.h:52
Double_t GetY() const
Get x position.
Definition: Nucleon.cxx:65
void SetRepulsionDistance(const Double_t repulsionDistance)
Default destructor.
void DoHardCoreCollision()
Default is OFF.
Double_t GetImpactParameter() const
Default destructor.
Int_t Close()
Fill header tree.
void SetMultiplicity(const Double_t val)
Increment Ncoll(x,y)
Definition: Nucleon.h:56
Double_t GetK() const
Get npp parameter.
void DoGaussianSmearing()
Default is OFF.
virtual ~StFastGlauberMcMaker()
Default constructor.
Double_t GetMultiplicity() const
Get Ncoll(x,y)
Definition: Nucleon.h:55
void SetB(const Double_t val)
Get entry for all branches.
void DoHardCoreSmearing()
Finish maker.
Double_t GetZ() const
Get y position.
Definition: Nucleon.cxx:71
Int_t Open(const TString filename)
Clear data members.
Bool_t IsConstEfficiency() const
Get mEfficiency (CAUTION: value has different meaning between constant and multiplicity dep...
Double_t GetXYZ(const Char_t *name="X") const
Get z position.
Definition: Nucleon.cxx:89
UInt_t GetNcoll() const
Get Npart(x,y)
Definition: Nucleon.h:54
Double_t GetEfficiency() const
Get x parameter.
UInt_t Version() const
Debug Mode ON.
const StNegativeBinomial * GetNegativeBinomial(const UInt_t id=0) const
Default destructor.
Definition: Nucleon.h:9
Int_t GetMultiplicity(const Double_t npart, const Double_t ncoll) const
Get multiplcity by convoluting NBD.
Int_t Fill()
Sort outputs in ROOT file.
Int_t Run(const UInt_t nevents)
Make one event.
Double_t GetPhi() const
Get phi (azimuthal angle)
Double_t GetX() const
Get k parameter.
Double_t GetTwoComponentMultiplicity(const Double_t npart, const Double_t ncoll) const
(1-x)*npart/2 + x*ncoll
void IncrementNpart()
Get weight factor to calculate average quantities.
Definition: Nucleon.h:51