Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ThermPtnSampler.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file ThermPtnSampler.cc
1 
2 #include "ThermPtnSampler.h"
3 #include "JetScapeXML.h"
4 #include "JetScapeLogger.h"
5 #include "JetScapeConstants.h"
6 #include <omp.h>
7 #include <vector>
8 #include <random>
9 #include <cmath>
10 #include <unordered_map>
11 #include <algorithm>
12 
13 using namespace Jetscape;
14 
15 ThermalPartonSampler::ThermalPartonSampler(unsigned int ran_seed){
16 
17  rng_engine.seed(ran_seed); // set the random seed from the framework
18 
19  // Gaussian weights for integration
20  GWeight[0] = 0.0621766166553473; GWeight[1] = 0.0619360674206832; GWeight[2] = 0.0614558995903167; GWeight[3] = 0.0607379708417702; GWeight[4] = 0.0597850587042655;
21  GWeight[5] = 0.0586008498132224; GWeight[6] = 0.0571899256477284; GWeight[7] = 0.0555577448062125; GWeight[8] = 0.0537106218889962; GWeight[9] = 0.0516557030695811;
22  GWeight[10] = 0.0494009384494663; GWeight[11] = 0.0469550513039484; GWeight[12] = 0.0443275043388033; GWeight[13] = 0.0415284630901477; GWeight[14] = 0.0385687566125877;
23  GWeight[15] = 0.0354598356151462; GWeight[16] = 0.0322137282235780; GWeight[17] = 0.0288429935805352; GWeight[18] = 0.0253606735700124; GWeight[19] = 0.0217802431701248;
24  GWeight[20] = 0.0181155607134894; GWeight[21] = 0.0143808227614856; GWeight[22] = 0.0105905483836510; GWeight[23] = 0.0067597991957454; GWeight[24] = 0.0029086225531551;
25  GWeight[25] = 0.0621766166553473; GWeight[26] = 0.0619360674206832; GWeight[27] = 0.0614558995903167; GWeight[28] = 0.0607379708417702; GWeight[29] = 0.0597850587042655;
26  GWeight[30] = 0.0586008498132224; GWeight[31] = 0.0571899256477284; GWeight[32] = 0.0555577448062125; GWeight[33] = 0.0537106218889962; GWeight[34] = 0.0516557030695811;
27  GWeight[35] = 0.0494009384494663; GWeight[36] = 0.0469550513039484; GWeight[37] = 0.0443275043388033; GWeight[38] = 0.0415284630901477; GWeight[39] = 0.0385687566125877;
28  GWeight[40] = 0.0354598356151462; GWeight[41] = 0.0322137282235780; GWeight[42] = 0.0288429935805352; GWeight[43] = 0.0253606735700124; GWeight[44] = 0.0217802431701248;
29  GWeight[45] = 0.0181155607134894; GWeight[46] = 0.0143808227614856; GWeight[47] = 0.0105905483836510; GWeight[48] = 0.0067597991957454; GWeight[49] = 0.0029086225531551;
30 
31  // Gaussian abscissas for integration
32  GAbs[0] = 0.0310983383271889; GAbs[1] = 0.0931747015600861; GAbs[2] = 0.1548905899981459; GAbs[3] = 0.2160072368760418; GAbs[4] = 0.2762881937795320;
33  GAbs[5] = 0.3355002454194373; GAbs[6] = 0.3934143118975651; GAbs[7] = 0.4498063349740388; GAbs[8] = 0.5044581449074642; GAbs[9] = 0.5571583045146501;
34  GAbs[10] = 0.6077029271849502; GAbs[11] = 0.6558964656854394; GAbs[12] = 0.7015524687068222; GAbs[13] = 0.7444943022260685; GAbs[14] = 0.7845558329003993;
35  GAbs[15] = 0.8215820708593360; GAbs[16] = 0.8554297694299461; GAbs[17] = 0.8859679795236131; GAbs[18] = 0.9130785566557919; GAbs[19] = 0.9366566189448780;
36  GAbs[20] = 0.9566109552428079; GAbs[21] = 0.9728643851066920; GAbs[22] = 0.9853540840480058; GAbs[23] = 0.9853540840480058; GAbs[24] = 0.9988664044200710;
37  GAbs[25] = -0.0310983383271889; GAbs[26] = -0.0931747015600861; GAbs[27] = -0.1548905899981459; GAbs[28] = -0.2160072368760418; GAbs[29] = -0.2762881937795320;
38  GAbs[30] = -0.3355002454194373; GAbs[31] = -0.3934143118975651; GAbs[32] = -0.4498063349740388; GAbs[33] = -0.5044581449074642; GAbs[34] = -0.5571583045146501;
39  GAbs[35] = -0.6077029271849502; GAbs[36] = -0.6558964656854394; GAbs[37] = -0.7015524687068222; GAbs[38] = -0.7444943022260685; GAbs[39] = -0.7845558329003993;
40  GAbs[40] = -0.8215820708593360; GAbs[41] = -0.8554297694299461; GAbs[42] = -0.8859679795236131; GAbs[43] = -0.9130785566557919; GAbs[44] = -0.9366566189448780;
41  GAbs[45] = -0.9566109552428079; GAbs[46] = -0.9728643851066920; GAbs[47] = -0.9853540840480058; GAbs[48] = -0.9853540840480058; GAbs[49] = -0.9988664044200710;
42 
43  // Adjustable parameters
44  xmq = 0.33/GEVFM; // use pythia values here
45  xms = 0.5/GEVFM; // use pythia values here
46  T = 0.165/GEVFM; // 165 MeV into fm^-1 // is set from outside with brick_temperature() or from hypersurface
47  NUMSTEP = 1048577; // 2^20+1, for steps of CDF Table, changes coarseness of momentum sampling
48 
49  // Adjustable params for 3+1d
50  CellDX = 0.2; CellDY = 0.2; CellDZ = 0.2; CellDT = 0.1; //these are the values used in SurfaceFinder.cc
51 
52  // Flags
53  SetNum = false; // Set 'true' to set number of particles by hand- !!!Statistics use above temperature!!!
54  SetNumLight = 1000000; //If SetNum == true, this many UD quarks are generated
55  SetNumStrange = 0 ; //If SetNum == true, this many S quarks are generated
56  ShuffleList = true; // Should list of particles be shuffled at the end
57 
58  // Brick Info
59  L = 4.0; // Thickness from box edge
60  W = 4.0; // x/y width of box
61  Time = 2.0; // Time of the brick partons
62 
63  Vx = 0.; // 'Uniform' no flow in x-dir
64  Vy = 0.; // 'Uniform' no flow in y-dir
65  Vz = 0.; // 'Uniform' no flow in z-dir
66 
67  surface.clear();
68 
69  num_ud = 0; num_s = 0;
70 
71  for(int iCDF=0; iCDF<NUMSTEP; ++iCDF){
72  std::vector<double> temp = {0., 0.};
73  CDFTabLight.push_back(temp);
74  CDFTabStrange.push_back(temp);
75  }
76 
77 }
78 
79 double ThermalPartonSampler::FermiPDF (double P, double M, double T, double mu){
80  return 1./( exp( (sqrt(M*M + P*P) - mu)/T ) + 1. );
81 }
82 
83 bool ThermalPartonSampler::SplitSort(double goal, int floor, int ceiling, int quark) {
84  std::vector<std::vector<double>>& CDFTab = (quark == 1) ? CDFTabLight : CDFTabStrange;
85 
86  int TargetPoint = ((floor + ceiling) / 2);
87  double TargetVal = CDFTab[TargetPoint][1];
88 
89  return (goal > TargetVal);
90 }
91 
92 void ThermalPartonSampler::CDFGenerator(double Temp, double M, int quark) {
93  double PStep;
94  double PMax = 10. * Temp; // CutOff for Integration
95  PStep = PMax / (NUMSTEP - 1); // Stepsize in P
96 
97  std::vector<std::vector<double>> CDFTab;
98  if (quark == 1) {
99  CDFTab = CDFTabLight;
100  } else if (quark == 2) {
101  CDFTab = CDFTabStrange;
102  } else {
103  JSWARN << "This is not a valid quark input for CDFGenerator()";
104  return;
105  }
106 
107  // Initialize tabulated results for CDF
108  CDFTab[0][0] = 0; // For zero momentum or less...
109  CDFTab[0][1] = 0; // There is zero chance
110 
111  // Precompute constant values used in the loop
112  double muPi0 = 0.0; // You need to provide the correct value of muPi0
113 
114  // Tabulate CDF(x) = int(0->x) PDF
115  for (int i = 1; i < NUMSTEP; i++) {
116  CDFTab[i][0] = CDFTab[i - 1][0] + PStep; // Calculate Momentum of the next step
117  double Fermi0 = FermiPDF(CDFTab[i - 1][0], M, Temp, muPi0); // PDF
118  double Fermi1 = FermiPDF(CDFTab[i][0], M, Temp, muPi0);
119  CDFTab[i][1] = CDFTab[i - 1][1] + (PStep / 2) * (Fermi0 * CDFTab[i - 1][0] * CDFTab[i - 1][0] + Fermi1 * CDFTab[i][0] * CDFTab[i][0]);
120  }
121 
122  // Normalize the CDF
123  for (int i = 0; i < NUMSTEP; i++) {
124  CDFTab[i][1] = CDFTab[i][1] / CDFTab[NUMSTEP - 1][1];
125  }
126 
127  if (quark == 1) {
128  CDFTabLight = CDFTab;
129  } else if (quark == 2) {
130  CDFTabStrange = CDFTab;
131  }
132 }
133 
134 void ThermalPartonSampler::MCSampler(double Temp, int quark) { // input temperature in fm^-1
135  double PMag, CosT, Phi, PStep;
136  double PMax = 10. * Temp; // CutOff for Integration
137  PStep = PMax / (NUMSTEP - 1); // Stepsize in P
138 
139  std::vector<std::vector<double>>& CDFTab = (quark == 1) ? CDFTabLight : CDFTabStrange;
140 
141  int floor = 0;
142  int ceiling = NUMSTEP - 1;
143  double PRoll;
144  double denominator;
145 
146  bool sample = true;
147  int test = 0;
148  while (sample) {
149  PRoll = ran();
150 
151  for (int i = 0; i < 25; i++) { // Use 25 iterations for both quark types
152  if (SplitSort(PRoll, floor, ceiling, quark)) {
153  floor = ((floor + ceiling) / 2);
154  } else {
155  ceiling = ((floor + ceiling) / 2);
156  }
157  }
158 
159  denominator = CDFTab[ceiling][1] - CDFTab[floor][1];
160  if (std::fabs(denominator) > 1e-16) {
161  PMag = PStep * (PRoll - CDFTab[floor][1]) / denominator + CDFTab[floor][0];
162  sample = false;
163  }
164  }
165 
166  CosT = (ran() - 0.5) * 2.0;
167  Phi = ran() * 2 * PI;
168 
169  NewX = PMag * sqrt(1 - CosT * CosT) * cos(Phi);
170  NewY = PMag * sqrt(1 - CosT * CosT) * sin(Phi);
171  NewZ = PMag * CosT;
172  NewP = PMag;
173 }
174 
175 void ThermalPartonSampler::samplebrick(){
176 
177  //preliminary parameter checks
178  if(L < 0.){L = -L; JSWARN << "Negative brick length - setting to positive " << L << " fm.";}
179  if(W < 0.){W = -W; JSWARN << "Negative brick width - setting to positive " << W << " fm.";}
180 
181  // Input read from cells
182  double LFSigma[4]; // LabFrame hypersurface (tau/t,x,y,eta/z=0)
183  double CMSigma[4]; // Center of mass hypersurface (tau/t,x,y,eta/z=0)
184  double Vel[4]; // Gamma & 3-velocity of cell (gamma, Vx, Vy, Vz) NOT FOUR VELOCITY
185 
186  // Calculated global quantities
187  int PartCount; // Total Count of Particles over ALL cells
188  double NumLight; // Number DENSITY of light quarks at set T
189  double NumStrange; // Number DENSITY of s quarks at set T
190  double UDDeg = 4.*6.; // Degeneracy of UD quarks
191  double OddDeg = 2.*6.; // Degeneracy of S quarks
192 
193  //counter for total number of light and strange quarks
194  int nL_tot = 0;
195  int nS_tot = 0;
196 
197  // Calculated quantities in cells
198  double LorBoost[4][4]; // Lorentz boost defined as used - form is always Lambda_u^v
199  int GeneratedParticles; // Number of particles to be generated this cell
200  double new_quark_energy; // store new quark energy for boost
201 
202  // End definition of static variables
203 
204  // Define hypersurface for brick here
205  // Default t=const hypersurface
206  // Vector must be covariant (negative signs on spatial components)
207  LFSigma[0] = 1.;
208  LFSigma[1] = 0.;
209  LFSigma[2] = 0.;
210  LFSigma[3] = 0.;
211 
212  Vel[1] = Vx;
213  Vel[2] = Vy;
214  Vel[3] = Vz;
215 
216  double vsquare = Vel[1]*Vel[1] + Vel[2]*Vel[2] + Vel[3]*Vel[3];
217 
218  if(vsquare > 1.){
219  JSWARN << "v^2 = " << vsquare;
220  JSWARN << "Unphysical velocity (brick flow)! Set to \"No Flow\" case";
221  Vel[1] = 0.;
222  Vel[2] = 0.;
223  Vel[3] = 0.;
224  }
225 
226  Vel[0] = 1. / std::sqrt(1.-vsquare); // gamma - Vel is not four velocity
227 
228  // Lambda_u ^v from (lab frame to rest frame with flow velocity)
229  if(vsquare == 0){
230  LorBoost[0][0]= Vel[0];
231  LorBoost[0][1]= Vel[0]*Vel[1];
232  LorBoost[0][2]= Vel[0]*Vel[2];
233  LorBoost[0][3]= Vel[0]*Vel[3];
234  LorBoost[1][0]= Vel[0]*Vel[1];
235  LorBoost[1][1]= 1.;
236  LorBoost[1][2]= 0.;
237  LorBoost[1][3]= 0.;
238  LorBoost[2][0]= Vel[0]*Vel[2];
239  LorBoost[2][1]= 0.;
240  LorBoost[2][2]= 1.;
241  LorBoost[2][3]= 0.;
242  LorBoost[3][0]= Vel[0]*Vel[3];
243  LorBoost[3][1]= 0.;
244  LorBoost[3][2]= 0.;
245  LorBoost[3][3]= 1.;
246  }else{
247  LorBoost[0][0]= Vel[0];
248  LorBoost[0][1]= Vel[0]*Vel[1];
249  LorBoost[0][2]= Vel[0]*Vel[2];
250  LorBoost[0][3]= Vel[0]*Vel[3];
251  LorBoost[1][0]= Vel[0]*Vel[1];
252  LorBoost[1][1]= (Vel[0] - 1.)*Vel[1]*Vel[1]/vsquare + 1.;
253  LorBoost[1][2]= (Vel[0] - 1.)*Vel[1]*Vel[2]/vsquare;
254  LorBoost[1][3]= (Vel[0] - 1.)*Vel[1]*Vel[3]/vsquare;
255  LorBoost[2][0]= Vel[0]*Vel[2];
256  LorBoost[2][1]= (Vel[0] - 1.)*Vel[1]*Vel[2]/vsquare;
257  LorBoost[2][2]= (Vel[0] - 1.)*Vel[2]*Vel[2]/vsquare + 1.;
258  LorBoost[2][3]= (Vel[0] - 1.)*Vel[2]*Vel[3]/vsquare;
259  LorBoost[3][0]= Vel[0]*Vel[3];
260  LorBoost[3][1]= (Vel[0] - 1.)*Vel[1]*Vel[3]/vsquare;
261  LorBoost[3][2]= (Vel[0] - 1.)*Vel[2]*Vel[3]/vsquare;
262  LorBoost[3][3]= (Vel[0] - 1.)*Vel[3]*Vel[3]/vsquare + 1.;
263  }
264 
265  // Code parity with hypersurface case
266  // Lambda_u^v Sigma_v = CMSigma_u
267  // Caution: CMSigma is a covariant vector
268  CMSigma[0] = (LorBoost[0][0]*LFSigma[0] + LorBoost[0][1]*LFSigma[1] + LorBoost[0][2]*LFSigma[2] + LorBoost[0][3]*LFSigma[3]);
269  CMSigma[1] = (LorBoost[1][0]*LFSigma[0] + LorBoost[1][1]*LFSigma[1] + LorBoost[1][2]*LFSigma[2] + LorBoost[1][3]*LFSigma[3]);
270  CMSigma[2] = (LorBoost[2][0]*LFSigma[0] + LorBoost[2][1]*LFSigma[1] + LorBoost[2][2]*LFSigma[2] + LorBoost[2][3]*LFSigma[3]);
271  CMSigma[3] = (LorBoost[3][0]*LFSigma[0] + LorBoost[3][1]*LFSigma[1] + LorBoost[3][2]*LFSigma[2] + LorBoost[3][3]*LFSigma[3]);
272 
273  /* Define Parton Densities */
274 
275  double cut = 10.*T; // Each coordinate of P is integrated this far
276  double E_light; // Store energy of light quarks
277  double E_strange; // Store energy of strange quarks
278  double GWeightProd; // Needed for integration
279  double pSpatialdSigma; // Needed for integration
280  PartCount = 0;
281  NumLight = 0; // Initialize density for light and strange quarks
282  NumStrange = 0;
283 
284  // Define the lambda function to compute the energy
285  auto computeEnergy = [](double mass, double cut, double GAbsL, double GAbsM, double GAbsK) {
286  return sqrt(mass * mass + (cut * GAbsL) * (cut * GAbsL) + (cut * GAbsM) * (cut * GAbsM) + (cut * GAbsK) * (cut * GAbsK));
287  };
288 
289  // Define the lambda function for the Fermi distribution
290  auto FermiDistribution = [](double energy, double temperature) {
291  return 1. / (exp(energy / temperature) + 1.);
292  };
293 
294  // GAUSSIAN INTEGRALS <n> = int f(p)d3p
295  for (int l=0; l<GPoints; l++){
296  for (int m=0; m<GPoints; m++){
297  for(int k=0; k<GPoints; k++){
298  GWeightProd = GWeight[l] * GWeight[m] * GWeight[k];
299  pSpatialdSigma = (cut * GAbs[l]) * CMSigma[1] + (cut * GAbs[m]) * CMSigma[2] + (cut * GAbs[k]) * CMSigma[3];
300 
301  // For UD quarks
302  E_light = computeEnergy(xmq, cut, GAbs[l], GAbs[m], GAbs[k]);
303  NumLight += GWeightProd * FermiDistribution(E_light, T) *
304  (E_light * CMSigma[0] + pSpatialdSigma) / E_light;
305 
306  // For S quarks
307  E_strange = computeEnergy(xms, cut, GAbs[l], GAbs[m], GAbs[k]);
308  NumStrange += GWeightProd * FermiDistribution(E_strange, T) *
309  (E_strange * CMSigma[0] + pSpatialdSigma) / E_strange;
310  }
311  }
312  }
313 
314  // Normalization factors: degeneracy, gaussian integration, and 1/(2Pi)^3
315  NumLight *= UDDeg*cut*cut*cut/(8.*PI*PI*PI);
316  NumStrange *= OddDeg*cut*cut*cut/(8.*PI*PI*PI);
317 
318  // Define rest-frame cumulative functions
319  CDFGenerator(T, xmq, 1); // for light quarks
320  CDFGenerator(T, xms, 2); // for strange quarks
321 
322  // U, D, UBAR, DBAR QUARKS
323  // <N> = V <n>
324  double NumHere = NumLight*L*W*W;
325 
326  // S, SBAR QUARKS
327  // <N> = V <n>
328  double NumOddHere = NumStrange*L*W*W;
329 
330  std::poisson_distribution<int> poisson_ud(NumHere);
331  std::poisson_distribution<int> poisson_s(NumOddHere);
332 
333  // In case of overwriting
334  if(SetNum){
335  NumHere = SetNumLight;
336  NumOddHere = SetNumStrange;
337  }
338 
339  // Generating light quarks
340  GeneratedParticles = poisson_ud(getRandomGenerator());
341 
342  // List of particles ( pos(x,y,z,t), mom(px,py,pz,E), species)
343  for(int partic = 0; partic < GeneratedParticles; partic++){
344  // adding space to PList for output quarks
345  std::vector<double> temp = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
346  Plist.push_back(temp);
347 
348  // Species - U,D,UBar,Dbar are equally likely
349  double SpecRoll = ran(); //Probability of species die roll
350  if (SpecRoll <= 0.25) { // UBar
351  Plist[PartCount][1] = -2;
352  } else if (SpecRoll <= 0.50) { // DBar
353  Plist[PartCount][1] = -1;
354  } else if (SpecRoll <= 0.75) { // D
355  Plist[PartCount][1] = 1;
356  } else { // U
357  Plist[PartCount][1] = 2;
358  }
359 
360  // Position
361  // Located at x,y pos of area element
362  double XRoll = ran()-0.5; // center at x=0
363  double YRoll = ran()-0.5; // center at y=0
364  double ZRoll = ran()-0.5; // center at z=0
365 
366  Plist[PartCount][7] = XRoll*L;
367  Plist[PartCount][8] = YRoll*W;
368  Plist[PartCount][9] = ZRoll*W;
369 
370  // Time
371  Plist[PartCount][10] = Time; // Tau = L/2.: assume jet at light speed
372 
373  // Momentum
374  // Sample rest frame momentum given T and mass of light quark
375  MCSampler(T, 1); // NewP=P, NewX=Px, ...
376 
377  // PLab^u = g^u^t Lambda_t ^v Pres^w g_w _v = Lambda ^u _w Pres_w (with velocity -v)
378  // (Lambda _u ^t with velocity v) == (Lambda ^u _t with velocity -v)
379  // This is boost as if particle sampled in rest frame
380  new_quark_energy = sqrt(xmq*xmq + NewP*NewP);
381  Plist[PartCount][6]= (LorBoost[0][0]*new_quark_energy + LorBoost[0][1]*NewX + LorBoost[0][2]*NewY + LorBoost[0][3]*NewZ)*GEVFM;
382  Plist[PartCount][3]= (LorBoost[1][0]*new_quark_energy + LorBoost[1][1]*NewX + LorBoost[1][2]*NewY + LorBoost[1][3]*NewZ)*GEVFM;
383  Plist[PartCount][4]= (LorBoost[2][0]*new_quark_energy + LorBoost[2][1]*NewX + LorBoost[2][2]*NewY + LorBoost[2][3]*NewZ)*GEVFM;
384  Plist[PartCount][5]= (LorBoost[3][0]*new_quark_energy + LorBoost[3][1]*NewX + LorBoost[3][2]*NewY + LorBoost[3][3]*NewZ)*GEVFM;
385  // Additional information
386  Plist[PartCount][0] = 1; // Event ID, to match jet formatting
387  Plist[PartCount][2] = 0; // Origin, to match jet formatting
388  Plist[PartCount][11] = 0; // Status - identifies as thermal quark
389  PartCount++;
390  }
391 
392  nL_tot += GeneratedParticles;
393 
394  // Generate strange quarks
395  GeneratedParticles = poisson_s(getRandomGenerator());
396 
397  nS_tot += GeneratedParticles;
398  //List of particles ( pos(x,y,z,t), mom(px,py,pz,E), species)
399  for(int partic = 0; partic < GeneratedParticles; partic++){
400  // adding space to PList for output quarks
401  std::vector<double> temp = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
402  Plist.push_back(temp);
403 
404  // Species - S,Sbar are equally likely
405  double SpecRoll = ran();
406  if(SpecRoll <= 0.5){ // SBar
407  Plist[PartCount][1] = -3;
408  } else { //S
409  Plist[PartCount][1] = 3;
410  }
411 
412  // Position
413  // Located at x,y pos of area element
414  double XRoll = ran()-0.5; // center at x=0
415  double YRoll = ran()-0.5; // center at y=0
416  double ZRoll = ran()-0.5; // center at z=0
417 
418  Plist[PartCount][7] = XRoll*L;
419  Plist[PartCount][8] = YRoll*W;
420  Plist[PartCount][9] = ZRoll*W;
421 
422  // Time
423  Plist[PartCount][10] = Time; // Tau = L/2.: assume jet at light speed
424 
425  // Momentum
426  // Sample rest frame momentum given T and mass of s quark
427  MCSampler(T, 2); // NewP=P, NewX=Px, ...
428 
429  // PLab^u = g^u^t Lambda_t ^v Pres^w g_w _v = Lambda ^u _w Pres_w (with velocity -v)
430  // (Lambda _u ^t with velocity v) == (Lambda ^u _t with velocity -v)
431  // This is boost as if particle sampled in rest frame
432  new_quark_energy = sqrt(xms*xms + NewP*NewP);
433  Plist[PartCount][6]= (LorBoost[0][0]*new_quark_energy + LorBoost[0][1]*NewX + LorBoost[0][2]*NewY + LorBoost[0][3]*NewZ)*GEVFM;
434  Plist[PartCount][3]= (LorBoost[1][0]*new_quark_energy + LorBoost[1][1]*NewX + LorBoost[1][2]*NewY + LorBoost[1][3]*NewZ)*GEVFM;
435  Plist[PartCount][4]= (LorBoost[2][0]*new_quark_energy + LorBoost[2][1]*NewX + LorBoost[2][2]*NewY + LorBoost[2][3]*NewZ)*GEVFM;
436  Plist[PartCount][5]= (LorBoost[3][0]*new_quark_energy + LorBoost[3][1]*NewX + LorBoost[3][2]*NewY + LorBoost[3][3]*NewZ)*GEVFM;
437 
438  // Additional information
439  Plist[PartCount][0] = 1; // Event ID, to match jet formatting
440  Plist[PartCount][2] = 0; // Origin, to match jet formatting
441  Plist[PartCount][11] = 0; // Status - identifies as thermal quark
442  PartCount++;
443  }
444 
445  JSDEBUG << "Light particles: " << nL_tot;
446  JSDEBUG << "Strange particles: " << nS_tot;
447 
448  num_ud = nL_tot;
449  num_s = nS_tot;
450 
451  //Shuffling PList
452  if(ShuffleList){
453  // Shuffle the Plist using the random engine
454  std::shuffle(&Plist[0], &Plist[PartCount], getRandomGenerator());
455  }
456 
457  //print Plist for testing
458  /*std::cout << std::setprecision(5);
459  std::ofstream thermalP;
460  thermalP.open("thermal_partons.dat", std::ios::app);
461  for(int p=0; p < Plist.size(); p++){
462  thermalP << Plist[p][0] << " " << Plist[p][1] << " " << Plist[p][2]
463  << " " << Plist[p][3] << " " << Plist[p][4] << " " << Plist[p][5]
464  << " " << Plist[p][6] << " " << Plist[p][7] << " " << Plist[p][8]
465  << " " << Plist[p][9] << " " << Plist[p][10] << " " << Plist[p][11]
466  << "\n";
467  }
468  thermalP.close();*/
469 }
470 
471 void ThermalPartonSampler::sample_3p1d(bool Cartesian_hydro){
472 
473  // Input read from cells
474  double CPos[4]; // Position of the current cell (tau/t, x, y , eta/z=0)
475  double LFSigma[4]; // LabFrame hypersurface (tau/t,x,y,eta/z=0), expect Sigma_mu
476  double CMSigma[4]; // Center of mass hypersurface (tau/t,x,y,eta/z=0)
477  double TRead; // Temperature of cell
478  double Vel[4]; // Gamma & 3-Velocity of cell (gamma, Vx, Vy, Vz) NOT FOUR VELOCITY
479  double tau_pos; // proper time from position of cell
480  double eta_pos; // eta from position of cell
481  double tau_sur; // proper time from normal vector of surface
482  double eta_sur; // eta from normal vector of surface
483 
484  // Calculated global quantities
485  int PartCount; // Total count of particles over ALL cells
486  double NumLight; // Number DENSITY of light quarks at set T
487  double NumStrange; // Number DENSITY of squarks at set T
488  double UDDeg = 4.*6.; // Degeneracy of UD quarks
489  double OddDeg = 2.*6.; // Degeneracy of S quarks
490 
491  // Calculated quantities in cells
492  double LorBoost[4][4]; // Lorentz boost defined as used - Form is always Lambda_u^v
493  int GeneratedParticles; // Number of particles to be generated this cell
494 
495  // Define parton densities
496  double cut; // Each coordinate of P is integrated this far
497  double E_light; // Store energy of light quarks
498  double E_strange; // Store energy of strange quarks
499  double GWeightProd; // Needed for integration
500  double pSpatialdSigma; // Needed for integration
501  PartCount = 0;
502  NumLight = 0; // Initialize density for light and strange quarks
503  NumStrange = 0;
504  GeneratedParticles = 0;
505  double new_quark_energy; // store new quark energy for boost
506 
507  //counter for total number of light and strange quarks
508  int nL_tot = 0;
509  int nS_tot = 0;
510 
511  // define the accuracy range for temperature values in which the CDFGenerator is executed (for the case the value is not in the cache yet)
512  const double accuracyRange = 0.003 / GEVFM; // for a usual hypersurface at const temperature there should be only one entry in the cache
513  // precompute CDFs and store them in a cache
514  std::unordered_map<double, std::vector<std::vector<double>>> cdfLightCache;
515  std::unordered_map<double, std::vector<std::vector<double>>> cdfStrangeCache;
516 
517  // lambda function to check if a temperature value is within the accuracy range of a cached temperature
518  auto withinAccuracyRange = [accuracyRange](double targetTemp, double cachedTemp) {
519  return std::fabs(targetTemp - cachedTemp) <= accuracyRange;
520  };
521 
522  // lambda function to retrieve the closest temperature from the cache within the accuracy range
523  auto getClosestCachedTemp = [](const std::unordered_map<double, std::vector<std::vector<double>>>& cache, double targetTemp) {
524  double closestTemp = -1.0;
525  double minDistance = std::numeric_limits<double>::max();
526  for (const auto& entry : cache) {
527  double cachedTemp = entry.first;
528  double distance = std::fabs(targetTemp - cachedTemp);
529  if (distance < minDistance) {
530  minDistance = distance;
531  closestTemp = cachedTemp;
532  }
533  }
534  return closestTemp;
535  };
536 
537  // Define the lambda function to compute the energy
538  auto computeEnergy = [](double mass, double cut, double GAbsL, double GAbsM, double GAbsK) {
539  return sqrt(mass * mass + (cut * GAbsL) * (cut * GAbsL) + (cut * GAbsM) * (cut * GAbsM) + (cut * GAbsK) * (cut * GAbsK));
540  };
541 
542  // Define the lambda function for the Fermi distribution
543  auto FermiDistribution = [](double energy, double temperature) {
544  return 1. / (exp(energy / temperature) + 1.);
545  };
546 
547  for(int iS=0; iS<surface.size(); ++iS){
548  tau_pos = surface[iS][0];
549  CPos[1] = surface[iS][1];
550  CPos[2] = surface[iS][2];
551  eta_pos = surface[iS][3]; //we need t,x,y,z
552  //this is also tau, x, y, eta
553  tau_sur = surface[iS][4];
554  LFSigma[1] = surface[iS][5];
555  LFSigma[2] = surface[iS][6];
556  eta_sur = surface[iS][7];
557  TRead = surface[iS][8] / GEVFM;
558  Vel[1] = surface[iS][9];
559  Vel[2] = surface[iS][10];
560  Vel[3] = surface[iS][11];
561 
562  cut = 10*TRead;
563 
564  // check if the CDFs for light quarks for this temperature are already in the cache and within the accuracy range
565  auto cdfLightIter = cdfLightCache.find(TRead);
566  if (cdfLightIter == cdfLightCache.end()) {
567  // if not found, find the closest temperature in the cache within the accuracy range
568  double closestTemp = getClosestCachedTemp(cdfLightCache,TRead);
569 
570  if (closestTemp >= 0. && withinAccuracyRange(TRead, closestTemp)) {
571  // use the CDFs from the closest temperature in the cache
572  CDFTabLight = cdfLightCache[closestTemp];
573  // set the current temperature value to the closest cached temperature
574  TRead = closestTemp;
575  } else {
576  // if no suitable entry is found, compute the CDFs and store them in the cache
577  CDFGenerator(TRead, xmq, 1); // for light quarks
578  cdfLightCache[TRead] = CDFTabLight;
579  }
580  } else {
581  // if found, use the precomputed CDFs from the cache
582  CDFTabLight = cdfLightIter->second;
583  }
584 
585  // check if the CDFs for strange quarks for this temperature are already in the cache and within the accuracy range
586  auto cdfStrangeIter = cdfStrangeCache.find(TRead);
587  if (cdfStrangeIter == cdfStrangeCache.end()) {
588  // if not found, find the closest temperature in the cache within the accuracy range
589  double closestTemp = getClosestCachedTemp(cdfStrangeCache,TRead);
590 
591  if (closestTemp >= 0 && withinAccuracyRange(TRead, closestTemp)) {
592  // use the CDFs from the closest temperature in the cache
593  CDFTabStrange = cdfStrangeCache[closestTemp];
594  // set the current temperature value to the closest cached temperature
595  TRead = closestTemp;
596  } else {
597  // if no suitable entry is found, compute the CDFs and store them in the cache
598  CDFGenerator(TRead, xms, 2); // for strange quarks
599  cdfStrangeCache[TRead] = CDFTabStrange;
600  }
601  } else {
602  // if found, use the precomputed CDFs from the cache
603  CDFTabStrange = cdfStrangeIter->second;
604  }
605 
606  if (Cartesian_hydro == false){
607  //getting t,z from tau, eta
608  CPos[0] = tau_pos*std::cosh(eta_pos);
609  CPos[3] = tau_pos*std::sinh(eta_pos);
610  //transform surface vector to Minkowski coordinates
611  double cosh_eta_pos = std::cosh(eta_pos);
612  double sinh_eta_pos = std::sinh(eta_pos);
613  LFSigma[0] = cosh_eta_pos*tau_sur - (sinh_eta_pos / tau_pos) * eta_sur;
614  LFSigma[3] = -sinh_eta_pos*tau_sur + (cosh_eta_pos / tau_pos) * eta_sur;
615  }else{ // check later!
616  CPos[0] = tau_pos;
617  CPos[3] = eta_pos;
618  LFSigma[0] = tau_sur;
619  LFSigma[3] = eta_sur;
620  }
621 
622  double vsquare = Vel[1]*Vel[1] + Vel[2]*Vel[2] + Vel[3]*Vel[3];
623  if(vsquare < 10e-16){
624  vsquare = 10e-16;
625  }
626 
627  // Deduced info
628  Vel[0] = 1. / sqrt(1-vsquare); // gamma - Vel is not four velocity
629 
630  // Lambda_u ^v
631  LorBoost[0][0]= Vel[0];
632  LorBoost[0][1]= Vel[0]*Vel[1];
633  LorBoost[0][2]= Vel[0]*Vel[2];
634  LorBoost[0][3]= Vel[0]*Vel[3];
635  LorBoost[1][0]= Vel[0]*Vel[1];
636  LorBoost[1][1]= (Vel[0] - 1.)*Vel[1]*Vel[1]/vsquare + 1.;
637  LorBoost[1][2]= (Vel[0] - 1.)*Vel[1]*Vel[2]/vsquare;
638  LorBoost[1][3]= (Vel[0] - 1.)*Vel[1]*Vel[3]/vsquare;
639  LorBoost[2][0]= Vel[0]*Vel[2];
640  LorBoost[2][1]= (Vel[0] - 1.)*Vel[1]*Vel[2]/vsquare;
641  LorBoost[2][2]= (Vel[0] - 1.)*Vel[2]*Vel[2]/vsquare + 1.;
642  LorBoost[2][3]= (Vel[0] - 1.)*Vel[2]*Vel[3]/vsquare;
643  LorBoost[3][0]= Vel[0]*Vel[3];
644  LorBoost[3][1]= (Vel[0] - 1.)*Vel[1]*Vel[3]/vsquare;
645  LorBoost[3][2]= (Vel[0] - 1.)*Vel[2]*Vel[3]/vsquare;
646  LorBoost[3][3]= (Vel[0] - 1.)*Vel[3]*Vel[3]/vsquare + 1.;
647 
648  if(vsquare == 0){
649  LorBoost[1][1]= 1.;
650  LorBoost[1][2]= 0;
651  LorBoost[1][3]= 0;
652  LorBoost[2][1]= 0;
653  LorBoost[2][2]= 1.;
654  LorBoost[2][3]= 0;
655  LorBoost[3][1]= 0;
656  LorBoost[3][2]= 0;
657  LorBoost[3][3]= 1.;
658  }
659  // Lambda_u^v Sigma_v = CMSigma_u
660  CMSigma[0] = (LorBoost[0][0]*LFSigma[0] + LorBoost[0][1]*LFSigma[1] + LorBoost[0][2]*LFSigma[2] + LorBoost[0][3]*LFSigma[3]);
661  CMSigma[1] = (LorBoost[1][0]*LFSigma[0] + LorBoost[1][1]*LFSigma[1] + LorBoost[1][2]*LFSigma[2] + LorBoost[1][3]*LFSigma[3]);
662  CMSigma[2] = (LorBoost[2][0]*LFSigma[0] + LorBoost[2][1]*LFSigma[1] + LorBoost[2][2]*LFSigma[2] + LorBoost[2][3]*LFSigma[3]);
663  CMSigma[3] = (LorBoost[3][0]*LFSigma[0] + LorBoost[3][1]*LFSigma[1] + LorBoost[3][2]*LFSigma[2] + LorBoost[3][3]*LFSigma[3]);
664 
665  // GAUSSIAN INTEGRALS <n> = int f(p)d3p
666  NumLight = 0.;
667  NumStrange = 0.;
668 
669  for (int l = 0; l < GPoints; l++) {
670  for (int m = 0; m < GPoints; m++) {
671  for (int k = 0; k < GPoints; k++) {
672  GWeightProd = GWeight[l] * GWeight[m] * GWeight[k];
673  pSpatialdSigma = (cut * GAbs[l]) * CMSigma[1] + (cut * GAbs[m]) * CMSigma[2] + (cut * GAbs[k]) * CMSigma[3];
674 
675  // For UD quarks
676  E_light = computeEnergy(xmq, cut, GAbs[l], GAbs[m], GAbs[k]);
677  NumLight += GWeightProd * FermiDistribution(E_light, TRead) *
678  (E_light * CMSigma[0] + pSpatialdSigma) / E_light;
679 
680  // For S quarks
681  E_strange = computeEnergy(xms, cut, GAbs[l], GAbs[m], GAbs[k]);
682  NumStrange += GWeightProd * FermiDistribution(E_strange, TRead) *
683  (E_strange * CMSigma[0] + pSpatialdSigma) / E_strange;
684  }
685  }
686  }
687 
688  // U, D, UBAR, DBAR QUARKS
689  // <N> = V <n>
690  double NumHere = NumLight*UDDeg*cut*cut*cut/(8.*PI*PI*PI);
691  std::poisson_distribution<int> poisson_ud(NumHere);
692 
693  // Generating light quarks
694  GeneratedParticles = poisson_ud(getRandomGenerator()); // Initialize particles created in this cell
695 
696  // List of particles ( pos(x,y,z,t), mom(px,py,pz,E), species)
697  for(int partic = 0; partic < GeneratedParticles; partic++){
698  // adding space to PList for output quarks
699  std::vector<double> temp = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
700  Plist.push_back(temp);
701 
702  // Species - U,D,UBar,Dbar are equally likely
703  double SpecRoll = ran(); //Probability of species die roll
704  if (SpecRoll <= 0.25) { // UBar
705  Plist[PartCount][1] = -2;
706  } else if (SpecRoll <= 0.50) { // DBar
707  Plist[PartCount][1] = -1;
708  } else if (SpecRoll <= 0.75) { // D
709  Plist[PartCount][1] = 1;
710  } else { // U
711  Plist[PartCount][1] = 2;
712  }
713 
714  // Position
715  // Located at x,y pos of area element
716  Plist[PartCount][10] = CPos[0] + (ran() - 0.5)*CellDT; // Tau
717  Plist[PartCount][7] = CPos[1] + (ran() - 0.5)*CellDX;
718  Plist[PartCount][8] = CPos[2] + (ran() - 0.5)*CellDY;
719  Plist[PartCount][9] = CPos[3] + (ran() - 0.5)*CellDZ;
720 
721  // Momentum
722  // Sample rest frame momentum given T and mass of light quark
723  MCSampler(TRead, 1); // NewP=P, NewX=Px, ...
724 
725  // USE THE SAME BOOST AS BEFORE
726  // PLab^u = g^u^t Lambda_t ^v pres^w g_w _v
727  // Returns P in GeV
728  new_quark_energy = sqrt(xmq*xmq + NewP*NewP);
729  Plist[PartCount][6] = (LorBoost[0][0]*new_quark_energy + LorBoost[0][1]*NewX + LorBoost[0][2]*NewY + LorBoost[0][3]*NewZ)*GEVFM;
730  Plist[PartCount][3] = (LorBoost[1][0]*new_quark_energy + LorBoost[1][1]*NewX + LorBoost[1][2]*NewY + LorBoost[1][3]*NewZ)*GEVFM;
731  Plist[PartCount][4] = (LorBoost[2][0]*new_quark_energy + LorBoost[2][1]*NewX + LorBoost[2][2]*NewY + LorBoost[2][3]*NewZ)*GEVFM;
732  Plist[PartCount][5] = (LorBoost[3][0]*new_quark_energy + LorBoost[3][1]*NewX + LorBoost[3][2]*NewY + LorBoost[3][3]*NewZ)*GEVFM;
733 
734  // Additional information
735  Plist[PartCount][0] = 1; // Event ID, to match jet formatting
736  Plist[PartCount][2] = 0; // Origin, to match jet formatting
737  Plist[PartCount][11] = 0; // Status - identifies as thermal quark
738  PartCount++;
739  }
740 
741  // S, SBAR QUARKS
742  // <N> = V <n>
743  double NumOddHere = NumStrange*OddDeg*cut*cut*cut/(8.*PI*PI*PI);
744  std::poisson_distribution<int> poisson_s(NumOddHere);
745 
746  // Generate s quarks
747  int nL = GeneratedParticles;
748  nL_tot += nL;
749 
750  GeneratedParticles = poisson_s(getRandomGenerator()); //Initialize particles created in this cell
751 
752  int nS = GeneratedParticles;
753  nS_tot += nS;
754  //List of particles ( pos(x,y,z,t), mom(px,py,pz,E), species)
755  for(int partic = 0; partic < GeneratedParticles; partic++){
756  // adding space to PList for output quarks
757  std::vector<double> temp = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
758  Plist.push_back(temp);
759 
760  // Species - S,Sbar are equally likely
761  double SpecRoll = ran();
762  if(SpecRoll <= 0.5){ // SBar
763  Plist[PartCount][1] = -3;
764  } else { //S
765  Plist[PartCount][1] = 3;
766  }
767 
768  // Position
769  // Located at x,y pos of area element
770  Plist[PartCount][10] = CPos[0] + (ran() - 0.5)*CellDT; // Tau
771  Plist[PartCount][7] = CPos[1] + (ran() - 0.5)*CellDX;
772  Plist[PartCount][8] = CPos[2] + (ran() - 0.5)*CellDY;
773  Plist[PartCount][9] = CPos[3] + (ran() - 0.5)*CellDZ;
774 
775  // Momentum
776  // Sample rest frame momentum given T and mass of s quark
777  MCSampler(TRead, 2); // NewP=P, NewX=Px, ...
778 
779  // USE THE SAME BOOST AS BEFORE
780  // PLab^u = g^u^t Lambda_t ^v pres^w g_w _v
781  // Returns P in GeV
782  new_quark_energy = sqrt(xms*xms + NewP*NewP);
783  Plist[PartCount][6] = (LorBoost[0][0]*new_quark_energy + LorBoost[0][1]*NewX + LorBoost[0][2]*NewY + LorBoost[0][3]*NewZ)*GEVFM;
784  Plist[PartCount][3] = (LorBoost[1][0]*new_quark_energy + LorBoost[1][1]*NewX + LorBoost[1][2]*NewY + LorBoost[1][3]*NewZ)*GEVFM;
785  Plist[PartCount][4] = (LorBoost[2][0]*new_quark_energy + LorBoost[2][1]*NewX + LorBoost[2][2]*NewY + LorBoost[2][3]*NewZ)*GEVFM;
786  Plist[PartCount][5] = (LorBoost[3][0]*new_quark_energy + LorBoost[3][1]*NewX + LorBoost[3][2]*NewY + LorBoost[3][3]*NewZ)*GEVFM;
787 
788  // Additional information
789  Plist[PartCount][0] = 1; // Event ID, to match jet formatting
790  Plist[PartCount][2] = 0; // Origin, to match jet formatting
791  Plist[PartCount][11] = 0; // Status - identifies as thermal quark
792  PartCount++;
793  }
794  }
795 
796  JSDEBUG << "Light particles: " << nL_tot;
797  JSDEBUG << "Strange particles: " << nS_tot;
798 
799  num_ud = nL_tot;
800  num_s = nS_tot;
801 
802  //Shuffling PList
803  if(ShuffleList){
804  // Shuffle the Plist using the random engine
805  std::shuffle(&Plist[0], &Plist[PartCount], getRandomGenerator());
806  }
807 
808  //print Plist for testing
809  /*std::cout << std::setprecision(5);
810  std::ofstream thermalP;
811  thermalP.open("thermal_partons.dat", std::ios::app);
812  for(int p=0; p < Plist.size(); p++){
813  thermalP << Plist[p][0] << " " << Plist[p][1] << " " << Plist[p][2]
814  << " " << Plist[p][3] << " " << Plist[p][4] << " " << Plist[p][5]
815  << " " << Plist[p][6] << " " << Plist[p][7] << " " << Plist[p][8]
816  << " " << Plist[p][9] << " " << Plist[p][10] << " " << Plist[p][11]
817  << "\n";
818  }
819  thermalP.close();*/
820 }
821 
822 void ThermalPartonSampler::sample_2p1d(double eta_max){
823 
824  // Input read from cells
825  double CPos[4]; // Position of the current cell (tau/t, x, y , eta/z=0)
826  double LFSigma[4]; // LabFrame hypersurface (tau/t,x,y,eta/z=0)
827  double CMSigma[4]; // Center of mass hypersurface (tau/t,x,y,eta/z=0)
828  double TRead; // Temperature of cell
829  double Vel[4]; // Gamma & 3-Velocity of cell (gamma, Vx, Vy, Vz) NOT FOUR VELOCITY
830  double tau_pos; // proper time from position of cell
831  double eta_pos; // eta from position of cell
832  double tau_sur; // proper time from normal vector of surface
833  double eta_sur; // eta from normal vector of surface
834 
835  // Calculated global quantities
836  int PartCount; // Total count of particles over ALL cells
837  double NumLight; // Number DENSITY of light quarks at set T
838  double NumStrange; // Number DENSITY of squarks at set T
839  double UDDeg = 4.*6.; // Degeneracy of UD quarks
840  double OddDeg = 2.*6.; // Degeneracy of S quarks
841 
842  // Calculated quantities in cells
843  double LorBoost[4][4]; // Lorentz boost defined as used - Form is always Lambda_u^v
844  int GeneratedParticles; // Number of particles to be generated this cell
845 
846  // Define parton densities
847  double cut; // Each coordinate of P is integrated this far
848  double E_light; // Store energy of light quarks
849  double E_strange; // Store energy of strange quarks
850  double GWeightProd; // Needed for integration
851  double pSpatialdSigma; // Needed for integration
852  PartCount = 0;
853  NumLight = 0; // Initialize density for light and strange quarks
854  NumStrange = 0;
855  GeneratedParticles = 0;
856  double new_quark_energy; // store new quark energy for boost
857 
858  //counter for total number of light and strange quarks
859  int nL_tot = 0;
860  int nS_tot = 0;
861 
862  // define the accuracy range for temperature values in which the CDFGenerator is executed (for the case the value is not in the cache yet)
863  const double accuracyRange = 0.003 / GEVFM; // for a usual hypersurface at const temperature there should be only one entry in the cache
864  // precompute CDFs and store them in a cache
865  std::unordered_map<double, std::vector<std::vector<double>>> cdfLightCache;
866  std::unordered_map<double, std::vector<std::vector<double>>> cdfStrangeCache;
867 
868  // lambda function to check if a temperature value is within the accuracy range of a cached temperature
869  auto withinAccuracyRange = [accuracyRange](double targetTemp, double cachedTemp) {
870  return std::fabs(targetTemp - cachedTemp) <= accuracyRange;
871  };
872 
873  // lambda function to retrieve the closest temperature from the cache within the accuracy range
874  auto getClosestCachedTemp = [](const std::unordered_map<double, std::vector<std::vector<double>>>& cache, double targetTemp) {
875  double closestTemp = -1.0;
876  double minDistance = std::numeric_limits<double>::max();
877  for (const auto& entry : cache) {
878  double cachedTemp = entry.first;
879  double distance = std::fabs(targetTemp - cachedTemp);
880  if (distance < minDistance) {
881  minDistance = distance;
882  closestTemp = cachedTemp;
883  }
884  }
885  return closestTemp;
886  };
887 
888  // Define the lambda function to compute the energy
889  auto computeEnergy = [](double mass, double cut, double GAbsL, double GAbsM, double GAbsK) {
890  return sqrt(mass * mass + (cut * GAbsL) * (cut * GAbsL) + (cut * GAbsM) * (cut * GAbsM) + (cut * GAbsK) * (cut * GAbsK));
891  };
892 
893  // Define the lambda function for the Fermi distribution
894  auto FermiDistribution = [](double energy, double temperature) {
895  return 1. / (exp(energy / temperature) + 1.);
896  };
897 
898  std::vector<double> NumLightList;
899  std::vector<double> NumStrangeList;
900 
901  double d_eta = CellDZ;
902 
903  int N_slices = std::floor((eta_max / d_eta)-0.5)+1;
904  for(int slice=1; slice <= (2*N_slices+1); slice++){
905  double eta_slice = (slice-N_slices-1)*d_eta;
906 
907  JSINFO << "Sampling thermal partons for pseudorapidity slice " << slice
908  << " (eta_min = " << eta_slice-(d_eta/2.) << ", eta_max = "
909  << eta_slice+(d_eta/2.) << ")";
910 
911  for(int iS=0; iS<surface.size(); ++iS){
912  tau_pos = surface[iS][0];
913  CPos[1] = surface[iS][1];
914  CPos[2] = surface[iS][2];
915  eta_pos = surface[iS][3]; //we need t,x,y,z
916  //this is also tau, x, y, eta
917  tau_sur = surface[iS][4];
918  LFSigma[1] = surface[iS][5];
919  LFSigma[2] = surface[iS][6];
920  eta_sur = surface[iS][7];
921  TRead = surface[iS][8] / GEVFM;
922  Vel[1] = surface[iS][9];
923  Vel[2] = surface[iS][10];
924  Vel[3] = surface[iS][11];
925 
926  cut = 10*TRead;
927 
928  // check if the CDFs for light quarks for this temperature are already in the cache and within the accuracy range
929  auto cdfLightIter = cdfLightCache.find(TRead);
930  if (cdfLightIter == cdfLightCache.end()) {
931  // if not found, find the closest temperature in the cache within the accuracy range
932  double closestTemp = getClosestCachedTemp(cdfLightCache,TRead);
933 
934  if (closestTemp >= 0. && withinAccuracyRange(TRead, closestTemp)) {
935  // use the CDFs from the closest temperature in the cache
936  CDFTabLight = cdfLightCache[closestTemp];
937  // set the current temperature value to the closest cached temperature
938  TRead = closestTemp;
939  } else {
940  // if no suitable entry is found, compute the CDFs and store them in the cache
941  CDFGenerator(TRead, xmq, 1); // for light quarks
942  cdfLightCache[TRead] = CDFTabLight;
943  }
944  } else {
945  // if found, use the precomputed CDFs from the cache
946  CDFTabLight = cdfLightIter->second;
947  }
948 
949  // check if the CDFs for strange quarks for this temperature are already in the cache and within the accuracy range
950  auto cdfStrangeIter = cdfStrangeCache.find(TRead);
951  if (cdfStrangeIter == cdfStrangeCache.end()) {
952  // if not found, find the closest temperature in the cache within the accuracy range
953  double closestTemp = getClosestCachedTemp(cdfStrangeCache,TRead);
954 
955  if (closestTemp >= 0 && withinAccuracyRange(TRead, closestTemp)) {
956  // use the CDFs from the closest temperature in the cache
957  CDFTabStrange = cdfStrangeCache[closestTemp];
958  // set the current temperature value to the closest cached temperature
959  TRead = closestTemp;
960  } else {
961  // if no suitable entry is found, compute the CDFs and store them in the cache
962  CDFGenerator(TRead, xms, 2); // for strange quarks
963  cdfStrangeCache[TRead] = CDFTabStrange;
964  }
965  } else {
966  // if found, use the precomputed CDFs from the cache
967  CDFTabStrange = cdfStrangeIter->second;
968  }
969 
970  //getting t,z from tau, eta
971  CPos[0] = tau_pos*std::cosh(eta_pos);
972  CPos[3] = tau_pos*std::sinh(eta_pos);
973  //transform surface vector to Minkowski coordinates
974  double cosh_eta_pos = std::cosh(eta_pos);
975  double sinh_eta_pos = std::sinh(eta_pos);
976  LFSigma[0] = cosh_eta_pos*tau_sur - (sinh_eta_pos / tau_pos) * eta_sur;
977  LFSigma[3] = -sinh_eta_pos*tau_sur + (cosh_eta_pos / tau_pos) * eta_sur;
978 
979  CellDZ = CPos[0] * 2. * std::sinh(d_eta/2.);
980 
981  double vsquare = Vel[1]*Vel[1] + Vel[2]*Vel[2] + Vel[3]*Vel[3];
982  if(vsquare < 10e-16){
983  vsquare = 10e-16;
984  }
985 
986  // Deduced info
987  Vel[0] = 1. / sqrt(1-vsquare); // gamma - Vel is not four velocity
988 
989  // Lambda_u ^v
990  LorBoost[0][0]= Vel[0];
991  LorBoost[0][1]= Vel[0]*Vel[1];
992  LorBoost[0][2]= Vel[0]*Vel[2];
993  LorBoost[0][3]= Vel[0]*Vel[3];
994  LorBoost[1][0]= Vel[0]*Vel[1];
995  LorBoost[1][1]= (Vel[0] - 1.)*Vel[1]*Vel[1]/vsquare + 1.;
996  LorBoost[1][2]= (Vel[0] - 1.)*Vel[1]*Vel[2]/vsquare;
997  LorBoost[1][3]= (Vel[0] - 1.)*Vel[1]*Vel[3]/vsquare;
998  LorBoost[2][0]= Vel[0]*Vel[2];
999  LorBoost[2][1]= (Vel[0] - 1.)*Vel[1]*Vel[2]/vsquare;
1000  LorBoost[2][2]= (Vel[0] - 1.)*Vel[2]*Vel[2]/vsquare + 1.;
1001  LorBoost[2][3]= (Vel[0] - 1.)*Vel[2]*Vel[3]/vsquare;
1002  LorBoost[3][0]= Vel[0]*Vel[3];
1003  LorBoost[3][1]= (Vel[0] - 1.)*Vel[1]*Vel[3]/vsquare;
1004  LorBoost[3][2]= (Vel[0] - 1.)*Vel[2]*Vel[3]/vsquare;
1005  LorBoost[3][3]= (Vel[0] - 1.)*Vel[3]*Vel[3]/vsquare + 1.;
1006 
1007  if(vsquare == 0){
1008  LorBoost[1][1]= 1.;
1009  LorBoost[1][2]= 0;
1010  LorBoost[1][3]= 0;
1011  LorBoost[2][1]= 0;
1012  LorBoost[2][2]= 1.;
1013  LorBoost[2][3]= 0;
1014  LorBoost[3][1]= 0;
1015  LorBoost[3][2]= 0;
1016  LorBoost[3][3]= 1.;
1017  }
1018 
1019  double NumHere;
1020  double NumOddHere;
1021  if(slice == 1){
1022  // Lambda_u^v Sigma_v = CMSigma_u
1023  CMSigma[0] = (LorBoost[0][0]*LFSigma[0] + LorBoost[0][1]*LFSigma[1] + LorBoost[0][2]*LFSigma[2] + LorBoost[0][3]*LFSigma[3]);
1024  CMSigma[1] = (LorBoost[1][0]*LFSigma[0] + LorBoost[1][1]*LFSigma[1] + LorBoost[1][2]*LFSigma[2] + LorBoost[1][3]*LFSigma[3]);
1025  CMSigma[2] = (LorBoost[2][0]*LFSigma[0] + LorBoost[2][1]*LFSigma[1] + LorBoost[2][2]*LFSigma[2] + LorBoost[2][3]*LFSigma[3]);
1026  CMSigma[3] = (LorBoost[3][0]*LFSigma[0] + LorBoost[3][1]*LFSigma[1] + LorBoost[3][2]*LFSigma[2] + LorBoost[3][3]*LFSigma[3]);
1027 
1028  // GAUSSIAN INTEGRALS <n> = int f(p)d3p
1029  NumLight = 0.;
1030  NumStrange = 0.;
1031 
1032  for (int l = 0; l < GPoints; l++) {
1033  for (int m = 0; m < GPoints; m++) {
1034  for (int k = 0; k < GPoints; k++) {
1035  GWeightProd = GWeight[l] * GWeight[m] * GWeight[k];
1036  pSpatialdSigma = (cut * GAbs[l]) * CMSigma[1] + (cut * GAbs[m]) * CMSigma[2] + (cut * GAbs[k]) * CMSigma[3];
1037 
1038  // For UD quarks
1039  E_light = computeEnergy(xmq, cut, GAbs[l], GAbs[m], GAbs[k]);
1040  NumLight += GWeightProd * FermiDistribution(E_light, TRead) *
1041  (E_light * CMSigma[0] + pSpatialdSigma) / E_light;
1042 
1043  // For S quarks
1044  E_strange = computeEnergy(xms, cut, GAbs[l], GAbs[m], GAbs[k]);
1045  NumStrange += GWeightProd * FermiDistribution(E_strange, TRead) *
1046  (E_strange * CMSigma[0] + pSpatialdSigma) / E_strange;
1047  }
1048  }
1049  }
1050  // U, D, UBAR, DBAR QUARKS
1051  // <N> = V <n>
1052  NumHere = NumLight*UDDeg*cut*cut*cut/(8.*PI*PI*PI);
1053 
1054  // S, SBAR QUARKS
1055  // <N> = V <n>
1056  NumOddHere = NumStrange*OddDeg*cut*cut*cut/(8.*PI*PI*PI);
1057 
1058  NumLightList.push_back(NumHere);
1059  NumStrangeList.push_back(NumOddHere);
1060  } else {
1061  NumHere = NumLightList[iS];
1062  NumOddHere = NumStrangeList[iS];
1063  }
1064 
1065  std::poisson_distribution<int> poisson_ud(NumHere);
1066 
1067  // Generating light quarks
1068  GeneratedParticles = poisson_ud(getRandomGenerator()); // Initialize particles created in this cell
1069 
1070  // List of particles ( pos(x,y,z,t), mom(px,py,pz,E), species)
1071  for(int partic = 0; partic < GeneratedParticles; partic++){
1072  // adding space to PList for output quarks
1073  std::vector<double> temp = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
1074  Plist.push_back(temp);
1075 
1076  // Species - U,D,UBar,Dbar are equally likely
1077  double SpecRoll = ran(); //Probability of species die roll
1078  if (SpecRoll <= 0.25) { // UBar
1079  Plist[PartCount][1] = -2;
1080  } else if (SpecRoll <= 0.50) { // DBar
1081  Plist[PartCount][1] = -1;
1082  } else if (SpecRoll <= 0.75) { // D
1083  Plist[PartCount][1] = 1;
1084  } else { // U
1085  Plist[PartCount][1] = 2;
1086  }
1087 
1088  // Position
1089  // Located at x,y pos of area element
1090  Plist[PartCount][10] = CPos[0] + (ran() - 0.5)*CellDT; // Tau
1091  Plist[PartCount][7] = CPos[1] + (ran() - 0.5)*CellDX;
1092  Plist[PartCount][8] = CPos[2] + (ran() - 0.5)*CellDY;
1093  Plist[PartCount][9] = CPos[3] + (ran() - 0.5)*CellDZ;
1094 
1095  if(std::abs(Plist[PartCount][9]) >= std::abs(Plist[PartCount][10])) {
1096  Plist[PartCount][10] = std::abs(Plist[PartCount][9]) + 10e-3;
1097  }
1098  double temp_t = std::sqrt(Plist[PartCount][10]*Plist[PartCount][10]-Plist[PartCount][9]*Plist[PartCount][9])
1099  * std::cosh(eta_slice+0.5*std::log((Plist[PartCount][10]+Plist[PartCount][9])/(Plist[PartCount][10]-Plist[PartCount][9])));
1100  double temp_z = std::sqrt(Plist[PartCount][10]*Plist[PartCount][10]-Plist[PartCount][9]*Plist[PartCount][9])
1101  * std::sinh(eta_slice+0.5*std::log((Plist[PartCount][10]+Plist[PartCount][9])/(Plist[PartCount][10]-Plist[PartCount][9])));
1102  Plist[PartCount][10] = temp_t;
1103  Plist[PartCount][9] = temp_z;
1104 
1105  // Momentum
1106  // Sample rest frame momentum given T and mass of light quark
1107  MCSampler(TRead, 1); // NewP=P, NewX=Px, ...
1108 
1109  Vel[3] = tanh(eta_slice);
1110  double vsquare = Vel[1]*Vel[1] + Vel[2]*Vel[2] + Vel[3]*Vel[3];
1111  if(vsquare < 10e-16){
1112  vsquare = 10e-16;
1113  }
1114 
1115  // Deduced info
1116  Vel[0] = 1. / sqrt(1-vsquare); // gamma - Vel is not four velocity
1117 
1118  // Lambda_u ^v
1119  LorBoost[0][0]= Vel[0];
1120  LorBoost[0][1]= Vel[0]*Vel[1];
1121  LorBoost[0][2]= Vel[0]*Vel[2];
1122  LorBoost[0][3]= Vel[0]*Vel[3];
1123  LorBoost[1][0]= Vel[0]*Vel[1];
1124  LorBoost[1][1]= (Vel[0] - 1.)*Vel[1]*Vel[1]/vsquare + 1.;
1125  LorBoost[1][2]= (Vel[0] - 1.)*Vel[1]*Vel[2]/vsquare;
1126  LorBoost[1][3]= (Vel[0] - 1.)*Vel[1]*Vel[3]/vsquare;
1127  LorBoost[2][0]= Vel[0]*Vel[2];
1128  LorBoost[2][1]= (Vel[0] - 1.)*Vel[1]*Vel[2]/vsquare;
1129  LorBoost[2][2]= (Vel[0] - 1.)*Vel[2]*Vel[2]/vsquare + 1.;
1130  LorBoost[2][3]= (Vel[0] - 1.)*Vel[2]*Vel[3]/vsquare;
1131  LorBoost[3][0]= Vel[0]*Vel[3];
1132  LorBoost[3][1]= (Vel[0] - 1.)*Vel[1]*Vel[3]/vsquare;
1133  LorBoost[3][2]= (Vel[0] - 1.)*Vel[2]*Vel[3]/vsquare;
1134  LorBoost[3][3]= (Vel[0] - 1.)*Vel[3]*Vel[3]/vsquare + 1.;
1135 
1136  if(vsquare == 0){
1137  LorBoost[1][1]= 1.;
1138  LorBoost[1][2]= 0;
1139  LorBoost[1][3]= 0;
1140  LorBoost[2][1]= 0;
1141  LorBoost[2][2]= 1.;
1142  LorBoost[2][3]= 0;
1143  LorBoost[3][1]= 0;
1144  LorBoost[3][2]= 0;
1145  LorBoost[3][3]= 1.;
1146  }
1147 
1148  // USE THE SAME BOOST AS BEFORE
1149  // PLab^u = g^u^t Lambda_t ^v pres^w g_w _v
1150  // Returns P in GeV
1151  new_quark_energy = sqrt(xmq*xmq + NewP*NewP);
1152  Plist[PartCount][6] = (LorBoost[0][0]*new_quark_energy + LorBoost[0][1]*NewX + LorBoost[0][2]*NewY + LorBoost[0][3]*NewZ)*GEVFM;
1153  Plist[PartCount][3] = (LorBoost[1][0]*new_quark_energy + LorBoost[1][1]*NewX + LorBoost[1][2]*NewY + LorBoost[1][3]*NewZ)*GEVFM;
1154  Plist[PartCount][4] = (LorBoost[2][0]*new_quark_energy + LorBoost[2][1]*NewX + LorBoost[2][2]*NewY + LorBoost[2][3]*NewZ)*GEVFM;
1155  Plist[PartCount][5] = (LorBoost[3][0]*new_quark_energy + LorBoost[3][1]*NewX + LorBoost[3][2]*NewY + LorBoost[3][3]*NewZ)*GEVFM;
1156 
1157  // Additional information
1158  Plist[PartCount][0] = 1; // Event ID, to match jet formatting
1159  Plist[PartCount][2] = 0; // Origin, to match jet formatting
1160  Plist[PartCount][11] = 0; // Status - identifies as thermal quark
1161  PartCount++;
1162  }
1163 
1164  std::poisson_distribution<int> poisson_s(NumOddHere);
1165 
1166  // Generate s quarks
1167  int nL = GeneratedParticles;
1168  nL_tot += nL;
1169 
1170  GeneratedParticles = poisson_s(getRandomGenerator()); //Initialize particles created in this cell
1171 
1172  int nS = GeneratedParticles;
1173  nS_tot += nS;
1174  //List of particles ( pos(x,y,z,t), mom(px,py,pz,E), species)
1175  for(int partic = 0; partic < GeneratedParticles; partic++){
1176  // adding space to PList for output quarks
1177  std::vector<double> temp = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
1178  Plist.push_back(temp);
1179 
1180  // Species - S,Sbar are equally likely
1181  double SpecRoll = ran();
1182  if(SpecRoll <= 0.5){ // SBar
1183  Plist[PartCount][1] = -3;
1184  } else { //S
1185  Plist[PartCount][1] = 3;
1186  }
1187 
1188  // Position
1189  // Located at x,y pos of area element
1190  Plist[PartCount][10] = CPos[0] + (ran() - 0.5)*CellDT; // Tau
1191  Plist[PartCount][7] = CPos[1] + (ran() - 0.5)*CellDX;
1192  Plist[PartCount][8] = CPos[2] + (ran() - 0.5)*CellDY;
1193  Plist[PartCount][9] = CPos[3] + (ran() - 0.5)*CellDZ;
1194 
1195  if(std::abs(Plist[PartCount][9]) >= std::abs(Plist[PartCount][10])) {
1196  Plist[PartCount][10] = std::abs(Plist[PartCount][9]) + 10e-3;
1197  }
1198  double temp_t = std::sqrt(Plist[PartCount][10]*Plist[PartCount][10]-Plist[PartCount][9]*Plist[PartCount][9])
1199  * std::cosh(eta_slice+0.5*std::log((Plist[PartCount][10]+Plist[PartCount][9])/(Plist[PartCount][10]-Plist[PartCount][9])));
1200  double temp_z = std::sqrt(Plist[PartCount][10]*Plist[PartCount][10]-Plist[PartCount][9]*Plist[PartCount][9])
1201  * std::sinh(eta_slice+0.5*std::log((Plist[PartCount][10]+Plist[PartCount][9])/(Plist[PartCount][10]-Plist[PartCount][9])));
1202  Plist[PartCount][10] = temp_t;
1203  Plist[PartCount][9] = temp_z;
1204 
1205  // Momentum
1206  // Sample rest frame momentum given T and mass of s quark
1207  MCSampler(TRead, 2); // NewP=P, NewX=Px, ...
1208 
1209  Vel[3] = tanh(eta_slice);
1210  double vsquare = Vel[1]*Vel[1] + Vel[2]*Vel[2] + Vel[3]*Vel[3];
1211  if(vsquare < 10e-16){
1212  vsquare = 10e-16;
1213  }
1214 
1215  // Deduced info
1216  Vel[0] = 1. / sqrt(1-vsquare); // gamma - Vel is not four velocity
1217 
1218  // Lambda_u ^v
1219  LorBoost[0][0]= Vel[0];
1220  LorBoost[0][1]= Vel[0]*Vel[1];
1221  LorBoost[0][2]= Vel[0]*Vel[2];
1222  LorBoost[0][3]= Vel[0]*Vel[3];
1223  LorBoost[1][0]= Vel[0]*Vel[1];
1224  LorBoost[1][1]= (Vel[0] - 1.)*Vel[1]*Vel[1]/vsquare + 1.;
1225  LorBoost[1][2]= (Vel[0] - 1.)*Vel[1]*Vel[2]/vsquare;
1226  LorBoost[1][3]= (Vel[0] - 1.)*Vel[1]*Vel[3]/vsquare;
1227  LorBoost[2][0]= Vel[0]*Vel[2];
1228  LorBoost[2][1]= (Vel[0] - 1.)*Vel[1]*Vel[2]/vsquare;
1229  LorBoost[2][2]= (Vel[0] - 1.)*Vel[2]*Vel[2]/vsquare + 1.;
1230  LorBoost[2][3]= (Vel[0] - 1.)*Vel[2]*Vel[3]/vsquare;
1231  LorBoost[3][0]= Vel[0]*Vel[3];
1232  LorBoost[3][1]= (Vel[0] - 1.)*Vel[1]*Vel[3]/vsquare;
1233  LorBoost[3][2]= (Vel[0] - 1.)*Vel[2]*Vel[3]/vsquare;
1234  LorBoost[3][3]= (Vel[0] - 1.)*Vel[3]*Vel[3]/vsquare + 1.;
1235 
1236  if(vsquare == 0){
1237  LorBoost[1][1]= 1.;
1238  LorBoost[1][2]= 0;
1239  LorBoost[1][3]= 0;
1240  LorBoost[2][1]= 0;
1241  LorBoost[2][2]= 1.;
1242  LorBoost[2][3]= 0;
1243  LorBoost[3][1]= 0;
1244  LorBoost[3][2]= 0;
1245  LorBoost[3][3]= 1.;
1246  }
1247 
1248  // USE THE SAME BOOST AS BEFORE
1249  // PLab^u = g^u^t Lambda_t ^v pres^w g_w _v
1250  // Returns P in GeV
1251  new_quark_energy = sqrt(xms*xms + NewP*NewP);
1252  Plist[PartCount][6] = (LorBoost[0][0]*new_quark_energy + LorBoost[0][1]*NewX + LorBoost[0][2]*NewY + LorBoost[0][3]*NewZ)*GEVFM;
1253  Plist[PartCount][3] = (LorBoost[1][0]*new_quark_energy + LorBoost[1][1]*NewX + LorBoost[1][2]*NewY + LorBoost[1][3]*NewZ)*GEVFM;
1254  Plist[PartCount][4] = (LorBoost[2][0]*new_quark_energy + LorBoost[2][1]*NewX + LorBoost[2][2]*NewY + LorBoost[2][3]*NewZ)*GEVFM;
1255  Plist[PartCount][5] = (LorBoost[3][0]*new_quark_energy + LorBoost[3][1]*NewX + LorBoost[3][2]*NewY + LorBoost[3][3]*NewZ)*GEVFM;
1256 
1257  // Additional information
1258  Plist[PartCount][0] = 1; // Event ID, to match jet formatting
1259  Plist[PartCount][2] = 0; // Origin, to match jet formatting
1260  Plist[PartCount][11] = 0; // Status - identifies as thermal quark
1261  PartCount++;
1262  }
1263  }
1264 
1265  }
1266 
1267  JSDEBUG << "Light particles: " << nL_tot;
1268  JSDEBUG << "Strange particles: " << nS_tot;
1269 
1270  num_ud = nL_tot;
1271  num_s = nS_tot;
1272 
1273  //Shuffling PList
1274  if(ShuffleList){
1275  // Shuffle the Plist using the random engine
1276  std::shuffle(&Plist[0], &Plist[PartCount], getRandomGenerator());
1277  }
1278 
1279  //print Plist for testing
1280  /*std::cout << std::setprecision(5);
1281  std::ofstream thermalP;
1282  thermalP.open("thermal_partons.dat", std::ios::app);
1283  for(int p=0; p < Plist.size(); p++){
1284  thermalP << Plist[p][0] << " " << Plist[p][1] << " " << Plist[p][2]
1285  << " " << Plist[p][3] << " " << Plist[p][4] << " " << Plist[p][5]
1286  << " " << Plist[p][6] << " " << Plist[p][7] << " " << Plist[p][8]
1287  << " " << Plist[p][9] << " " << Plist[p][10] << " " << Plist[p][11]
1288  << "\n";
1289  }
1290  thermalP.close();*/
1291 }