DEMPgen/src/MatterEffects.cxx

0001 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0002   Rory Evans, Feb. 2018.
0003   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
0004 #include <stdio.h>
0005 #include <iostream>
0006 #include <string.h>
0007
0008 #include "TMath.h"
0009 #include "TRandom3.h"
0010 #include "Particle.hxx"
0011 #include "MatterEffects.hxx"
0012 #include "Constants.hxx"
0013
0014 using namespace TMath;
0015
0016 void MatterEffects::MultiScatter(Particle* P, double radlen)
0017 {
0018
0019   Particle * Pms = new Particle(P->M(), 0, 0, P->P());
0020   //std::cout << P->GetPid()<<"\t"<<P->Theta() << std::endl;
0021   beta = P->P()/Sqrt(Power(P->M(),2)+Power(P->P(),2));
0022   //std::cout << beta << std::endl;
0023   mstheta0 = 13.6/(beta*(P->P()))*Abs(P->GetCharge())*Sqrt(radlen)
0024     * (1+0.038*Log(radlen));
0025
0026   //std::cout << mstheta0 << std::endl;
0027   mstheta = gRandom->Gaus(0, mstheta0);
0028   //std::cout << mstheta << std::endl;
0029
0030
0031   msphi = Pi()*gRandom->Rndm() - Pi()/2;
0032
0033   Pms->RotateY(mstheta);
0034   Pms->RotateZ(msphi);
0035
0036   Pms->RotateY(P->Theta());
0037   Pms->RotateZ(P->Phi());
0038
0039   *P = *Pms;
0040
0041   delete Pms;
0042
0043 }
0044
0045 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0046   Following functions adapted from eloss.h from
0047   Zafar Ahmed's Generator (in turn adapted from SIMC)
0048   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
0049
0050 double MatterEffects::eta(double aZ)
0051 {
0052   //Phys.Rev.D 12,1884 A46
0053   return log( 1440 * pow( aZ , -2. / 3. ) ) / log( 183 * pow( aZ , -1. / 3. ) );
0054 }
0055
0056 double MatterEffects::b(double aZ)
0057 {
0058   //Phys.Rev.D 12,1884 A45
0059   if ( aZ!=0 )
0060     return 4. / 3. * ( 1 + 1. / 9. * ( aZ + 1 ) / ( aZ + eta( aZ ) ) / log( 183 * pow( aZ , -1/3. ) ) );
0061   else
0062     return 0;
0063 }
0064
0065 double MatterEffects::X0(double aZ, double aA)
0066 {
0067   if ( aZ==1 )
0068     {
0069       Lrad=5.31;
0070       Lrad_prime=6.144;
0071     }
0072   else if ( aZ==2 )
0073     {
0074       Lrad=4.79;
0075       Lrad_prime=5.621;
0076     }
0077   else if ( aZ==3 )
0078     {
0079       Lrad=4.74;
0080       Lrad_prime=5.805;
0081     }
0082   else if ( aZ==4 )
0083     {
0084       Lrad=4.71;
0085       Lrad_prime=5.924;
0086     }
0087   else
0088     {
0089       Lrad       = log( 184.15 * pow( aZ , -1./3.));
0090       Lrad_prime = log( 1194.  * pow( aZ , -2./3.));
0091     }
0092
0093   a = constants::alpha * aZ;
0094
0095   a = a * a;
0096
0097   f_Z = a * ( 1./(1+a) + 0.20206 - 0.0369 * a + 0.0083 * a * a - 0.002 * a * a * a );
0098
0099   if ( aZ!=0 && aA!=0 ) {
0100     return 716.408 * aA / ( aZ * aZ * ( Lrad - f_Z ) + aZ * Lrad_prime );
0101   }
0102   else
0103     return 0;
0104
0105 }
0106
0107 void MatterEffects::BremLoss( Particle* P, Double_t abt)
0108 {
0109   aE0 = P->E();
0110
0111   //std::cout << aE0 << "\t";
0112
0113   result = 0;
0114   if ( abt!=0 )
0115     result = aE0 * pow( gRandom->Rndm() * 0.999 , 1. / abt );
0116
0117   if ( result > ( aE0 - constants::electron_mass_mev) )
0118     result = aE0 - constants::electron_mass_mev;
0119
0120   if ( result < 0 )
0121     result = 0;
0122
0123   //std::cout << result << std::endl;
0124
0125   P->SetThetaPhiE(P->Theta(),P->Phi(),P->E()-result);
0126
0127 }
0128
0129 void MatterEffects::IonLoss(Particle* P, double a, double z, double rho, double t )
0130 {
0131
0132   mass = P->GetMass();
0133   aE0 = P->E();
0134
0135   //std::cout << aE0 << "\t";
0136
0137   lK = 0.307075;  // Page 398. cm^2/g for A=1 g/mol
0138
0139   lbetasq = 1 - mass * mass / ( aE0 * aE0 );
0140
0141   lxi = lK / 2 * z / a * t / lbetasq;  //aT: g/cm^2
0142
0143   lhbarwsq = 28.816 * 28.816 * rho * z / a * 1e-12;  //Page 398 MeV arho is density of absorber
0144
0145   j = 0.200;
0146
0147   Delta_p = lxi * ( log( 2 * mass * lxi / lhbarwsq ) + j ); // Equation 32.12
0148
0149   lw = 4 * lxi;
0150
0151   result = 0;
0152
0153   if ( z != 0 && a != 0 && t != 0 && rho != 0 )
0154
0155     result = gRandom->Landau( Delta_p , lw );
0156
0157   if ( result > ( aE0 - mass ) )
0158     result = aE0 - mass;
0159
0160   if ( result < 0 )
0161     result = 0;
0162
0163   //td::cout << result << std::endl;
0164
0165   P->SetThetaPhiE(P->Theta(),P->Phi(),P->E()-result);
0166
0167   /*
0168     char *p_name = P->GetName();
0169     char ms_name[100];
0170
0171     strcpy(ms_name, p_name);
0172     strcat(ms_name, "_Ion");
0173
0174     //Pms->SetName(ms_name);
0175     */
0176
0177 }