Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-03 08:20:02 

0001 
0002 #include "ThermPtnSampler.h"
0003 #include "JetScapeXML.h"
0004 #include "JetScapeLogger.h"
0005 #include "JetScapeConstants.h"
0006 #include <omp.h>
0007 #include <vector>
0008 #include <random>
0009 #include <cmath>
0010 #include <unordered_map>
0011 #include <algorithm>
0012 
0013 using namespace Jetscape;
0014 
0015 ThermalPartonSampler::ThermalPartonSampler(unsigned int ran_seed){
0016 
0017     rng_engine.seed(ran_seed); // set the random seed from the framework
0018 
0019     // Gaussian weights for integration
0020     GWeight[0]  = 0.0621766166553473; GWeight[1]  = 0.0619360674206832; GWeight[2]  = 0.0614558995903167; GWeight[3]  = 0.0607379708417702; GWeight[4]  = 0.0597850587042655;
0021     GWeight[5]  = 0.0586008498132224; GWeight[6]  = 0.0571899256477284; GWeight[7]  = 0.0555577448062125; GWeight[8]  = 0.0537106218889962; GWeight[9]  = 0.0516557030695811;
0022     GWeight[10] = 0.0494009384494663; GWeight[11] = 0.0469550513039484; GWeight[12] = 0.0443275043388033; GWeight[13] = 0.0415284630901477; GWeight[14] = 0.0385687566125877;
0023     GWeight[15] = 0.0354598356151462; GWeight[16] = 0.0322137282235780; GWeight[17] = 0.0288429935805352; GWeight[18] = 0.0253606735700124; GWeight[19] = 0.0217802431701248;
0024     GWeight[20] = 0.0181155607134894; GWeight[21] = 0.0143808227614856; GWeight[22] = 0.0105905483836510; GWeight[23] = 0.0067597991957454; GWeight[24] = 0.0029086225531551;
0025     GWeight[25] = 0.0621766166553473; GWeight[26] = 0.0619360674206832; GWeight[27] = 0.0614558995903167; GWeight[28] = 0.0607379708417702; GWeight[29] = 0.0597850587042655;
0026     GWeight[30] = 0.0586008498132224; GWeight[31] = 0.0571899256477284; GWeight[32] = 0.0555577448062125; GWeight[33] = 0.0537106218889962; GWeight[34] = 0.0516557030695811;
0027     GWeight[35] = 0.0494009384494663; GWeight[36] = 0.0469550513039484; GWeight[37] = 0.0443275043388033; GWeight[38] = 0.0415284630901477; GWeight[39] = 0.0385687566125877;
0028     GWeight[40] = 0.0354598356151462; GWeight[41] = 0.0322137282235780; GWeight[42] = 0.0288429935805352; GWeight[43] = 0.0253606735700124; GWeight[44] = 0.0217802431701248;
0029     GWeight[45] = 0.0181155607134894; GWeight[46] = 0.0143808227614856; GWeight[47] = 0.0105905483836510; GWeight[48] = 0.0067597991957454; GWeight[49] = 0.0029086225531551;
0030 
0031     // Gaussian abscissas for integration
0032     GAbs[0]  =  0.0310983383271889; GAbs[1]  =  0.0931747015600861; GAbs[2]  =  0.1548905899981459; GAbs[3]  =  0.2160072368760418; GAbs[4]  =  0.2762881937795320;
0033     GAbs[5]  =  0.3355002454194373; GAbs[6]  =  0.3934143118975651; GAbs[7]  =  0.4498063349740388; GAbs[8]  =  0.5044581449074642; GAbs[9]  =  0.5571583045146501;
0034     GAbs[10] =  0.6077029271849502; GAbs[11] =  0.6558964656854394; GAbs[12] =  0.7015524687068222; GAbs[13] =  0.7444943022260685; GAbs[14] =  0.7845558329003993;
0035     GAbs[15] =  0.8215820708593360; GAbs[16] =  0.8554297694299461; GAbs[17] =  0.8859679795236131; GAbs[18] =  0.9130785566557919; GAbs[19] =  0.9366566189448780;
0036     GAbs[20] =  0.9566109552428079; GAbs[21] =  0.9728643851066920; GAbs[22] =  0.9853540840480058; GAbs[23] =  0.9853540840480058; GAbs[24] =  0.9988664044200710;
0037     GAbs[25] = -0.0310983383271889; GAbs[26] = -0.0931747015600861; GAbs[27] = -0.1548905899981459; GAbs[28] = -0.2160072368760418; GAbs[29] = -0.2762881937795320;
0038     GAbs[30] = -0.3355002454194373; GAbs[31] = -0.3934143118975651; GAbs[32] = -0.4498063349740388; GAbs[33] = -0.5044581449074642; GAbs[34] = -0.5571583045146501;
0039     GAbs[35] = -0.6077029271849502; GAbs[36] = -0.6558964656854394; GAbs[37] = -0.7015524687068222; GAbs[38] = -0.7444943022260685; GAbs[39] = -0.7845558329003993;
0040     GAbs[40] = -0.8215820708593360; GAbs[41] = -0.8554297694299461; GAbs[42] = -0.8859679795236131; GAbs[43] = -0.9130785566557919; GAbs[44] = -0.9366566189448780;
0041     GAbs[45] = -0.9566109552428079; GAbs[46] = -0.9728643851066920; GAbs[47] = -0.9853540840480058; GAbs[48] = -0.9853540840480058; GAbs[49] = -0.9988664044200710;
0042 
0043     // Adjustable parameters
0044     xmq = 0.33/GEVFM; // use pythia values here
0045     xms = 0.5/GEVFM; // use pythia values here
0046     T = 0.165/GEVFM;  // 165 MeV into fm^-1 // is set from outside with brick_temperature() or from hypersurface
0047     NUMSTEP = 1048577;  // 2^20+1, for steps of CDF Table, changes coarseness of momentum sampling
0048 
0049     // Adjustable params for 3+1d
0050     CellDX = 0.2; CellDY = 0.2; CellDZ = 0.2; CellDT = 0.1; //these are the values used in SurfaceFinder.cc
0051 
0052     // Flags
0053     SetNum = false; // Set 'true' to set number of particles by hand- !!!Statistics use above temperature!!!
0054     SetNumLight = 1000000; //If SetNum == true, this many UD quarks are generated
0055     SetNumStrange = 0 ; //If SetNum == true, this many S quarks are generated
0056     ShuffleList = true; // Should list of particles be shuffled at the end
0057 
0058     // Brick Info
0059     L = 4.0; // Thickness from box edge
0060     W = 4.0; // x/y width of box
0061     Time = 2.0; // Time of the brick partons
0062 
0063     Vx = 0.; // 'Uniform' no flow in x-dir
0064     Vy = 0.; // 'Uniform' no flow in y-dir
0065     Vz = 0.; // 'Uniform' no flow in z-dir
0066 
0067     surface.clear();
0068 
0069     num_ud = 0; num_s = 0;
0070 
0071     for(int iCDF=0; iCDF<NUMSTEP; ++iCDF){
0072         std::vector<double> temp = {0., 0.};
0073         CDFTabLight.push_back(temp);
0074         CDFTabStrange.push_back(temp);
0075     }
0076 
0077 }
0078 
0079 double ThermalPartonSampler::FermiPDF (double P, double M, double T, double mu){
0080     return 1./( exp( (sqrt(M*M + P*P) - mu)/T ) + 1. );
0081 }
0082 
0083 bool ThermalPartonSampler::SplitSort(double goal, int floor, int ceiling, int quark) {
0084     std::vector<std::vector<double>>& CDFTab = (quark == 1) ? CDFTabLight : CDFTabStrange;
0085 
0086     int TargetPoint = ((floor + ceiling) / 2);
0087     double TargetVal = CDFTab[TargetPoint][1];
0088 
0089     return (goal > TargetVal);
0090 }
0091 
0092 void ThermalPartonSampler::CDFGenerator(double Temp, double M, int quark) {
0093     double PStep;
0094     double PMax = 10. * Temp;  // CutOff for Integration
0095     PStep = PMax / (NUMSTEP - 1); // Stepsize in P
0096 
0097     std::vector<std::vector<double>> CDFTab;
0098     if (quark == 1) {
0099         CDFTab = CDFTabLight;
0100     } else if (quark == 2) {
0101         CDFTab = CDFTabStrange;
0102     } else {
0103         JSWARN << "This is not a valid quark input for CDFGenerator()";
0104         return;
0105     }
0106 
0107     // Initialize tabulated results for CDF
0108     CDFTab[0][0] = 0; // For zero momentum or less...
0109     CDFTab[0][1] = 0; // There is zero chance
0110 
0111     // Precompute constant values used in the loop
0112     double muPi0 = 0.0; // You need to provide the correct value of muPi0
0113 
0114     // Tabulate CDF(x) = int(0->x) PDF
0115     for (int i = 1; i < NUMSTEP; i++) {
0116         CDFTab[i][0] = CDFTab[i - 1][0] + PStep; // Calculate Momentum of the next step
0117         double Fermi0 = FermiPDF(CDFTab[i - 1][0], M, Temp, muPi0); // PDF
0118         double Fermi1 = FermiPDF(CDFTab[i][0], M, Temp, muPi0);
0119         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]);
0120     }
0121 
0122     // Normalize the CDF
0123     for (int i = 0; i < NUMSTEP; i++) {
0124         CDFTab[i][1] = CDFTab[i][1] / CDFTab[NUMSTEP - 1][1];
0125     }
0126 
0127     if (quark == 1) {
0128         CDFTabLight = CDFTab;
0129     } else if (quark == 2) {
0130         CDFTabStrange = CDFTab;
0131     }
0132 }
0133 
0134 void ThermalPartonSampler::MCSampler(double Temp, int quark) { // input temperature in fm^-1
0135     double PMag, CosT, Phi, PStep;
0136     double PMax = 10. * Temp;  // CutOff for Integration
0137     PStep = PMax / (NUMSTEP - 1); // Stepsize in P
0138 
0139     std::vector<std::vector<double>>& CDFTab = (quark == 1) ? CDFTabLight : CDFTabStrange;
0140 
0141     int floor = 0;
0142     int ceiling = NUMSTEP - 1;
0143     double PRoll;
0144     double denominator;
0145 
0146     bool sample = true;
0147     int test = 0;
0148     while (sample) {
0149         PRoll = ran();
0150 
0151         for (int i = 0; i < 25; i++) { // Use 25 iterations for both quark types
0152             if (SplitSort(PRoll, floor, ceiling, quark)) {
0153                 floor = ((floor + ceiling) / 2);
0154             } else {
0155                 ceiling = ((floor + ceiling) / 2);
0156             }
0157         }
0158 
0159         denominator = CDFTab[ceiling][1] - CDFTab[floor][1];
0160         if (std::fabs(denominator) > 1e-16) {
0161             PMag = PStep * (PRoll - CDFTab[floor][1]) / denominator + CDFTab[floor][0];
0162             sample = false;
0163         }
0164     }
0165 
0166     CosT = (ran() - 0.5) * 2.0;
0167     Phi = ran() * 2 * PI;
0168 
0169     NewX = PMag * sqrt(1 - CosT * CosT) * cos(Phi);
0170     NewY = PMag * sqrt(1 - CosT * CosT) * sin(Phi);
0171     NewZ = PMag * CosT;
0172     NewP = PMag;
0173 }
0174 
0175 void ThermalPartonSampler::samplebrick(){
0176 
0177     //preliminary parameter checks
0178     if(L < 0.){L = -L; JSWARN << "Negative brick length - setting to positive " << L << " fm.";}
0179     if(W < 0.){W = -W; JSWARN << "Negative brick width - setting to positive " << W << " fm.";}
0180 
0181     // Input read from cells
0182     double LFSigma[4];      // LabFrame hypersurface (tau/t,x,y,eta/z=0)
0183     double CMSigma[4];      // Center of mass hypersurface (tau/t,x,y,eta/z=0)
0184     double Vel[4];          // Gamma & 3-velocity of cell (gamma, Vx, Vy, Vz) NOT FOUR VELOCITY
0185 
0186     // Calculated global quantities
0187     int PartCount;          // Total Count of Particles over ALL cells
0188     double NumLight;        // Number DENSITY of light quarks at set T
0189     double NumStrange;      // Number DENSITY of s quarks at set T
0190     double UDDeg = 4.*6.;   // Degeneracy of UD quarks
0191     double OddDeg = 2.*6.;  // Degeneracy of S quarks
0192 
0193     //counter for total number of light and strange quarks
0194     int nL_tot = 0;
0195     int nS_tot = 0;
0196 
0197     // Calculated quantities in cells
0198     double LorBoost[4][4];  // Lorentz boost defined as used - form is always Lambda_u^v
0199     int GeneratedParticles; // Number of particles to be generated this cell
0200     double new_quark_energy; // store new quark energy for boost
0201 
0202     // End definition of static variables
0203 
0204     // Define hypersurface for brick here
0205     // Default t=const hypersurface
0206     // Vector must be covariant (negative signs on spatial components)
0207     LFSigma[0] = 1.;
0208     LFSigma[1] = 0.;
0209     LFSigma[2] = 0.;
0210     LFSigma[3] = 0.;
0211 
0212     Vel[1] = Vx;
0213     Vel[2] = Vy;
0214     Vel[3] = Vz;
0215 
0216     double vsquare = Vel[1]*Vel[1] + Vel[2]*Vel[2] + Vel[3]*Vel[3];
0217 
0218     if(vsquare > 1.){
0219         JSWARN << "v^2 = " << vsquare;
0220         JSWARN << "Unphysical velocity (brick flow)! Set to \"No Flow\" case";
0221         Vel[1] = 0.;
0222         Vel[2] = 0.;
0223         Vel[3] = 0.;
0224     }
0225 
0226     Vel[0] = 1. / std::sqrt(1.-vsquare);   // gamma - Vel is not four velocity
0227 
0228     // Lambda_u ^v from (lab frame to rest frame with flow velocity)
0229     if(vsquare == 0){
0230         LorBoost[0][0]= Vel[0];
0231         LorBoost[0][1]= Vel[0]*Vel[1];
0232         LorBoost[0][2]= Vel[0]*Vel[2];
0233         LorBoost[0][3]= Vel[0]*Vel[3];
0234         LorBoost[1][0]= Vel[0]*Vel[1];
0235         LorBoost[1][1]= 1.;
0236         LorBoost[1][2]= 0.;
0237         LorBoost[1][3]= 0.;
0238         LorBoost[2][0]= Vel[0]*Vel[2];
0239         LorBoost[2][1]= 0.;
0240         LorBoost[2][2]= 1.;
0241         LorBoost[2][3]= 0.;
0242         LorBoost[3][0]= Vel[0]*Vel[3];
0243         LorBoost[3][1]= 0.;
0244         LorBoost[3][2]= 0.;
0245         LorBoost[3][3]= 1.;
0246     }else{
0247         LorBoost[0][0]= Vel[0];
0248         LorBoost[0][1]= Vel[0]*Vel[1];
0249         LorBoost[0][2]= Vel[0]*Vel[2];
0250         LorBoost[0][3]= Vel[0]*Vel[3];
0251         LorBoost[1][0]= Vel[0]*Vel[1];
0252         LorBoost[1][1]= (Vel[0] - 1.)*Vel[1]*Vel[1]/vsquare + 1.;
0253         LorBoost[1][2]= (Vel[0] - 1.)*Vel[1]*Vel[2]/vsquare;
0254         LorBoost[1][3]= (Vel[0] - 1.)*Vel[1]*Vel[3]/vsquare;
0255         LorBoost[2][0]= Vel[0]*Vel[2];
0256         LorBoost[2][1]= (Vel[0] - 1.)*Vel[1]*Vel[2]/vsquare;
0257         LorBoost[2][2]= (Vel[0] - 1.)*Vel[2]*Vel[2]/vsquare + 1.;
0258         LorBoost[2][3]= (Vel[0] - 1.)*Vel[2]*Vel[3]/vsquare;
0259         LorBoost[3][0]= Vel[0]*Vel[3];
0260         LorBoost[3][1]= (Vel[0] - 1.)*Vel[1]*Vel[3]/vsquare;
0261         LorBoost[3][2]= (Vel[0] - 1.)*Vel[2]*Vel[3]/vsquare;
0262         LorBoost[3][3]= (Vel[0] - 1.)*Vel[3]*Vel[3]/vsquare + 1.;
0263     }
0264 
0265     // Code parity with hypersurface case
0266     // Lambda_u^v Sigma_v = CMSigma_u
0267     // Caution: CMSigma is a covariant vector
0268     CMSigma[0] = (LorBoost[0][0]*LFSigma[0] + LorBoost[0][1]*LFSigma[1] + LorBoost[0][2]*LFSigma[2] + LorBoost[0][3]*LFSigma[3]);
0269     CMSigma[1] = (LorBoost[1][0]*LFSigma[0] + LorBoost[1][1]*LFSigma[1] + LorBoost[1][2]*LFSigma[2] + LorBoost[1][3]*LFSigma[3]);
0270     CMSigma[2] = (LorBoost[2][0]*LFSigma[0] + LorBoost[2][1]*LFSigma[1] + LorBoost[2][2]*LFSigma[2] + LorBoost[2][3]*LFSigma[3]);
0271     CMSigma[3] = (LorBoost[3][0]*LFSigma[0] + LorBoost[3][1]*LFSigma[1] + LorBoost[3][2]*LFSigma[2] + LorBoost[3][3]*LFSigma[3]);
0272 
0273     /* Define Parton Densities */
0274 
0275     double cut = 10.*T; // Each coordinate of P is integrated this far
0276     double E_light; // Store energy of light quarks
0277     double E_strange; // Store energy of strange quarks
0278     double GWeightProd; // Needed for integration
0279     double pSpatialdSigma; // Needed for integration
0280     PartCount = 0;
0281     NumLight = 0;       // Initialize density for light and strange quarks
0282     NumStrange = 0;
0283 
0284     // Define the lambda function to compute the energy
0285     auto computeEnergy = [](double mass, double cut, double GAbsL, double GAbsM, double GAbsK) {
0286         return sqrt(mass * mass + (cut * GAbsL) * (cut * GAbsL) + (cut * GAbsM) * (cut * GAbsM) + (cut * GAbsK) * (cut * GAbsK));
0287     };
0288 
0289     // Define the lambda function for the Fermi distribution
0290     auto FermiDistribution = [](double energy, double temperature) {
0291         return 1. / (exp(energy / temperature) + 1.);
0292     };
0293 
0294     // GAUSSIAN INTEGRALS <n> = int f(p)d3p
0295     for (int l=0; l<GPoints; l++){
0296         for (int m=0; m<GPoints; m++){
0297             for(int k=0; k<GPoints; k++){
0298                 GWeightProd = GWeight[l] * GWeight[m] * GWeight[k];
0299                 pSpatialdSigma = (cut * GAbs[l]) * CMSigma[1] + (cut * GAbs[m]) * CMSigma[2] + (cut * GAbs[k]) * CMSigma[3];
0300 
0301                 // For UD quarks
0302                 E_light = computeEnergy(xmq, cut, GAbs[l], GAbs[m], GAbs[k]);
0303                 NumLight += GWeightProd * FermiDistribution(E_light, T) *
0304                             (E_light * CMSigma[0] + pSpatialdSigma) / E_light;
0305 
0306                 // For S quarks
0307                 E_strange = computeEnergy(xms, cut, GAbs[l], GAbs[m], GAbs[k]);
0308                 NumStrange += GWeightProd * FermiDistribution(E_strange, T) *
0309                             (E_strange * CMSigma[0] + pSpatialdSigma) / E_strange;
0310             }
0311         }
0312     }
0313 
0314     // Normalization factors: degeneracy, gaussian integration, and 1/(2Pi)^3
0315     NumLight *= UDDeg*cut*cut*cut/(8.*PI*PI*PI);
0316     NumStrange *= OddDeg*cut*cut*cut/(8.*PI*PI*PI);
0317 
0318     // Define rest-frame cumulative functions
0319     CDFGenerator(T, xmq, 1); // for light quarks
0320     CDFGenerator(T, xms, 2); // for strange quarks
0321 
0322     // U, D, UBAR, DBAR QUARKS
0323     // <N> = V <n>
0324     double NumHere = NumLight*L*W*W;
0325 
0326     // S, SBAR QUARKS
0327     // <N> = V <n>
0328     double NumOddHere = NumStrange*L*W*W;
0329 
0330     std::poisson_distribution<int> poisson_ud(NumHere);
0331     std::poisson_distribution<int> poisson_s(NumOddHere);
0332 
0333     // In case of overwriting
0334     if(SetNum){
0335         NumHere = SetNumLight;
0336         NumOddHere = SetNumStrange;
0337     }
0338 
0339     // Generating light quarks
0340     GeneratedParticles = poisson_ud(getRandomGenerator());
0341 
0342     // List of particles ( pos(x,y,z,t), mom(px,py,pz,E), species)
0343     for(int partic = 0; partic < GeneratedParticles; partic++){
0344         // adding space to PList for output quarks
0345         std::vector<double> temp = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
0346         Plist.push_back(temp);
0347 
0348         // Species - U,D,UBar,Dbar are equally likely
0349         double SpecRoll = ran(); //Probability of species die roll
0350         if (SpecRoll <= 0.25) { // UBar
0351             Plist[PartCount][1] = -2;
0352         } else if (SpecRoll <= 0.50) { // DBar
0353             Plist[PartCount][1] = -1;
0354         } else if (SpecRoll <= 0.75) { // D
0355             Plist[PartCount][1] = 1;
0356         } else { // U
0357             Plist[PartCount][1] = 2;
0358         }
0359 
0360         // Position
0361         // Located at x,y pos of area element
0362         double XRoll = ran()-0.5; // center at x=0
0363         double YRoll = ran()-0.5; // center at y=0
0364         double ZRoll = ran()-0.5; // center at z=0
0365 
0366         Plist[PartCount][7] = XRoll*L;
0367         Plist[PartCount][8] = YRoll*W;
0368         Plist[PartCount][9] = ZRoll*W;
0369 
0370         // Time
0371         Plist[PartCount][10] = Time; // Tau = L/2.: assume jet at light speed
0372 
0373         // Momentum
0374         // Sample rest frame momentum given T and mass of light quark
0375         MCSampler(T, 1); // NewP=P, NewX=Px, ...
0376 
0377         // PLab^u = g^u^t Lambda_t ^v Pres^w g_w _v  = Lambda ^u _w Pres_w (with velocity -v)
0378         // (Lambda _u ^t with velocity v) == (Lambda ^u _t with velocity -v)
0379         // This is boost as if particle sampled in rest frame
0380         new_quark_energy = sqrt(xmq*xmq + NewP*NewP);
0381         Plist[PartCount][6]= (LorBoost[0][0]*new_quark_energy + LorBoost[0][1]*NewX + LorBoost[0][2]*NewY + LorBoost[0][3]*NewZ)*GEVFM;
0382         Plist[PartCount][3]= (LorBoost[1][0]*new_quark_energy + LorBoost[1][1]*NewX + LorBoost[1][2]*NewY + LorBoost[1][3]*NewZ)*GEVFM;
0383         Plist[PartCount][4]= (LorBoost[2][0]*new_quark_energy + LorBoost[2][1]*NewX + LorBoost[2][2]*NewY + LorBoost[2][3]*NewZ)*GEVFM;
0384         Plist[PartCount][5]= (LorBoost[3][0]*new_quark_energy + LorBoost[3][1]*NewX + LorBoost[3][2]*NewY + LorBoost[3][3]*NewZ)*GEVFM;
0385         // Additional information
0386         Plist[PartCount][0]  = 1; // Event ID, to match jet formatting
0387         Plist[PartCount][2]  = 0; // Origin, to match jet formatting
0388         Plist[PartCount][11] = 0; // Status - identifies as thermal quark
0389         PartCount++;
0390     }
0391 
0392     nL_tot += GeneratedParticles;
0393 
0394     // Generate strange quarks
0395     GeneratedParticles = poisson_s(getRandomGenerator());
0396 
0397     nS_tot += GeneratedParticles;
0398     //List of particles ( pos(x,y,z,t), mom(px,py,pz,E), species)
0399     for(int partic = 0; partic < GeneratedParticles; partic++){
0400         // adding space to PList for output quarks
0401         std::vector<double> temp = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
0402         Plist.push_back(temp);
0403 
0404         // Species - S,Sbar are equally likely
0405         double SpecRoll = ran();
0406         if(SpecRoll <= 0.5){ // SBar
0407             Plist[PartCount][1] = -3;
0408         } else { //S
0409             Plist[PartCount][1] = 3;
0410         }
0411 
0412         // Position
0413         // Located at x,y pos of area element
0414         double XRoll = ran()-0.5; // center at x=0
0415         double YRoll = ran()-0.5; // center at y=0
0416         double ZRoll = ran()-0.5; // center at z=0
0417 
0418         Plist[PartCount][7] = XRoll*L;
0419         Plist[PartCount][8] = YRoll*W;
0420         Plist[PartCount][9] = ZRoll*W;
0421 
0422         // Time
0423         Plist[PartCount][10] = Time; // Tau = L/2.: assume jet at light speed
0424 
0425         // Momentum
0426         // Sample rest frame momentum given T and mass of s quark
0427         MCSampler(T, 2); // NewP=P, NewX=Px, ...
0428 
0429         // PLab^u = g^u^t Lambda_t ^v Pres^w g_w _v  = Lambda ^u _w Pres_w (with velocity -v)
0430         // (Lambda _u ^t with velocity v) == (Lambda ^u _t with velocity -v)
0431         // This is boost as if particle sampled in rest frame
0432         new_quark_energy = sqrt(xms*xms + NewP*NewP);
0433         Plist[PartCount][6]= (LorBoost[0][0]*new_quark_energy + LorBoost[0][1]*NewX + LorBoost[0][2]*NewY + LorBoost[0][3]*NewZ)*GEVFM;
0434         Plist[PartCount][3]= (LorBoost[1][0]*new_quark_energy + LorBoost[1][1]*NewX + LorBoost[1][2]*NewY + LorBoost[1][3]*NewZ)*GEVFM;
0435         Plist[PartCount][4]= (LorBoost[2][0]*new_quark_energy + LorBoost[2][1]*NewX + LorBoost[2][2]*NewY + LorBoost[2][3]*NewZ)*GEVFM;
0436         Plist[PartCount][5]= (LorBoost[3][0]*new_quark_energy + LorBoost[3][1]*NewX + LorBoost[3][2]*NewY + LorBoost[3][3]*NewZ)*GEVFM;
0437 
0438         // Additional information
0439         Plist[PartCount][0]  = 1; // Event ID, to match jet formatting
0440         Plist[PartCount][2]  = 0; // Origin, to match jet formatting
0441         Plist[PartCount][11] = 0; // Status - identifies as thermal quark
0442         PartCount++;
0443     }
0444 
0445     JSDEBUG << "Light particles: " << nL_tot;
0446     JSDEBUG << "Strange particles: " << nS_tot;
0447 
0448     num_ud = nL_tot;
0449     num_s = nS_tot;
0450 
0451     //Shuffling PList
0452     if(ShuffleList){
0453         // Shuffle the Plist using the random engine
0454         std::shuffle(&Plist[0], &Plist[PartCount], getRandomGenerator());
0455     }
0456 
0457     //print Plist for testing
0458     /*std::cout << std::setprecision(5);
0459     std::ofstream thermalP;
0460     thermalP.open("thermal_partons.dat", std::ios::app);
0461     for(int p=0; p < Plist.size(); p++){
0462         thermalP << Plist[p][0] << " " << Plist[p][1] << " " << Plist[p][2]
0463         << " " << Plist[p][3] << " " << Plist[p][4] << " " << Plist[p][5]
0464         << " " << Plist[p][6] << " " << Plist[p][7] << " " << Plist[p][8]
0465         << " " << Plist[p][9] << " " << Plist[p][10] << " " << Plist[p][11]
0466         << "\n";
0467     }
0468     thermalP.close();*/
0469 }
0470 
0471 void ThermalPartonSampler::sample_3p1d(bool Cartesian_hydro){
0472 
0473     // Input read from cells
0474     double CPos[4];     // Position of the current cell (tau/t, x, y , eta/z=0)
0475     double LFSigma[4];      // LabFrame hypersurface (tau/t,x,y,eta/z=0), expect Sigma_mu
0476     double CMSigma[4];      // Center of mass hypersurface (tau/t,x,y,eta/z=0)
0477     double TRead;           // Temperature of cell
0478     double Vel[4];          // Gamma & 3-Velocity of cell (gamma, Vx, Vy, Vz) NOT FOUR VELOCITY
0479     double tau_pos;         // proper time from position of cell
0480     double eta_pos;         // eta from position of cell
0481     double tau_sur;         // proper time from normal vector of surface
0482     double eta_sur;         // eta from normal vector of surface
0483 
0484     // Calculated global quantities
0485     int PartCount;          // Total count of particles over ALL cells
0486     double NumLight;        // Number DENSITY of light quarks at set T
0487     double NumStrange;      // Number DENSITY of squarks at set T
0488     double UDDeg = 4.*6.;   // Degeneracy of UD quarks
0489     double OddDeg = 2.*6.;  // Degeneracy of S quarks
0490 
0491     // Calculated quantities in cells
0492     double LorBoost[4][4];  // Lorentz boost defined as used - Form is always Lambda_u^v
0493     int GeneratedParticles; // Number of particles to be generated this cell
0494 
0495     // Define parton densities
0496     double cut; // Each coordinate of P is integrated this far
0497     double E_light; // Store energy of light quarks
0498     double E_strange; // Store energy of strange quarks
0499     double GWeightProd; // Needed for integration
0500     double pSpatialdSigma; // Needed for integration
0501     PartCount = 0;
0502     NumLight = 0;       // Initialize density for light and strange quarks
0503     NumStrange = 0;
0504     GeneratedParticles = 0;
0505     double new_quark_energy; // store new quark energy for boost
0506 
0507     //counter for total number of light and strange quarks
0508     int nL_tot = 0;
0509     int nS_tot = 0;
0510 
0511     // define the accuracy range for temperature values in which the CDFGenerator is executed (for the case the value is not in the cache yet)
0512     const double accuracyRange = 0.003 / GEVFM; // for a usual hypersurface at const temperature there should be only one entry in the cache
0513     // precompute CDFs and store them in a cache
0514     std::unordered_map<double, std::vector<std::vector<double>>> cdfLightCache;
0515     std::unordered_map<double, std::vector<std::vector<double>>> cdfStrangeCache;
0516 
0517     // lambda function to check if a temperature value is within the accuracy range of a cached temperature
0518     auto withinAccuracyRange = [accuracyRange](double targetTemp, double cachedTemp) {
0519         return std::fabs(targetTemp - cachedTemp) <= accuracyRange;
0520     };
0521 
0522     // lambda function to retrieve the closest temperature from the cache within the accuracy range
0523     auto getClosestCachedTemp = [](const std::unordered_map<double, std::vector<std::vector<double>>>& cache, double targetTemp) {
0524         double closestTemp = -1.0;
0525         double minDistance = std::numeric_limits<double>::max();
0526         for (const auto& entry : cache) {
0527             double cachedTemp = entry.first;
0528             double distance = std::fabs(targetTemp - cachedTemp);
0529             if (distance < minDistance) {
0530                 minDistance = distance;
0531                 closestTemp = cachedTemp;
0532             }
0533         }
0534         return closestTemp;
0535     };
0536 
0537     // Define the lambda function to compute the energy
0538     auto computeEnergy = [](double mass, double cut, double GAbsL, double GAbsM, double GAbsK) {
0539         return sqrt(mass * mass + (cut * GAbsL) * (cut * GAbsL) + (cut * GAbsM) * (cut * GAbsM) + (cut * GAbsK) * (cut * GAbsK));
0540     };
0541 
0542     // Define the lambda function for the Fermi distribution
0543     auto FermiDistribution = [](double energy, double temperature) {
0544         return 1. / (exp(energy / temperature) + 1.);
0545     };
0546 
0547     for(int iS=0; iS<surface.size(); ++iS){
0548         tau_pos = surface[iS][0];
0549         CPos[1] = surface[iS][1];
0550         CPos[2] = surface[iS][2];
0551         eta_pos = surface[iS][3];  //we need t,x,y,z
0552         //this is also tau, x, y, eta
0553         tau_sur = surface[iS][4];
0554         LFSigma[1] = surface[iS][5];
0555         LFSigma[2] = surface[iS][6];
0556         eta_sur = surface[iS][7];
0557         TRead = surface[iS][8] / GEVFM;
0558         Vel[1] = surface[iS][9];
0559         Vel[2] = surface[iS][10];
0560         Vel[3] = surface[iS][11];
0561 
0562         cut = 10*TRead;
0563 
0564         // check if the CDFs for light quarks for this temperature are already in the cache and within the accuracy range
0565         auto cdfLightIter = cdfLightCache.find(TRead);
0566         if (cdfLightIter == cdfLightCache.end()) {
0567             // if not found, find the closest temperature in the cache within the accuracy range
0568             double closestTemp = getClosestCachedTemp(cdfLightCache,TRead);
0569 
0570             if (closestTemp >= 0. && withinAccuracyRange(TRead, closestTemp)) {
0571                 // use the CDFs from the closest temperature in the cache
0572                 CDFTabLight = cdfLightCache[closestTemp];
0573                 // set the current temperature value to the closest cached temperature
0574                 TRead = closestTemp;
0575             } else {
0576                 // if no suitable entry is found, compute the CDFs and store them in the cache
0577                 CDFGenerator(TRead, xmq, 1); // for light quarks
0578                 cdfLightCache[TRead] = CDFTabLight;
0579             }
0580         } else {
0581             // if found, use the precomputed CDFs from the cache
0582             CDFTabLight = cdfLightIter->second;
0583         }
0584 
0585         // check if the CDFs for strange quarks for this temperature are already in the cache and within the accuracy range
0586         auto cdfStrangeIter = cdfStrangeCache.find(TRead);
0587         if (cdfStrangeIter == cdfStrangeCache.end()) {
0588             // if not found, find the closest temperature in the cache within the accuracy range
0589             double closestTemp = getClosestCachedTemp(cdfStrangeCache,TRead);
0590 
0591             if (closestTemp >= 0 && withinAccuracyRange(TRead, closestTemp)) {
0592                 // use the CDFs from the closest temperature in the cache
0593                 CDFTabStrange = cdfStrangeCache[closestTemp];
0594                 // set the current temperature value to the closest cached temperature
0595                 TRead = closestTemp;
0596             } else {
0597                 // if no suitable entry is found, compute the CDFs and store them in the cache
0598                 CDFGenerator(TRead, xms, 2); // for strange quarks
0599                 cdfStrangeCache[TRead] = CDFTabStrange;
0600             }
0601         } else {
0602             // if found, use the precomputed CDFs from the cache
0603             CDFTabStrange = cdfStrangeIter->second;
0604         }
0605 
0606         if (Cartesian_hydro == false){
0607             //getting t,z from tau, eta
0608             CPos[0] = tau_pos*std::cosh(eta_pos);
0609             CPos[3] = tau_pos*std::sinh(eta_pos);
0610             //transform surface vector to Minkowski coordinates
0611             double cosh_eta_pos = std::cosh(eta_pos);
0612             double sinh_eta_pos = std::sinh(eta_pos);
0613             LFSigma[0] = cosh_eta_pos*tau_sur - (sinh_eta_pos / tau_pos) * eta_sur;
0614             LFSigma[3] = -sinh_eta_pos*tau_sur + (cosh_eta_pos / tau_pos) * eta_sur;
0615         }else{ // check later!
0616             CPos[0] = tau_pos;
0617             CPos[3] = eta_pos;
0618             LFSigma[0] = tau_sur;
0619             LFSigma[3] = eta_sur;
0620         }
0621 
0622         double vsquare = Vel[1]*Vel[1] + Vel[2]*Vel[2] + Vel[3]*Vel[3];
0623         if(vsquare < 10e-16){
0624             vsquare = 10e-16;
0625         }
0626 
0627         // Deduced info
0628         Vel[0] = 1. / sqrt(1-vsquare); // gamma - Vel is not four velocity
0629 
0630         // Lambda_u ^v
0631         LorBoost[0][0]= Vel[0];
0632         LorBoost[0][1]= Vel[0]*Vel[1];
0633         LorBoost[0][2]= Vel[0]*Vel[2];
0634         LorBoost[0][3]= Vel[0]*Vel[3];
0635         LorBoost[1][0]= Vel[0]*Vel[1];
0636         LorBoost[1][1]= (Vel[0] - 1.)*Vel[1]*Vel[1]/vsquare + 1.;
0637         LorBoost[1][2]= (Vel[0] - 1.)*Vel[1]*Vel[2]/vsquare;
0638         LorBoost[1][3]= (Vel[0] - 1.)*Vel[1]*Vel[3]/vsquare;
0639         LorBoost[2][0]= Vel[0]*Vel[2];
0640         LorBoost[2][1]= (Vel[0] - 1.)*Vel[1]*Vel[2]/vsquare;
0641         LorBoost[2][2]= (Vel[0] - 1.)*Vel[2]*Vel[2]/vsquare + 1.;
0642         LorBoost[2][3]= (Vel[0] - 1.)*Vel[2]*Vel[3]/vsquare;
0643         LorBoost[3][0]= Vel[0]*Vel[3];
0644         LorBoost[3][1]= (Vel[0] - 1.)*Vel[1]*Vel[3]/vsquare;
0645         LorBoost[3][2]= (Vel[0] - 1.)*Vel[2]*Vel[3]/vsquare;
0646         LorBoost[3][3]= (Vel[0] - 1.)*Vel[3]*Vel[3]/vsquare + 1.;
0647 
0648         if(vsquare == 0){
0649             LorBoost[1][1]= 1.;
0650             LorBoost[1][2]= 0;
0651             LorBoost[1][3]= 0;
0652             LorBoost[2][1]= 0;
0653             LorBoost[2][2]= 1.;
0654             LorBoost[2][3]= 0;
0655             LorBoost[3][1]= 0;
0656             LorBoost[3][2]= 0;
0657             LorBoost[3][3]= 1.;
0658         }
0659         // Lambda_u^v Sigma_v = CMSigma_u
0660         CMSigma[0] = (LorBoost[0][0]*LFSigma[0] + LorBoost[0][1]*LFSigma[1] + LorBoost[0][2]*LFSigma[2] + LorBoost[0][3]*LFSigma[3]);
0661         CMSigma[1] = (LorBoost[1][0]*LFSigma[0] + LorBoost[1][1]*LFSigma[1] + LorBoost[1][2]*LFSigma[2] + LorBoost[1][3]*LFSigma[3]);
0662         CMSigma[2] = (LorBoost[2][0]*LFSigma[0] + LorBoost[2][1]*LFSigma[1] + LorBoost[2][2]*LFSigma[2] + LorBoost[2][3]*LFSigma[3]);
0663         CMSigma[3] = (LorBoost[3][0]*LFSigma[0] + LorBoost[3][1]*LFSigma[1] + LorBoost[3][2]*LFSigma[2] + LorBoost[3][3]*LFSigma[3]);
0664 
0665         // GAUSSIAN INTEGRALS <n> = int f(p)d3p
0666         NumLight = 0.;
0667         NumStrange = 0.;
0668 
0669         for (int l = 0; l < GPoints; l++) {
0670             for (int m = 0; m < GPoints; m++) {
0671                 for (int k = 0; k < GPoints; k++) {
0672                     GWeightProd = GWeight[l] * GWeight[m] * GWeight[k];
0673                     pSpatialdSigma = (cut * GAbs[l]) * CMSigma[1] + (cut * GAbs[m]) * CMSigma[2] + (cut * GAbs[k]) * CMSigma[3];
0674 
0675                     // For UD quarks
0676                     E_light = computeEnergy(xmq, cut, GAbs[l], GAbs[m], GAbs[k]);
0677                     NumLight += GWeightProd * FermiDistribution(E_light, TRead) *
0678                                (E_light * CMSigma[0] + pSpatialdSigma) / E_light;
0679 
0680                     // For S quarks
0681                     E_strange = computeEnergy(xms, cut, GAbs[l], GAbs[m], GAbs[k]);
0682                     NumStrange += GWeightProd * FermiDistribution(E_strange, TRead) *
0683                                  (E_strange * CMSigma[0] + pSpatialdSigma) / E_strange;
0684                 }
0685             }
0686         }
0687 
0688         // U, D, UBAR, DBAR QUARKS
0689         // <N> = V <n>
0690         double NumHere = NumLight*UDDeg*cut*cut*cut/(8.*PI*PI*PI);
0691         std::poisson_distribution<int> poisson_ud(NumHere);
0692 
0693         // Generating light quarks
0694         GeneratedParticles = poisson_ud(getRandomGenerator()); // Initialize particles created in this cell
0695 
0696         // List of particles ( pos(x,y,z,t), mom(px,py,pz,E), species)
0697         for(int partic = 0; partic < GeneratedParticles; partic++){
0698             // adding space to PList for output quarks
0699             std::vector<double> temp = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
0700             Plist.push_back(temp);
0701 
0702             // Species - U,D,UBar,Dbar are equally likely
0703             double SpecRoll = ran(); //Probability of species die roll
0704             if (SpecRoll <= 0.25) { // UBar
0705                 Plist[PartCount][1] = -2;
0706             } else if (SpecRoll <= 0.50) { // DBar
0707                 Plist[PartCount][1] = -1;
0708             } else if (SpecRoll <= 0.75) { // D
0709                 Plist[PartCount][1] = 1;
0710             } else { // U
0711                 Plist[PartCount][1] = 2;
0712             }
0713 
0714             // Position
0715             // Located at x,y pos of area element
0716             Plist[PartCount][10] = CPos[0] + (ran() - 0.5)*CellDT; // Tau
0717             Plist[PartCount][7] = CPos[1] + (ran() - 0.5)*CellDX;
0718             Plist[PartCount][8] = CPos[2] + (ran() - 0.5)*CellDY;
0719             Plist[PartCount][9] = CPos[3] + (ran() - 0.5)*CellDZ;
0720 
0721             // Momentum
0722             // Sample rest frame momentum given T and mass of light quark
0723             MCSampler(TRead, 1); // NewP=P, NewX=Px, ...
0724 
0725             // USE THE SAME BOOST AS BEFORE
0726             // PLab^u = g^u^t Lambda_t ^v pres^w g_w _v
0727             // Returns P in GeV
0728             new_quark_energy = sqrt(xmq*xmq + NewP*NewP);
0729             Plist[PartCount][6] = (LorBoost[0][0]*new_quark_energy + LorBoost[0][1]*NewX + LorBoost[0][2]*NewY + LorBoost[0][3]*NewZ)*GEVFM;
0730             Plist[PartCount][3] = (LorBoost[1][0]*new_quark_energy + LorBoost[1][1]*NewX + LorBoost[1][2]*NewY + LorBoost[1][3]*NewZ)*GEVFM;
0731             Plist[PartCount][4] = (LorBoost[2][0]*new_quark_energy + LorBoost[2][1]*NewX + LorBoost[2][2]*NewY + LorBoost[2][3]*NewZ)*GEVFM;
0732             Plist[PartCount][5] = (LorBoost[3][0]*new_quark_energy + LorBoost[3][1]*NewX + LorBoost[3][2]*NewY + LorBoost[3][3]*NewZ)*GEVFM;
0733 
0734             // Additional information
0735             Plist[PartCount][0]  = 1; // Event ID, to match jet formatting
0736             Plist[PartCount][2]  = 0; // Origin, to match jet formatting
0737             Plist[PartCount][11] = 0; // Status - identifies as thermal quark
0738             PartCount++;
0739         }
0740 
0741         // S, SBAR QUARKS
0742         // <N> = V <n>
0743         double NumOddHere = NumStrange*OddDeg*cut*cut*cut/(8.*PI*PI*PI);
0744         std::poisson_distribution<int> poisson_s(NumOddHere);
0745 
0746         // Generate s quarks
0747         int nL = GeneratedParticles;
0748         nL_tot += nL;
0749 
0750         GeneratedParticles = poisson_s(getRandomGenerator()); //Initialize particles created in this cell
0751 
0752         int nS = GeneratedParticles;
0753         nS_tot += nS;
0754         //List of particles ( pos(x,y,z,t), mom(px,py,pz,E), species)
0755         for(int partic = 0; partic < GeneratedParticles; partic++){
0756             // adding space to PList for output quarks
0757             std::vector<double> temp = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
0758             Plist.push_back(temp);
0759 
0760             // Species - S,Sbar are equally likely
0761             double SpecRoll = ran();
0762             if(SpecRoll <= 0.5){ // SBar
0763                 Plist[PartCount][1] = -3;
0764             } else { //S
0765                 Plist[PartCount][1] = 3;
0766             }
0767 
0768             // Position
0769             // Located at x,y pos of area element
0770             Plist[PartCount][10] = CPos[0] + (ran() - 0.5)*CellDT; // Tau
0771             Plist[PartCount][7] = CPos[1] + (ran() - 0.5)*CellDX;
0772             Plist[PartCount][8] = CPos[2] + (ran() - 0.5)*CellDY;
0773             Plist[PartCount][9] = CPos[3] + (ran() - 0.5)*CellDZ;
0774 
0775             // Momentum
0776             // Sample rest frame momentum given T and mass of s quark
0777             MCSampler(TRead, 2); // NewP=P, NewX=Px, ...
0778 
0779             // USE THE SAME BOOST AS BEFORE
0780             // PLab^u = g^u^t Lambda_t ^v pres^w g_w _v
0781             // Returns P in GeV
0782             new_quark_energy = sqrt(xms*xms + NewP*NewP);
0783             Plist[PartCount][6] = (LorBoost[0][0]*new_quark_energy + LorBoost[0][1]*NewX + LorBoost[0][2]*NewY + LorBoost[0][3]*NewZ)*GEVFM;
0784             Plist[PartCount][3] = (LorBoost[1][0]*new_quark_energy + LorBoost[1][1]*NewX + LorBoost[1][2]*NewY + LorBoost[1][3]*NewZ)*GEVFM;
0785             Plist[PartCount][4] = (LorBoost[2][0]*new_quark_energy + LorBoost[2][1]*NewX + LorBoost[2][2]*NewY + LorBoost[2][3]*NewZ)*GEVFM;
0786             Plist[PartCount][5] = (LorBoost[3][0]*new_quark_energy + LorBoost[3][1]*NewX + LorBoost[3][2]*NewY + LorBoost[3][3]*NewZ)*GEVFM;
0787 
0788             // Additional information
0789             Plist[PartCount][0]  = 1; // Event ID, to match jet formatting
0790             Plist[PartCount][2]  = 0; // Origin, to match jet formatting
0791             Plist[PartCount][11] = 0; // Status - identifies as thermal quark
0792             PartCount++;
0793         }
0794     }
0795 
0796     JSDEBUG << "Light particles: " << nL_tot;
0797     JSDEBUG << "Strange particles: " << nS_tot;
0798 
0799     num_ud = nL_tot;
0800     num_s = nS_tot;
0801 
0802     //Shuffling PList
0803     if(ShuffleList){
0804         // Shuffle the Plist using the random engine
0805         std::shuffle(&Plist[0], &Plist[PartCount], getRandomGenerator());
0806     }
0807 
0808     //print Plist for testing
0809     /*std::cout << std::setprecision(5);
0810     std::ofstream thermalP;
0811     thermalP.open("thermal_partons.dat", std::ios::app);
0812     for(int p=0; p < Plist.size(); p++){
0813         thermalP << Plist[p][0] << " " << Plist[p][1] << " " << Plist[p][2]
0814         << " " << Plist[p][3] << " " << Plist[p][4] << " " << Plist[p][5]
0815         << " " << Plist[p][6] << " " << Plist[p][7] << " " << Plist[p][8]
0816         << " " << Plist[p][9] << " " << Plist[p][10] << " " << Plist[p][11]
0817         << "\n";
0818     }
0819     thermalP.close();*/
0820 }
0821 
0822 void ThermalPartonSampler::sample_2p1d(double eta_max){
0823 
0824     // Input read from cells
0825     double CPos[4];     // Position of the current cell (tau/t, x, y , eta/z=0)
0826     double LFSigma[4];      // LabFrame hypersurface (tau/t,x,y,eta/z=0)
0827     double CMSigma[4];      // Center of mass hypersurface (tau/t,x,y,eta/z=0)
0828     double TRead;           // Temperature of cell
0829     double Vel[4];          // Gamma & 3-Velocity of cell (gamma, Vx, Vy, Vz) NOT FOUR VELOCITY
0830     double tau_pos;         // proper time from position of cell
0831     double eta_pos;         // eta from position of cell
0832     double tau_sur;         // proper time from normal vector of surface
0833     double eta_sur;         // eta from normal vector of surface
0834 
0835     // Calculated global quantities
0836     int PartCount;          // Total count of particles over ALL cells
0837     double NumLight;        // Number DENSITY of light quarks at set T
0838     double NumStrange;      // Number DENSITY of squarks at set T
0839     double UDDeg = 4.*6.;   // Degeneracy of UD quarks
0840     double OddDeg = 2.*6.;  // Degeneracy of S quarks
0841 
0842     // Calculated quantities in cells
0843     double LorBoost[4][4];  // Lorentz boost defined as used - Form is always Lambda_u^v
0844     int GeneratedParticles; // Number of particles to be generated this cell
0845 
0846     // Define parton densities
0847     double cut; // Each coordinate of P is integrated this far
0848     double E_light; // Store energy of light quarks
0849     double E_strange; // Store energy of strange quarks
0850     double GWeightProd; // Needed for integration
0851     double pSpatialdSigma; // Needed for integration
0852     PartCount = 0;
0853     NumLight = 0;       // Initialize density for light and strange quarks
0854     NumStrange = 0;
0855     GeneratedParticles = 0;
0856     double new_quark_energy; // store new quark energy for boost
0857 
0858     //counter for total number of light and strange quarks
0859     int nL_tot = 0;
0860     int nS_tot = 0;
0861 
0862     // define the accuracy range for temperature values in which the CDFGenerator is executed (for the case the value is not in the cache yet)
0863     const double accuracyRange = 0.003 / GEVFM; // for a usual hypersurface at const temperature there should be only one entry in the cache
0864     // precompute CDFs and store them in a cache
0865     std::unordered_map<double, std::vector<std::vector<double>>> cdfLightCache;
0866     std::unordered_map<double, std::vector<std::vector<double>>> cdfStrangeCache;
0867 
0868     // lambda function to check if a temperature value is within the accuracy range of a cached temperature
0869     auto withinAccuracyRange = [accuracyRange](double targetTemp, double cachedTemp) {
0870         return std::fabs(targetTemp - cachedTemp) <= accuracyRange;
0871     };
0872 
0873     // lambda function to retrieve the closest temperature from the cache within the accuracy range
0874     auto getClosestCachedTemp = [](const std::unordered_map<double, std::vector<std::vector<double>>>& cache, double targetTemp) {
0875         double closestTemp = -1.0;
0876         double minDistance = std::numeric_limits<double>::max();
0877         for (const auto& entry : cache) {
0878             double cachedTemp = entry.first;
0879             double distance = std::fabs(targetTemp - cachedTemp);
0880             if (distance < minDistance) {
0881                 minDistance = distance;
0882                 closestTemp = cachedTemp;
0883             }
0884         }
0885         return closestTemp;
0886     };
0887 
0888     // Define the lambda function to compute the energy
0889     auto computeEnergy = [](double mass, double cut, double GAbsL, double GAbsM, double GAbsK) {
0890         return sqrt(mass * mass + (cut * GAbsL) * (cut * GAbsL) + (cut * GAbsM) * (cut * GAbsM) + (cut * GAbsK) * (cut * GAbsK));
0891     };
0892 
0893     // Define the lambda function for the Fermi distribution
0894     auto FermiDistribution = [](double energy, double temperature) {
0895         return 1. / (exp(energy / temperature) + 1.);
0896     };
0897 
0898     std::vector<double> NumLightList;
0899     std::vector<double> NumStrangeList;
0900 
0901     double d_eta = CellDZ;
0902 
0903     int N_slices = std::floor((eta_max / d_eta)-0.5)+1;
0904     for(int slice=1; slice <= (2*N_slices+1); slice++){
0905         double eta_slice = (slice-N_slices-1)*d_eta;
0906 
0907         JSINFO << "Sampling thermal partons for pseudorapidity slice " << slice
0908         << " (eta_min = " << eta_slice-(d_eta/2.) << ", eta_max = "
0909         << eta_slice+(d_eta/2.) << ")";
0910 
0911         for(int iS=0; iS<surface.size(); ++iS){
0912             tau_pos = surface[iS][0];
0913             CPos[1] = surface[iS][1];
0914             CPos[2] = surface[iS][2];
0915             eta_pos = surface[iS][3];  //we need t,x,y,z
0916             //this is also tau, x, y, eta
0917             tau_sur = surface[iS][4];
0918             LFSigma[1] = surface[iS][5];
0919             LFSigma[2] = surface[iS][6];
0920             eta_sur = surface[iS][7];
0921             TRead = surface[iS][8] / GEVFM;
0922             Vel[1] = surface[iS][9];
0923             Vel[2] = surface[iS][10];
0924             Vel[3] = surface[iS][11];
0925 
0926             cut = 10*TRead;
0927 
0928             // check if the CDFs for light quarks for this temperature are already in the cache and within the accuracy range
0929             auto cdfLightIter = cdfLightCache.find(TRead);
0930             if (cdfLightIter == cdfLightCache.end()) {
0931                 // if not found, find the closest temperature in the cache within the accuracy range
0932                 double closestTemp = getClosestCachedTemp(cdfLightCache,TRead);
0933 
0934                 if (closestTemp >= 0. && withinAccuracyRange(TRead, closestTemp)) {
0935                     // use the CDFs from the closest temperature in the cache
0936                     CDFTabLight = cdfLightCache[closestTemp];
0937                     // set the current temperature value to the closest cached temperature
0938                     TRead = closestTemp;
0939                 } else {
0940                     // if no suitable entry is found, compute the CDFs and store them in the cache
0941                     CDFGenerator(TRead, xmq, 1); // for light quarks
0942                     cdfLightCache[TRead] = CDFTabLight;
0943                 }
0944             } else {
0945                 // if found, use the precomputed CDFs from the cache
0946                 CDFTabLight = cdfLightIter->second;
0947             }
0948 
0949             // check if the CDFs for strange quarks for this temperature are already in the cache and within the accuracy range
0950             auto cdfStrangeIter = cdfStrangeCache.find(TRead);
0951             if (cdfStrangeIter == cdfStrangeCache.end()) {
0952                 // if not found, find the closest temperature in the cache within the accuracy range
0953                 double closestTemp = getClosestCachedTemp(cdfStrangeCache,TRead);
0954 
0955                 if (closestTemp >= 0 && withinAccuracyRange(TRead, closestTemp)) {
0956                     // use the CDFs from the closest temperature in the cache
0957                     CDFTabStrange = cdfStrangeCache[closestTemp];
0958                     // set the current temperature value to the closest cached temperature
0959                     TRead = closestTemp;
0960                 } else {
0961                     // if no suitable entry is found, compute the CDFs and store them in the cache
0962                     CDFGenerator(TRead, xms, 2); // for strange quarks
0963                     cdfStrangeCache[TRead] = CDFTabStrange;
0964                 }
0965             } else {
0966                 // if found, use the precomputed CDFs from the cache
0967                 CDFTabStrange = cdfStrangeIter->second;
0968             }
0969 
0970             //getting t,z from tau, eta
0971             CPos[0] = tau_pos*std::cosh(eta_pos);
0972             CPos[3] = tau_pos*std::sinh(eta_pos);
0973             //transform surface vector to Minkowski coordinates
0974             double cosh_eta_pos = std::cosh(eta_pos);
0975             double sinh_eta_pos = std::sinh(eta_pos);
0976             LFSigma[0] = cosh_eta_pos*tau_sur - (sinh_eta_pos / tau_pos) * eta_sur;
0977             LFSigma[3] = -sinh_eta_pos*tau_sur + (cosh_eta_pos / tau_pos) * eta_sur;
0978 
0979             CellDZ = CPos[0] * 2. * std::sinh(d_eta/2.);
0980 
0981             double vsquare = Vel[1]*Vel[1] + Vel[2]*Vel[2] + Vel[3]*Vel[3];
0982             if(vsquare < 10e-16){
0983                 vsquare = 10e-16;
0984             }
0985 
0986             // Deduced info
0987             Vel[0] = 1. / sqrt(1-vsquare); // gamma - Vel is not four velocity
0988 
0989             // Lambda_u ^v
0990             LorBoost[0][0]= Vel[0];
0991             LorBoost[0][1]= Vel[0]*Vel[1];
0992             LorBoost[0][2]= Vel[0]*Vel[2];
0993             LorBoost[0][3]= Vel[0]*Vel[3];
0994             LorBoost[1][0]= Vel[0]*Vel[1];
0995             LorBoost[1][1]= (Vel[0] - 1.)*Vel[1]*Vel[1]/vsquare + 1.;
0996             LorBoost[1][2]= (Vel[0] - 1.)*Vel[1]*Vel[2]/vsquare;
0997             LorBoost[1][3]= (Vel[0] - 1.)*Vel[1]*Vel[3]/vsquare;
0998             LorBoost[2][0]= Vel[0]*Vel[2];
0999             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 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!