estarlight/src/beambeamsystem.cpp

0001 ///////////////////////////////////////////////////////////////////////////
0002 //
0003 //    Copyright 2010
0004 //
0005 //    This file is part of starlight.
0006 //
0007 //    starlight is free software: you can redistribute it and/or modify
0008 //    it under the terms of the GNU General Public License as published by
0009 //    the Free Software Foundation, either version 3 of the License, or
0010 //    (at your option) any later version.
0011 //
0012 //    starlight is distributed in the hope that it will be useful,
0013 //    but WITHOUT ANY WARRANTY; without even the implied warranty of
0014 //    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0015 //    GNU General Public License for more details.
0016 //
0017 //    You should have received a copy of the GNU General Public License
0018 //    along with starlight. If not, see <http://www.gnu.org/licenses/>.
0019 //
0020 ///////////////////////////////////////////////////////////////////////////
0021 //
0022 // File and Version Information:
0023 // $Rev:: 260                         $: revision of last commit
0024 // $Author:: butter                   $: author of last commit
0025 // $Date:: 2016-05-03 04:07:34 +0100 #$: date of last commit
0026 //
0027 // Description:
0028 //
0029 //
0030 //
0031 ///////////////////////////////////////////////////////////////////////////
0032
0033
0034 #include <iostream>
0035 #include <fstream>
0036 #include <cmath>
0037
0038 #include "inputParameters.h"
0039 #include "reportingUtils.h"
0040 #include "starlightconstants.h"
0041 #include "bessel.h"
0042 #include "beambeamsystem.h"
0043
0044
0045 using namespace std;
0046 using namespace starlightConstants;
0047
0048
0049 //______________________________________________________________________________
0050 beamBeamSystem::beamBeamSystem(const inputParameters& inputParametersInstance,
0051                    const beam&            electronBeam,
0052                                const beam&            targetBeam)
0053   : _beamLorentzGamma(inputParametersInstance.beamLorentzGamma()),
0054     _electronBeamLorentzGamma(inputParametersInstance.electronBeamLorentzGamma()),
0055     _targetBeamLorentzGamma(inputParametersInstance.targetBeamLorentzGamma()),
0056     _beamBreakupMode (inputParametersInstance.beamBreakupMode()),
0057     _electronBeam           (electronBeam),
0058     _targetBeam           (targetBeam),
0059     _breakupProbabilities(0)
0060 //    _breakupImpactParameterStep(1.007),
0061 //    _breakupCutOff(10e-6)
0062 {
0063   init();
0064 }
0065
0066
0067
0068
0069 //______________________________________________________________________________
0070 beamBeamSystem::beamBeamSystem(const inputParameters& inputParametersInstance)
0071     : _beamLorentzGamma(inputParametersInstance.beamLorentzGamma()),
0072           _electronBeamLorentzGamma(inputParametersInstance.electronBeamLorentzGamma()),
0073           _targetBeamLorentzGamma(inputParametersInstance.targetBeamLorentzGamma()),
0074       _beamBreakupMode (inputParametersInstance.beamBreakupMode()),
0075       _electronBeam           (1,
0076                    0,
0077                 inputParametersInstance.productionMode(),
0078                             inputParametersInstance.electronBeamLorentzGamma()),
0079       _targetBeam           (inputParametersInstance.targetBeamZ(),
0080                         inputParametersInstance.targetBeamA(),
0081                 inputParametersInstance.productionMode(),
0082                             inputParametersInstance.targetBeamLorentzGamma()),
0083       _breakupProbabilities(0)
0084 //    _breakupImpactParameterStep(1.007),
0085 //    _breakupCutOff(10e-10)
0086 {
0087   init();
0088 }
0089
0090
0091
0092 //______________________________________________________________________________
0093 beamBeamSystem::~beamBeamSystem()
0094 { }
0095
0096 void beamBeamSystem::init()
0097 {
0098    // Calculate beam gamma in CMS frame
0099    double rap1 = acosh(_electronBeamLorentzGamma);
0100    double rap2 = -acosh(_targetBeamLorentzGamma);
0101
0102    _cmsBoost = (rap1+rap2)/2.;
0103
0104    _beamLorentzGamma = cosh((rap1-rap2)/2);
0105    _electronBeam.setBeamLorentzGamma(_beamLorentzGamma);
0106    _targetBeam.setBeamLorentzGamma(_beamLorentzGamma);
0107
0108    generateBreakupProbabilities();
0109 }
0110 //______________________________________________________________________________
0111 double
0112 beamBeamSystem::probabilityOfBreakup(const double D) const
0113 {
0114   cout << "D=" << D << " NOT SUPPORTED IN eX COLLISIONS" << endl;
0115   return 1;
0116 }
0117
0118
0119 void
0120 beamBeamSystem::generateBreakupProbabilities()
0121 {
0122
0123   /*
0124     double bMin = (_beam1.nuclearRadius()+_beam2.nuclearRadius())/2.;
0125
0126     // Do this only for nucleus-nucleus collisions.
0127     // pp and pA are handled directly in probabilityOfBreakup
0128 //    if ((_beam1.Z() != 1) && (_beam1.A() != 1) && (_beam2.Z() != 1) && _beam2.A() != 1) {
0129 //  this change allows deuterium and tritium to be handled here
0130     if ((_beam1.Z() != 1) && (_beam1.A() != 1) && _beam2.A() != 1) {
0131
0132         if (_beamBreakupMode == 1)
0133             printInfo << "Hard Sphere Break criteria. b > " << 2. * _beam1.nuclearRadius() << endl;
0134         if (_beamBreakupMode == 2)
0135             printInfo << "Requiring XnXn [Coulomb] breakup. " << endl;
0136         if (_beamBreakupMode == 3)
0137             printInfo << "Requiring 1n1n [Coulomb only] breakup. " << endl;
0138         if (_beamBreakupMode == 4)
0139             printInfo << "Requiring both nuclei to remain intact. " << endl;
0140         if (_beamBreakupMode == 5)
0141             printInfo << "Requiring no hadronic interactions. " << endl;
0142         if (_beamBreakupMode == 6)
0143             printInfo << "Requiring breakup of one or both nuclei. " << endl;
0144         if (_beamBreakupMode == 7)
0145             printInfo << "Requiring breakup of one nucleus (Xn,0n). " << endl;
0146
0147     double pOfB = 0;
0148     double b = bMin;
0149         double totRad = _beam1.nuclearRadius()+_beam2.nuclearRadius();
0150
0151     while(1)
0152     {
0153
0154             if(_beamBreakupMode != 5)
0155             {
0156                 if(b > (totRad*1.5))
0157                 {
0158                     if(pOfB<_breakupCutOff)
0159                     {
0160 //                         std::cout << "Break off b: " << b << std::endl;
0161 //                         std::cout << "Number of PofB bins: " << _breakupProbabilities.size() << std::endl;
0162                         break;
0163                     }
0164                 }
0165             }
0166             else
0167             {
0168                 if((1-pOfB)<_breakupCutOff)
0169                 {
0170 //                         std::cout << "Break off b: " << b << std::endl;
0171 //                         std::cout << "Number of PofB bins: " << _breakupProbabilities.size() << std::endl;
0172                         break;
0173                 }
0174             }
0175 //             std::cout << 1-pOfBreakup << std::endl;
0176
0177             probabilityOfHadronBreakup(b);
0178             probabilityOfPhotonBreakup(b, _beamBreakupMode);
0179
0180             //What was probability of photonbreakup depending upon mode selection,
0181             // is now done in the photonbreakupfunction
0182             if (_beamBreakupMode == 1) {
0183                 if (b >_beam1.nuclearRadius()+_beam2.nuclearRadius())  // symmetry
0184                     _pHadronBreakup = 0;
0185                 else
0186                     _pHadronBreakup = 999.;
0187             }
0188
0189             b *= _breakupImpactParameterStep;
0190         pOfB = exp(-1 * _pHadronBreakup) * _pPhotonBreakup;
0191             _breakupProbabilities.push_back(pOfB);
0192         } // End while(1)
0193     }
0194   */
0195 }
0196
0197 //______________________________________________________________________________
0198 double
0199 beamBeamSystem::probabilityOfHadronBreakup(const double impactparameter)
0200 {
0201     //probability of hadron breakup,
0202     //this is what is returned when the function is called
0203     double gamma = _beamLorentzGamma;
0204     //input for gamma_em
0205     double b = impactparameter;
0206     int a1 = 0;
0207     int a2 = _targetBeam.A();
0208
0209     static int IFIRSTH = 0;
0210     static double DELL=0., DELR=0., SIGNN=0., R1=0., A1=0., A2=0., R2=0., RHO1=0.;
0211     static double RHO2=0., NZ1=0., NZ2=0., NR1=0., NR2=0.,RR1=0., RR2=0., NY=0., NX=0.;
0212     static double AN1=0., AN2=0.;
0213     double RSQ=0.,Z1=0.,Z2=0.,Y=0.,X=0.,XB=0.,RPU=0.,IRUP=0.,RTU=0.;
0214     double IRUT=0.,T1=0.,T2=0.;
0215     static double DEN1[20002], DEN2[20002];
0216     double energy,sigmainmb;
0217     if (IFIRSTH != 0) goto L100;
0218     //Initialize
0219     //Integration delta x, delta z
0220     IFIRSTH = 1;
0221     DELL   = .05;
0222     DELR   = .01;
0223
0224     // replace this with a parameterization from the particle data book  SRK 4/1025
0225     //use two sigma_NN's. 52mb at rhic 100gev/beam, 88mb at LHC 2.9tev/beam, gamma is in cm system
0226     //SIGNN = 5.2;
0227     //if ( gamma > 500. ) SIGNN = 8.8;
0228     energy=2*gamma*0.938;   // center of mass energy, in GeV
0229       // This equation is from section 50 of the particle data book, the subsection on "Total Hadronic Cross-Sections, using the parameterization for sqrt{s} > 7 GeV.
0230       // only the first and second terms contribute significantly, but leave them all here for good measure
0231       sigmainmb = 0.2838*pow(log(energy),2)+33.73+13.67*pow(energy,-0.412)-7.77*pow(energy,-0.5626);
0232       SIGNN=sigmainmb/10.;
0233
0234     //use parameter from Constants
0235     R1 = ( 0 );
0236         R2 = ( _targetBeam.nuclearRadius());
0237     A1 = 0.535; //This is woodsaxonskindepth
0238         A2 = 0.535;
0239     //write(6,12)r1,a1,signn  Here is where we could probably set this up asymmetrically R2=_beam2.nuclearRadius() and RHO2=ap2=_beam2.A()
0240     // R2 = R1;
0241     RHO1 = a1;
0242     RHO2 = a2;
0243     NZ1  = ((R1+5.)/DELR);
0244     NR1  = NZ1;
0245     NZ2  = ((R2+5.)/DELR);
0246     NR2  = NZ2;
0247     RR1  = -DELR;
0248         RR2  = -DELR;
0249     NY   = ((R1+5.)/DELL);
0250     NX   = 2*NY;
0251     // This calculates T_A(b) for beam 1 and stores it in DEN1[IR1]
0252     for ( int IR1 = 1; IR1 <= NR1; IR1++) {
0253         DEN1[IR1] = 0.;
0254         RR1       = RR1+DELR;
0255         Z1        = -DELR/2;
0256
0257         for ( int IZ1 = 1; IZ1 <= NZ1; IZ1++) {
0258             Z1  = Z1+DELR;
0259             RSQ = RR1*RR1+Z1*Z1;
0260             DEN1[IR1] = DEN1[IR1]+1./(1.+exp((sqrt(RSQ)-R1)/A1));
0261         }
0262
0263         DEN1[IR1] = DEN1[IR1]*2.*DELR;
0264     }
0265
0266     // This calculates T_A(b) for beam 2 and stores it in DEN2[IR2]
0267     for ( int IR2 = 1; IR2 <= NR2; IR2++) {
0268         DEN2[IR2] = 0.;
0269         RR2       = RR2+DELR;
0270         Z2        = -DELR/2;
0271
0272         for ( int IZ2 = 1; IZ2 <= NZ2; IZ2++) {
0273             Z2  = Z2+DELR;
0274             RSQ = RR2*RR2+Z2*Z2;
0275             DEN2[IR2] = DEN2[IR2]+1./(1.+exp((sqrt(RSQ)-R2)/A2));
0276         }
0277
0278         DEN2[IR2] = DEN2[IR2]*2.*DELR;
0279     }
0280
0281     AN1 = 0.;
0282     RR1 = 0.;
0283         RR2 = 0.;
0284
0285     for ( int IR1 =1; IR1 <= NR1; IR1++) {
0286         RR1 = RR1+DELR;
0287         AN1 = AN1+RR1*DEN1[IR1]*DELR*2.*starlightConstants::pi;
0288     }
0289     for ( int IR2 =1; IR2 <= NR2; IR2++) {
0290         RR2 = RR2+DELR;
0291         AN2 = AN2+RR2*DEN2[IR2]*DELR*2.*starlightConstants::pi;
0292     }
0293
0294
0295     //.1 to turn mb into fm^2
0296     //Calculate breakup probability here
0297  L100:
0298     _pHadronBreakup = 0.;
0299     if ( b > 25. ) return _pHadronBreakup;
0300     Y = -.5*DELL;
0301     for ( int IY = 1; IY <= NY; IY++) {
0302         Y = Y+DELL;
0303         X = -DELL*float(NY+1);
0304
0305         for ( int IX = 1; IX <=NX; IX++) {
0306             X = X+DELL;
0307             XB = b-X;
0308             RPU = sqrt(X*X+Y*Y);
0309             IRUP = (RPU/DELR)+1;
0310             RTU  = sqrt(XB*XB+Y*Y);
0311             IRUT = (RTU/DELR)+1;
0312             T1   = DEN2[(int)IRUT]*RHO2/AN2;
0313             T2   = DEN1[(int)IRUP]*RHO1/AN1;
0314             //Eq.6 BCW, Baltz, Chasman, White, Nucl. Inst. & Methods A 417, 1 (1998)
0315             _pHadronBreakup=_pHadronBreakup+2.*T1*(1.-exp(-SIGNN*T2))*DELL*DELL;
0316         }//for(IX)
0317     }//for(IY)
0318
0319     return _pHadronBreakup;
0320 }
0321
0322
0323 //______________________________________________________________________________
0324 double
0325 beamBeamSystem::probabilityOfPhotonBreakup(const double impactparameter, const int mode)
0326 {
0327     static double ee[10001], eee[162], se[10001];
0328
0329     _pPhotonBreakup =0.;   //Might default the probability with a different value?
0330     double b = impactparameter;
0331     int zp = 1;  //What about _beam2? Generic approach?
0332     int ap = 0;
0333
0334     //Was initialized at the start of the function originally, been moved inward.
0335     double pxn=0.;
0336     double p1n=0.;
0337
0338     //Used to be done prior to entering the function. Done properly for assymetric?
0339     double gammatarg = 2.*_beamLorentzGamma*_beamLorentzGamma-1.;
0340     double omaxx =0.;
0341     //This was done prior entering the function as well
0342     if (_beamLorentzGamma > 500.){
0343         omaxx=1.E10;
0344     }
0345     else{
0346         omaxx=1.E7;
0347     }
0348
0349
0350     double e1[23]= {0.,103.,106.,112.,119.,127.,132.,145.,171.,199.,230.,235.,
0351                     254.,280.,300.,320.,330.,333.,373.,390.,420.,426.,440.};
0352     double s1[23]= {0.,12.0,11.5,12.0,12.0,12.0,15.0,17.0,28.0,33.0,
0353                     52.0,60.0,70.0,76.0,85.0,86.0,89.0,89.0,75.0,76.0,69.0,59.0,61.0};
0354     double e2[12]={0.,2000.,3270.,4100.,4810.,6210.,6600.,
0355                    7790.,8400.,9510.,13600.,16400.};
0356     double s2[12]={0.,.1266,.1080,.0805,.1017,.0942,.0844,.0841,.0755,.0827,
0357                    .0626,.0740};
0358     double e3[29]={0.,26.,28.,30.,32.,34.,36.,38.,40.,44.,46.,48.,50.,52.,55.,
0359                    57.,62.,64.,66.,69.,72.,74.,76.,79.,82.,86.,92.,98.,103.};
0360     double s3[29]={0.,30.,21.5,22.5,18.5,17.5,15.,14.5,19.,17.5,16.,14.,
0361                    20.,16.5,17.5,17.,15.5,18.,15.5,15.5,15.,13.5,18.,14.5,15.5,12.5,13.,
0362                    13.,12.};
0363     static double sa[161]={0.,0.,.004,.008,.013,.017,.021,.025,.029,.034,.038,.042,.046,
0364                            .051,.055,.059,.063,.067,.072,.076,.08,.085,.09,.095,.1,.108,.116,
0365                            .124,.132,.14,.152,.164,.176,.188,.2,.22,.24,.26,.28,.3,.32,.34,
0366                            .36,.38,.4,.417,.433,.450,.467,.483,.5,.51,.516,.52,.523,.5245,
0367                            .525,.5242,
0368                            .5214,.518,.512,.505,.495,.482,.469,.456,.442,.428,.414,.4,.386,
0369                            .370,.355,.34,.325,.310,.295,.280,.265,.25,.236,.222,.208,.194,
0370                            .180,.166,
0371                            .152,.138,.124,.11,.101,.095,.09,.085,.08,.076,.072,.069,.066,
0372                            .063,.06,.0575,.055,.0525,.05,.04875,.0475,.04625,.045,.04375,
0373                            .0425,.04125,.04,.03875,.0375,.03625,.035,.03375,.0325,.03125,.03,
0374                            .02925,.0285,.02775,.027,.02625,.0255,.02475,.024,.02325,.0225,
0375                            .02175,.021,.02025,.0195,.01875,.018,.01725,.0165,.01575,.015,
0376                            .01425,.0135,.01275,.012,.01125,.0105,.00975,.009,.00825,.0075,
0377                            .00675,.006,.00525,.0045,.00375,.003,.00225,.0015,.00075,0.};
0378
0379
0380
0381     double sen[161]={0.,0.,.012,.025,.038,.028,.028,.038,.035,.029,.039,.035,
0382                      .038,.032,.038,.041,.041,.049,.055,.061,.072,.076,.070,.067,
0383                      .080,.103,.125,.138,.118,.103,.129,.155,.170,.180,.190,.200,
0384                      .215,.250,.302,.310,.301,.315,.330,.355,.380,.400,.410,.420,
0385                      .438,.456,.474,.492,.510,.533,.556,.578,.6,.62,.63,.638,
0386                      .640,.640,.637,.631,.625,.618,.610,.600,.580,.555,.530,.505,
0387                      .480,.455,.435,.410,.385,.360,.340,.320,.300,.285,.270,.255,
0388                      .240,.225,.210,.180,.165,.150,.140,.132,.124,.116,.108,.100,
0389                      .092,.084,.077,.071,.066,.060,.055,.051,.048,.046,.044,.042,
0390                      .040,.038,.036,.034,.032,.030,.028,.027,.026,.025,.025,.025,
0391                      .024,.024,.024,.024,.024,.023,.023,.023,.023,.023,.022,.022,
0392                      .022,.022,.022,.021,.021,.021,.020,.020,
0393                      .020,.019,.018,.017,.016,.015,.014,.013,.012,.011,.010,.009,
0394                      .008,.007,.006,.005,.004,.003,.002,.001,0.};
0395
0396     // gammay,p gamma,n of Armstrong begin at 265 incr 25
0397
0398
0399     double sigt[160]={0.,.4245,.4870,.5269,.4778,.4066,.3341,.2444,.2245,.2005,
0400                       .1783,.1769,.1869,.1940,.2117,.2226,.2327,.2395,.2646,.2790,.2756,
0401                       .2607,.2447,.2211,.2063,.2137,.2088,.2017,.2050,.2015,.2121,.2175,
0402                       .2152,.1917,.1911,.1747,.1650,.1587,.1622,.1496,.1486,.1438,.1556,
0403                       .1468,.1536,.1544,.1536,.1468,.1535,.1442,.1515,.1559,.1541,.1461,
0404                       .1388,.1565,.1502,.1503,.1454,.1389,.1445,.1425,.1415,.1424,.1432,
0405                       .1486,.1539,.1354,.1480,.1443,.1435,.1491,.1435,.1380,.1317,.1445,
0406                       .1375,.1449,.1359,.1383,.1390,.1361,.1286,.1359,.1395,.1327,.1387,
0407                       .1431,.1403,.1404,.1389,.1410,.1304,.1363,.1241,.1284,.1299,.1325,
0408                       .1343,.1387,.1328,.1444,.1334,.1362,.1302,.1338,.1339,.1304,.1314,
0409                       .1287,.1404,.1383,.1292,.1436,.1280,.1326,.1321,.1268,.1278,.1243,
0410                       .1239,.1271,.1213,.1338,.1287,.1343,.1231,.1317,.1214,.1370,.1232,
0411                       .1301,.1348,.1294,.1278,.1227,.1218,.1198,.1193,.1342,.1323,.1248,
0412                       .1220,.1139,.1271,.1224,.1347,.1249,.1163,.1362,.1236,.1462,.1356,
0413                       .1198,.1419,.1324,.1288,.1336,.1335,.1266};
0414
0415
0416     double sigtn[160]={0.,.3125,.3930,.4401,.4582,.3774,.3329,.2996,.2715,.2165,
0417                        .2297,.1861,.1551,.2020,.2073,.2064,.2193,.2275,.2384,.2150,.2494,
0418                        .2133,.2023,.1969,.1797,.1693,.1642,.1463,.1280,.1555,.1489,.1435,
0419                        .1398,.1573,.1479,.1493,.1417,.1403,.1258,.1354,.1394,.1420,.1364,
0420                        .1325,.1455,.1326,.1397,.1286,.1260,.1314,.1378,.1353,.1264,.1471,
0421                        .1650,.1311,.1261,.1348,.1277,.1518,.1297,.1452,.1453,.1598,.1323,
0422                        .1234,.1212,.1333,.1434,.1380,.1330,.12,.12,.12,.12,.12,.12,.12,.12,
0423                        .12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,
0424                        .12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,
0425                        .12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,
0426                        .12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,
0427                        .12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12};
0428
0429
0430
0431     static int IFIRSTP=0;
0432
0433
0434     double si1=0, g1 =0,   o1=0;
0435     int   ne = 0, ij =0;
0436     double delo=0, omax =0, gk1m=0;
0437     static double scon=0., zcon=0.,o0=0.;
0438
0439
0440     double x=0,y=0,eps=0,eta=0,em=0,exx=0,s=0,ictr=0,pom=0,vec=0,gk1=0;
0441
0442     //  maximum energy for GDR dissocation (in target frame, in MeV)
0443
0444     double omax1n=24.01;
0445
0446     if (IFIRSTP != 0) goto L100;
0447
0448     IFIRSTP=1;
0449
0450
0451     //This is dependenant on gold or lead....Might need to expand
0452     if (zp == 79)
0453         {
0454
0455
0456             ap=197;
0457             si1=540.;
0458             g1=4.75;
0459
0460             // peak and minimum energies for GDR excitation (in MeV)
0461             o1=13.70;
0462             o0=8.1;
0463         }
0464     else
0465         {
0466             zp=82;   //assumed to be lead
0467             ap=208;
0468             si1=640.;
0469             g1=4.05;
0470             o1=13.42;
0471             o0=7.4;
0472             for(int j=1;j<=160;j++)
0473                 {
0474
0475                     sa[j]=sen[j];
0476                 }
0477         }
0478     //Part II of initialization
0479     delo = .05;
0480     //.1 to turn mb into fm^2
0481     scon = .1*g1*g1*si1;
0482     zcon = zp/(gammatarg*( pi)*(
0483                                                     hbarcmev))*zp/(gammatarg*( pi)*
0484                      ( hbarcmev))/137.04;//alpha?
0485
0486     //single neutron from GDR, Veyssiere et al. Nucl. Phys. A159, 561 (1970)
0487     for ( int i = 1; i <= 160; i++) {
0488         eee[i] = o0+.1*(i-1);
0489         sa[i]  = 100.*sa[i];
0490     }
0491     //See Baltz, Rhoades-Brown, and Weneser, Phys. Rev. E 54, 4233 (1996)
0492     //for details of the following photo cross-sections
0493     eee[161]=24.1;
0494     ne=int((25.-o0)/delo)+1;
0495     //GDR any number of neutrons, Veyssiere et al., Nucl. Phys. A159, 561 (1970)
0496     for ( int i = 1; i <= ne; i++ ) {
0497         ee[i] = o0+(i-1)*delo;
0498         //cout<<" ee 1 "<<ee[i]<<"  "<<i<<endl;
0499
0500         se[i] = scon*ee[i]*ee[i]/(((o1*o1-ee[i]*ee[i])*(o1*o1-ee[i]*ee[i]))
0501                                   +ee[i]*ee[i]*g1*g1);
0502     }
0503     ij = ne;   //Risky?
0504     //25-103 MeV, Lepretre, et al., Nucl. Phys. A367, 237 (1981)
0505     for ( int j = 1; j <= 27; j++ ) {
0506         ij = ij+1;
0507         ee[ij] = e3[j];
0508         //cout<<" ee 2 "<<ee[ij]<<"  "<<ij<<endl;
0509
0510         se[ij] = .1*ap*s3[j]/208.;
0511     }
0512     //103-440 MeV, Carlos, et al., Nucl. Phys. A431, 573 (1984)
0513     for ( int j = 1; j <= 22; j++ ) {
0514         ij = ij+1;
0515         ee[ij] = e1[j];
0516         //cout<<" ee 3 "<<ee[ij]<<"  "<<ij<<endl;
0517         se[ij] = .1*ap*s1[j]/208.;
0518     }
0519     //440 MeV-2 GeV Armstrong et al.
0520     for ( int j = 9; j <= 70; j++) {
0521         ij = ij+1;
0522         ee[ij] = ee[ij-1]+25.;
0523         //cout<<" ee 4 "<<ee[ij]<<"  "<<ij<<endl;
0524         se[ij] = .1*(zp*sigt[j]+(ap-zp)*sigtn[j]);
0525     }
0526     //2-16.4 GeV Michalowski; Caldwell
0527     for ( int j = 1; j <= 11; j++) {
0528         ij = ij+1;
0529         ee[ij] = e2[j];
0530         //cout<<" ee 5 "<<ee[ij]<<"   "<<ij<<endl;
0531         se[ij] = .1*ap*s2[j];
0532     }
0533     //Regge paramteres
0534     x = .0677;
0535     y = .129;
0536     eps = .0808;
0537     eta = .4525;
0538     em = .94;
0539     exx = pow(10,.05);
0540
0541     //Regge model for high energy
0542     s = .002*em*ee[ij];
0543     //make sure we reach LHC energies
0544     ictr = 100;
0545     if ( gammatarg > (2.*150.*150.)) ictr = 150;
0546     for ( int j = 1; j <= ictr; j++ ) {
0547         ij = ij+1;
0548         s = s*exx;
0549         ee[ij] = 1000.*.5*(s-em*em)/em;
0550         //cout<<" ee 6 "<<ee[ij]<<"   "<<ij<<endl;
0551         pom = x*pow(s,eps);
0552         vec = y*pow(s,(-eta));
0553         se[ij] = .1*.65*ap*(pom+vec);
0554     }
0555     ee[ij+1] = 99999999999.;
0556     //done with initaliation
0557     //clear counters for 1N, XN
0558  L100:
0559
0560     p1n = 0.;
0561     pxn = 0.;
0562     //start XN calculation
0563     //what's the b-dependent highest energy of interest?
0564
0565     omax = min(omaxx,4.*gammatarg*( hbarcmev)/b);
0566     if ( omax < o0 ) return _pPhotonBreakup;
0567     gk1m = bessel::dbesk1(ee[1]*b/(( hbarcmev)*gammatarg));
0568     int k = 2;
0569  L212:
0570     if (ee[k] < omax ) {
0571         gk1 = bessel::dbesk1(ee[k]*b/(( hbarcmev)*gammatarg));
0572         //Eq. 3 of BCW--NIM in Physics Research A 417 (1998) pp1-8:
0573         pxn=pxn+zcon*(ee[k]-ee[k-1])*.5*(se[k-1]*ee[k-1]*gk1m*gk1m+se[k]*ee[k]*gk1*gk1);
0574         k = k + 1;
0575         gk1m = gk1;
0576         goto L212;
0577     }
0578     //one neutron dissociation
0579     omax = min(omax1n,4.*gammatarg*( hbarcmev)/b);
0580     gk1m = bessel::dbesk1(eee[1]*b/(( hbarcmev)*gammatarg));
0581     k = 2;
0582  L102:
0583     if (eee[k] < omax ) {
0584         gk1 = bessel::dbesk1(eee[k]*b/(( hbarcmev)*gammatarg));
0585         //Like Eq3 but with only the one neutron out GDR photo cross section input
0586         p1n = p1n+zcon*(eee[k]-eee[k-1])*.5*(sa[k-1]*eee[k-1]*gk1m*gk1m+sa[k]*eee[k]*gk1*gk1);
0587         k = k+1;
0588         gk1m = gk1;
0589         goto L102;
0590     }
0591
0592
0593     if (( mode) == 1) _pPhotonBreakup = 1.;
0594     if (( mode) == 2) _pPhotonBreakup = (1-exp(-1*pxn))*(1-exp(-1*pxn));
0595     if (( mode) == 3) _pPhotonBreakup = (p1n*exp(-1*pxn))*(p1n*exp(-1*pxn));
0596     if (( mode) == 4) _pPhotonBreakup = exp(-2*pxn);
0597     if (( mode) == 5) _pPhotonBreakup = 1.;
0598     if (( mode) == 6) _pPhotonBreakup = (1. - exp(-2.*pxn));
0599     if (( mode) == 7) _pPhotonBreakup = 2.*exp(-pxn)*(1.-exp(-pxn));
0600
0601     //cout<<pxn<<" "<<zcon<<" "<<ee[k]<<" "<<se[k-1]<<" "<<gk1m<<"  "<<gk1<<"  "<<k<<"  "<<ee[k+1]<< "  "<<b<< endl;
0602
0603     return _pPhotonBreakup;
0604 }