Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // 26 // 27 // 27 // 28 // 28 // 29 // G4 Physics class: G4QuasiElRatios for N+A e 29 // G4 Physics class: G4QuasiElRatios for N+A elastic cross sections 30 // Created: M.V. Kossov, CERN/ITEP(Moscow), 10 30 // Created: M.V. Kossov, CERN/ITEP(Moscow), 10-OCT-01 31 // The last update: M.V. Kossov, CERN/ITEP (Mo 31 // The last update: M.V. Kossov, CERN/ITEP (Moscow) 15-Oct-06 32 // 32 // 33 // ------------------------------------------- 33 // ---------------------------------------------------------------------- 34 // This class has been extracted from the CHIP 34 // This class has been extracted from the CHIPS model. 35 // All the dependencies on CHIPS classes have 35 // All the dependencies on CHIPS classes have been removed. 36 // Short description: Provides percentage of q 36 // Short description: Provides percentage of quasi-free and quasi-elastic 37 // reactions in the inelastic reactions. 37 // reactions in the inelastic reactions. 38 // ------------------------------------------- 38 // ---------------------------------------------------------------------- 39 39 40 40 41 #include "G4QuasiElRatios.hh" 41 #include "G4QuasiElRatios.hh" 42 #include "G4PhysicalConstants.hh" 42 #include "G4PhysicalConstants.hh" 43 #include "G4SystemOfUnits.hh" 43 #include "G4SystemOfUnits.hh" 44 #include "G4Proton.hh" 44 #include "G4Proton.hh" 45 #include "G4Neutron.hh" 45 #include "G4Neutron.hh" 46 #include "G4Deuteron.hh" 46 #include "G4Deuteron.hh" 47 #include "G4Triton.hh" 47 #include "G4Triton.hh" 48 #include "G4He3.hh" 48 #include "G4He3.hh" 49 #include "G4Alpha.hh" 49 #include "G4Alpha.hh" 50 #include "G4ThreeVector.hh" 50 #include "G4ThreeVector.hh" 51 #include "G4CrossSectionDataSetRegistry.hh" 51 #include "G4CrossSectionDataSetRegistry.hh" 52 #include "G4Pow.hh" 52 #include "G4Pow.hh" 53 #include "G4Log.hh" 53 #include "G4Log.hh" 54 #include "G4Exp.hh" 54 #include "G4Exp.hh" 55 55 56 namespace { 56 namespace { 57 const G4int nps=150; // Number o 57 const G4int nps=150; // Number of steps in the R(s) LinTable 58 const G4int mps=nps+1; // Number o 58 const G4int mps=nps+1; // Number of elements in the R(s) LinTable 59 const G4double sma=150.; // The firs 59 const G4double sma=150.; // The first LinTabEl(s=0)=1., s>sma -> logTab 60 const G4double ds=sma/nps; // Step of 60 const G4double ds=sma/nps; // Step of the linear Table 61 const G4int nls=100; // Number o 61 const G4int nls=100; // Number of steps in the R(lns) logTable 62 const G4int mls=nls+1; // Number o 62 const G4int mls=nls+1; // Number of elements in the R(lns) logTable 63 const G4double lsi=5.; // The min 63 const G4double lsi=5.; // The min ln(s) logTabEl(s=148.4 < sma=150.) 64 const G4double lsa=9.; // The max 64 const G4double lsa=9.; // The max ln(s) logTabEl(s=148.4 - 8103. mb) 65 const G4double mi=G4Exp(lsi);// The min s 65 const G4double mi=G4Exp(lsi);// The min s of logTabEl(~ 148.4 mb) 66 const G4double min_s=G4Exp(lsa);// The max 66 const G4double min_s=G4Exp(lsa);// The max s of logTabEl(~ 8103. mb) 67 const G4double dls=(lsa-lsi)/nls;// Step o 67 const G4double dls=(lsa-lsi)/nls;// Step of the logarithmic Table 68 const G4double edls=G4Exp(dls);// Multipli 68 const G4double edls=G4Exp(dls);// Multiplication step of the logarithmic Table 69 const G4double toler=.01; // The tola 69 const G4double toler=.01; // The tolarence mb defining the same cross-section 70 const G4double C=1.246; 70 const G4double C=1.246; 71 const G4double lmi=3.5; // min of (l 71 const G4double lmi=3.5; // min of (lnP-lmi)^2 parabola 72 const G4double pbe=.0557; // elastic ( 72 const G4double pbe=.0557; // elastic (lnP-lmi)^2 parabola coefficient 73 const G4double pbt=.3; // total (ln 73 const G4double pbt=.3; // total (lnP-lmi)^2 parabola coefficient 74 const G4double pmi=.1; // Below tha 74 const G4double pmi=.1; // Below that fast LE calculation is made 75 const G4double pma=1000.; // Above tha 75 const G4double pma=1000.; // Above that fast HE calculation is made 76 const G4int nlp=300; // Number 76 const G4int nlp=300; // Number of steps in the S(lnp) logTable(5% step) 77 const G4int mlp=nlp+1; // Number 77 const G4int mlp=nlp+1; // Number of elements in the S(lnp) logTable 78 const G4double lpi=-5.; // The min 78 const G4double lpi=-5.; // The min ln(p) logTabEl(p=6.7 MeV/c - 22. TeV/c) 79 const G4double lpa=10.; // The max 79 const G4double lpa=10.; // The max ln(p) logTabEl(p=6.7 MeV/c - 22. TeV/c) 80 const G4double mip=G4Exp(lpi);// The min p 80 const G4double mip=G4Exp(lpi);// The min p of logTabEl(~ 6.7 MeV/c) 81 const G4double map=G4Exp(lpa);// The max p 81 const G4double map=G4Exp(lpa);// The max p of logTabEl(~ 22. TeV) 82 const G4double dlp=(lpa-lpi)/nlp;// Step o 82 const G4double dlp=(lpa-lpi)/nlp;// Step of the logarithmic Table 83 const G4double edlp=G4Exp(dlp);// Multipli 83 const G4double edlp=G4Exp(dlp);// Multiplication step of the logarithmic Table 84 } 84 } 85 85 86 86 87 G4QuasiElRatios::G4QuasiElRatios() 87 G4QuasiElRatios::G4QuasiElRatios() 88 { 88 { 89 vT = new std::vector<G4double*>; 89 vT = new std::vector<G4double*>; 90 vL = new std::vector<G4double*>; 90 vL = new std::vector<G4double*>; 91 vX = new std::vector<std::pair<G4double,G4 91 vX = new std::vector<std::pair<G4double,G4double>*>; 92 92 93 lastSRatio=0.; // The last sig 93 lastSRatio=0.; // The last sigma value for which R was calculated 94 lastRRatio=0.; // The last rat 94 lastRRatio=0.; // The last ratio R which was calculated 95 lastARatio=0; // theLast of ca 95 lastARatio=0; // theLast of calculated A 96 lastHRatio=0.; // theLast of ma 96 lastHRatio=0.; // theLast of max s initialized in the LinTable 97 lastNRatio=0; // theLast of to 97 lastNRatio=0; // theLast of topBin number initialized in LinTable 98 lastMRatio=0.; // theLast of re 98 lastMRatio=0.; // theLast of rel max ln(s) initialized in LogTable 99 lastKRatio=0; // theLast of to 99 lastKRatio=0; // theLast of topBin number initialized in LogTable 100 lastTRatio=0; // theLast of po 100 lastTRatio=0; // theLast of pointer to LinTable in the C++ heap 101 lastLRatio=0; // theLast of po 101 lastLRatio=0; // theLast of pointer to LogTable in the C++ heap 102 lastPtot=0.; // The last mome 102 lastPtot=0.; // The last momentum for which XS was calculated 103 lastHtot=0; // The last proj 103 lastHtot=0; // The last projPDG for which XS was calculated 104 lastFtot=true; // The last nucl 104 lastFtot=true; // The last nucleon for which XS was calculated 105 lastItot=0; // The Last index 105 lastItot=0; // The Last index for which XS was calculated 106 lastMtot=0.; // The Last rel m 106 lastMtot=0.; // The Last rel max ln(p) initialized in LogTable 107 lastKtot=0; // The Last topBin 107 lastKtot=0; // The Last topBin number initialized in LogTable 108 lastXtot = nullptr; // The Last ETPoin 108 lastXtot = nullptr; // The Last ETPointers to LogTable in heap 109 109 110 PCSmanager=(G4ChipsProtonElasticXS*)G4Cros 110 PCSmanager=(G4ChipsProtonElasticXS*)G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4ChipsProtonElasticXS::Default_Name()); 111 111 112 NCSmanager=(G4ChipsNeutronElasticXS*)G4Cro 112 NCSmanager=(G4ChipsNeutronElasticXS*)G4CrossSectionDataSetRegistry::Instance()->GetCrossSectionDataSet(G4ChipsNeutronElasticXS::Default_Name()); 113 } 113 } 114 114 115 G4QuasiElRatios::~G4QuasiElRatios() 115 G4QuasiElRatios::~G4QuasiElRatios() 116 { 116 { 117 std::vector<G4double*>::iterator pos; 117 std::vector<G4double*>::iterator pos; 118 for(pos=vT->begin(); pos<vT->end(); pos++) 118 for(pos=vT->begin(); pos<vT->end(); pos++) 119 { delete [] *pos; } 119 { delete [] *pos; } 120 vT->clear(); 120 vT->clear(); 121 delete vT; 121 delete vT; 122 122 123 for(pos=vL->begin(); pos<vL->end(); pos++) 123 for(pos=vL->begin(); pos<vL->end(); pos++) 124 { delete [] *pos; } 124 { delete [] *pos; } 125 vL->clear(); 125 vL->clear(); 126 delete vL; 126 delete vL; 127 127 128 std::vector<std::pair<G4double,G4double>*> 128 std::vector<std::pair<G4double,G4double>*>::iterator pos2; 129 for(pos2=vX->begin(); pos2<vX->end(); pos2 129 for(pos2=vX->begin(); pos2<vX->end(); pos2++) 130 { delete [] *pos2; } 130 { delete [] *pos2; } 131 vX->clear(); 131 vX->clear(); 132 delete vX; 132 delete vX; 133 } 133 } 134 134 135 // Calculation of pair(QuasiFree/Inelastic,Qua 135 // Calculation of pair(QuasiFree/Inelastic,QuasiElastic/QuasiFree) 136 std::pair<G4double,G4double> G4QuasiElRatios:: 136 std::pair<G4double,G4double> G4QuasiElRatios::GetRatios(G4double pIU, G4int pPDG, 137 137 G4int tgZ, G4int tgN) 138 { 138 { 139 G4double R=0.; 139 G4double R=0.; 140 G4double QF2In=1.; 140 G4double QF2In=1.; // Prototype of QuasiFree/Inel ratio for hN_tot 141 G4int tgA=tgZ+tgN; 141 G4int tgA=tgZ+tgN; 142 if(tgA<2) return std::make_pair(QF2In,R); 142 if(tgA<2) return std::make_pair(QF2In,R); // No quasi-elastic on the only nucleon 143 std::pair<G4double,G4double> ElTot=GetElTo 143 std::pair<G4double,G4double> ElTot=GetElTot(pIU, pPDG, tgZ, tgN); // mean hN El&Tot(IU) 144 //if( ( (pPDG>999 && pIU<227.) || pIU<27.) 144 //if( ( (pPDG>999 && pIU<227.) || pIU<27.) && tgA>1) R=1.; // @@ TMP to accelerate @lowE 145 if(pPDG>999 && pIU<227. && tgZ+tgN>1) R=1. 145 if(pPDG>999 && pIU<227. && tgZ+tgN>1) R=1.; // To accelerate @lowE 146 else if(ElTot.second>0.) 146 else if(ElTot.second>0.) 147 { 147 { 148 R=ElTot.first/ElTot.second; 148 R=ElTot.first/ElTot.second; // El/Total ratio (does not depend on units 149 QF2In=GetQF2IN_Ratio(ElTot.second/mill 149 QF2In=GetQF2IN_Ratio(ElTot.second/millibarn, tgZ+tgN); // QuasiFree/Inelastic ratio 150 } 150 } 151 return std::make_pair(QF2In,R); 151 return std::make_pair(QF2In,R); 152 } 152 } 153 153 154 // Calculatio QasiFree/Inelastic Ratio as a fu 154 // Calculatio QasiFree/Inelastic Ratio as a function of total hN cross-section (mb) and A 155 G4double G4QuasiElRatios::GetQF2IN_Ratio(G4dou 155 G4double G4QuasiElRatios::GetQF2IN_Ratio(G4double m_s, G4int A) 156 { 156 { 157 // LogTable is created only if necessary. 157 // LogTable is created only if necessary. The ratio R(s>8100 mb) = 0 for any nuclei 158 if(m_s<toler || A<2) return 1.; 158 if(m_s<toler || A<2) return 1.; 159 if(m_s>min_s) return 0.; 159 if(m_s>min_s) return 0.; 160 160 161 //if(A>238) 161 //if(A>238) 162 //{ 162 //{ 163 // G4cout<<"-Warning-G4QuasiElRatio::Ge 163 // G4cout<<"-Warning-G4QuasiElRatio::GetQF2IN_Ratio:A="<<A<<">238, return zero"<<G4endl; 164 // return 0.; 164 // return 0.; 165 //} 165 //} 166 G4int nDB=(G4int)vARatio.size(); // A num << 166 G4int nDB=vARatio.size(); // A number of nuclei already initialized in AMDB 167 if(nDB && lastARatio==A && m_s==lastSRatio 167 if(nDB && lastARatio==A && m_s==lastSRatio) return lastRRatio; // VI do not use tolerance 168 G4bool found=false; 168 G4bool found=false; 169 G4int i=-1; 169 G4int i=-1; 170 if(nDB) for (i=0; i<nDB; i++) if(A==vARati 170 if(nDB) for (i=0; i<nDB; i++) if(A==vARatio[i]) // Search for this A in AMDB 171 { 171 { 172 found=true; // 172 found=true; // The A value is found 173 break; 173 break; 174 } 174 } 175 if(!nDB || !found) // C 175 if(!nDB || !found) // Create new line in the AMDB 176 { 176 { 177 lastARatio = A; 177 lastARatio = A; 178 lastTRatio = new G4double[mps]; 178 lastTRatio = new G4double[mps]; // Create the linear Table 179 lastNRatio = static_cast<int>(m_s/ds)+ 179 lastNRatio = static_cast<int>(m_s/ds)+1; // MaxBin to be initialized 180 if(lastNRatio>nps) 180 if(lastNRatio>nps) 181 { 181 { 182 lastNRatio=nps; 182 lastNRatio=nps; 183 lastHRatio=sma; 183 lastHRatio=sma; 184 } 184 } 185 else lastHRatio = lastNRatio*ds; 185 else lastHRatio = lastNRatio*ds; // Calculate max initialized s for LinTab 186 G4double sv=0; 186 G4double sv=0; 187 lastTRatio[0]=1.; 187 lastTRatio[0]=1.; 188 for(G4int j=1; j<=lastNRatio; j++) 188 for(G4int j=1; j<=lastNRatio; j++) // Calculate LogTab values 189 { 189 { 190 sv+=ds; 190 sv+=ds; 191 lastTRatio[j]=CalcQF2IN_Ratio(sv,A 191 lastTRatio[j]=CalcQF2IN_Ratio(sv,A); 192 } 192 } 193 lastLRatio=new G4double[mls]; 193 lastLRatio=new G4double[mls]; // Create the logarithmic Table 194 // Initialize the logarithmic Table << 194 if(m_s>sma) // Initialize the logarithmic Table 195 for(G4int j=0; j<mls; ++j) lastLRatio[j]=0.0 << 196 if(m_s>sma) << 197 { 195 { 198 G4double ls=G4Log(m_s); 196 G4double ls=G4Log(m_s); 199 lastKRatio = static_cast<int>((ls- 197 lastKRatio = static_cast<int>((ls-lsi)/dls)+1; // MaxBin to be initialized in LogTaB 200 if(lastKRatio>nls) 198 if(lastKRatio>nls) 201 { 199 { 202 lastKRatio=nls; 200 lastKRatio=nls; 203 lastMRatio=lsa-lsi; 201 lastMRatio=lsa-lsi; 204 } 202 } 205 else lastMRatio = lastKRatio*dls; 203 else lastMRatio = lastKRatio*dls; // Calculate max initialized ln(s)-lsi for LogTab 206 sv=mi; 204 sv=mi; 207 for(G4int j=0; j<=lastKRatio; j++) 205 for(G4int j=0; j<=lastKRatio; j++) // Calculate LogTab values 208 { 206 { 209 lastLRatio[j]=CalcQF2IN_Ratio( 207 lastLRatio[j]=CalcQF2IN_Ratio(sv,A); 210 if(j!=lastKRatio) sv*=edls; 208 if(j!=lastKRatio) sv*=edls; 211 } 209 } 212 } 210 } 213 else // 211 else // LogTab is not initialized 214 { 212 { 215 lastKRatio = 0; 213 lastKRatio = 0; 216 lastMRatio = 0.; 214 lastMRatio = 0.; 217 } 215 } 218 i++; // 216 i++; // Make a new record to AMDB and position on it 219 vARatio.push_back(lastARatio); 217 vARatio.push_back(lastARatio); 220 vHRatio.push_back(lastHRatio); 218 vHRatio.push_back(lastHRatio); 221 vNRatio.push_back(lastNRatio); 219 vNRatio.push_back(lastNRatio); 222 vMRatio.push_back(lastMRatio); 220 vMRatio.push_back(lastMRatio); 223 vKRatio.push_back(lastKRatio); 221 vKRatio.push_back(lastKRatio); 224 vT->push_back(lastTRatio); 222 vT->push_back(lastTRatio); 225 vL->push_back(lastLRatio); 223 vL->push_back(lastLRatio); 226 } 224 } 227 else // T 225 else // The A value was found in AMDB 228 { 226 { 229 lastARatio=vARatio[i]; 227 lastARatio=vARatio[i]; 230 lastHRatio=vHRatio[i]; 228 lastHRatio=vHRatio[i]; 231 lastNRatio=vNRatio[i]; 229 lastNRatio=vNRatio[i]; 232 lastMRatio=vMRatio[i]; 230 lastMRatio=vMRatio[i]; 233 lastKRatio=vKRatio[i]; 231 lastKRatio=vKRatio[i]; 234 lastTRatio=(*vT)[i]; 232 lastTRatio=(*vT)[i]; 235 lastLRatio=(*vL)[i]; 233 lastLRatio=(*vL)[i]; 236 if(m_s>lastHRatio) 234 if(m_s>lastHRatio) // At least LinTab must be updated 237 { 235 { 238 G4int nextN=lastNRatio+1; 236 G4int nextN=lastNRatio+1; // The next bin to be initialized 239 if(lastNRatio<nps) 237 if(lastNRatio<nps) 240 { 238 { 241 G4double sv=lastHRatio; // bug fix by 239 G4double sv=lastHRatio; // bug fix by WP 242 240 243 lastNRatio = static_cast<int>( 241 lastNRatio = static_cast<int>(m_s/ds)+1;// MaxBin to be initialized 244 if(lastNRatio>nps) 242 if(lastNRatio>nps) 245 { 243 { 246 lastNRatio=nps; 244 lastNRatio=nps; 247 lastHRatio=sma; 245 lastHRatio=sma; 248 } 246 } 249 else lastHRatio = lastNRatio*d 247 else lastHRatio = lastNRatio*ds; // Calculate max initialized s for LinTab 250 248 251 for(G4int j=nextN; j<=lastNRat 249 for(G4int j=nextN; j<=lastNRatio; j++)// Calculate LogTab values 252 { 250 { 253 sv+=ds; 251 sv+=ds; 254 lastTRatio[j]=CalcQF2IN_Ra 252 lastTRatio[j]=CalcQF2IN_Ratio(sv,A); 255 } 253 } 256 } // End of LinTab update 254 } // End of LinTab update 257 if(lastNRatio>=nextN) 255 if(lastNRatio>=nextN) 258 { 256 { 259 vHRatio[i]=lastHRatio; 257 vHRatio[i]=lastHRatio; 260 vNRatio[i]=lastNRatio; 258 vNRatio[i]=lastNRatio; 261 } 259 } 262 G4int nextK=lastKRatio+1; 260 G4int nextK=lastKRatio+1; 263 if(!lastKRatio) nextK=0; 261 if(!lastKRatio) nextK=0; 264 if(m_s>sma && lastKRatio<nls) 262 if(m_s>sma && lastKRatio<nls) // LogTab must be updated 265 { 263 { 266 G4double sv=G4Exp(lastMRatio+l 264 G4double sv=G4Exp(lastMRatio+lsi); // Define starting poit (lastM will be changed) 267 G4double ls=G4Log(m_s); 265 G4double ls=G4Log(m_s); 268 lastKRatio = static_cast<int>( 266 lastKRatio = static_cast<int>((ls-lsi)/dls)+1; // MaxBin to be initialized in LogTaB 269 if(lastKRatio>nls) 267 if(lastKRatio>nls) 270 { 268 { 271 lastKRatio=nls; 269 lastKRatio=nls; 272 lastMRatio=lsa-lsi; 270 lastMRatio=lsa-lsi; 273 } 271 } 274 else lastMRatio = lastKRatio*d 272 else lastMRatio = lastKRatio*dls; // Calculate max initialized ln(s)-lsi for LogTab 275 for(G4int j=nextK; j<=lastKRat 273 for(G4int j=nextK; j<=lastKRatio; j++)// Calculate LogTab values 276 { 274 { 277 sv*=edls; 275 sv*=edls; 278 lastLRatio[j]=CalcQF2IN_Ra 276 lastLRatio[j]=CalcQF2IN_Ratio(sv,A); 279 } 277 } 280 } // End of LogTab update 278 } // End of LogTab update 281 if(lastKRatio>=nextK) 279 if(lastKRatio>=nextK) 282 { 280 { 283 vMRatio[i]=lastMRatio; 281 vMRatio[i]=lastMRatio; 284 vKRatio[i]=lastKRatio; 282 vKRatio[i]=lastKRatio; 285 } 283 } 286 } 284 } 287 } 285 } 288 // Now one can use tabeles to calculate th 286 // Now one can use tabeles to calculate the value 289 if(m_s<sma) // 287 if(m_s<sma) // Use linear table 290 { 288 { 291 G4int n=static_cast<int>(m_s/ds); 289 G4int n=static_cast<int>(m_s/ds); // Low edge number of the bin 292 G4double d=m_s-n*ds; 290 G4double d=m_s-n*ds; // Linear shift 293 G4double v=lastTRatio[n]; 291 G4double v=lastTRatio[n]; // Base 294 lastRRatio=v+d*(lastTRatio[n+1]-v)/ds; 292 lastRRatio=v+d*(lastTRatio[n+1]-v)/ds; // Result 295 } 293 } 296 else // U 294 else // Use log table 297 { 295 { 298 G4double ls=G4Log(m_s)-lsi; // 296 G4double ls=G4Log(m_s)-lsi; // ln(s)-l_min 299 G4int n=static_cast<int>(ls/dls); / 297 G4int n=static_cast<int>(ls/dls); // Low edge number of the bin 300 G4double d=ls-n*dls; / 298 G4double d=ls-n*dls; // Log shift 301 G4double v=lastLRatio[n]; 299 G4double v=lastLRatio[n]; // Base 302 lastRRatio=v+d*(lastLRatio[n+1]-v)/dls 300 lastRRatio=v+d*(lastLRatio[n+1]-v)/dls; // Result 303 } 301 } 304 if(lastRRatio<0.) lastRRatio=0.; 302 if(lastRRatio<0.) lastRRatio=0.; 305 if(lastRRatio>1.) lastRRatio=1.; 303 if(lastRRatio>1.) lastRRatio=1.; 306 return lastRRatio; 304 return lastRRatio; 307 } // End of CalcQF2IN_Ratio 305 } // End of CalcQF2IN_Ratio 308 306 309 // Calculatio QasiFree/Inelastic Ratio as a fu 307 // Calculatio QasiFree/Inelastic Ratio as a function of total hN cross-section and A 310 G4double G4QuasiElRatios::CalcQF2IN_Ratio(G4do 308 G4double G4QuasiElRatios::CalcQF2IN_Ratio(G4double m_s, G4int A) 311 { 309 { 312 G4double s2=m_s*m_s; 310 G4double s2=m_s*m_s; 313 G4double s4=s2*s2; 311 G4double s4=s2*s2; 314 G4double ss=std::sqrt(std::sqrt(m_s)); 312 G4double ss=std::sqrt(std::sqrt(m_s)); 315 G4double P=7.48e-5*s2/(1.+8.77e12/s4/s4/s2 313 G4double P=7.48e-5*s2/(1.+8.77e12/s4/s4/s2); 316 G4double E=.2644+.016/(1.+G4Exp((29.54-m_s 314 G4double E=.2644+.016/(1.+G4Exp((29.54-m_s)/2.49)); 317 G4double F=ss*.1526*G4Exp(-s2*ss*.0000859) 315 G4double F=ss*.1526*G4Exp(-s2*ss*.0000859); 318 return C*G4Exp(-E*G4Pow::GetInstance()->po 316 return C*G4Exp(-E*G4Pow::GetInstance()->powA(G4double(A-1.),F))/G4Pow::GetInstance()->powA(G4double(A),P); 319 } // End of CalcQF2IN_Ratio 317 } // End of CalcQF2IN_Ratio 320 318 321 // Calculatio pair(hN_el,hN_tot) (mb): p in Ge 319 // Calculatio pair(hN_el,hN_tot) (mb): p in GeV/c, index(PDG,F) (see FetchElTot) 322 std::pair<G4double,G4double> G4QuasiElRatios:: 320 std::pair<G4double,G4double> G4QuasiElRatios::CalcElTot(G4double p, G4int I) 323 { 321 { 324 // ---------> Each parameter set can have 322 // ---------> Each parameter set can have not more than nPoints=128 parameters 325 323 326 G4double El=0.; // pr 324 G4double El=0.; // prototype of the elastic hN cross-section 327 G4double To=0.; // pr 325 G4double To=0.; // prototype of the total hN cross-section 328 if(p<=0.) 326 if(p<=0.) 329 { 327 { 330 G4cout<<"-Warning-G4QuasiElRatios::Cal 328 G4cout<<"-Warning-G4QuasiElRatios::CalcElTot: p="<<p<<" is zero or negative"<<G4endl; 331 return std::make_pair(El,To); 329 return std::make_pair(El,To); 332 } 330 } 333 if (!I) // pp 331 if (!I) // pp/nn 334 { 332 { 335 if(p<pmi) 333 if(p<pmi) 336 { 334 { 337 G4double p2=p*p; 335 G4double p2=p*p; 338 El=1./(.00012+p2*.2); 336 El=1./(.00012+p2*.2); 339 To=El; 337 To=El; 340 } 338 } 341 else if(p>pma) 339 else if(p>pma) 342 { 340 { 343 G4double lp=G4Log(p)-lmi; 341 G4double lp=G4Log(p)-lmi; 344 G4double lp2=lp*lp; 342 G4double lp2=lp*lp; 345 El=pbe*lp2+6.72; 343 El=pbe*lp2+6.72; 346 To=pbt*lp2+38.2; 344 To=pbt*lp2+38.2; 347 } 345 } 348 else 346 else 349 { 347 { 350 G4double p2=p*p; 348 G4double p2=p*p; 351 G4double LE=1./(.00012+p2*.2); 349 G4double LE=1./(.00012+p2*.2); 352 G4double lp=G4Log(p)-lmi; 350 G4double lp=G4Log(p)-lmi; 353 G4double lp2=lp*lp; 351 G4double lp2=lp*lp; 354 G4double rp2=1./p2; 352 G4double rp2=1./p2; 355 El=LE+(pbe*lp2+6.72+32.6/p)/(1.+rp 353 El=LE+(pbe*lp2+6.72+32.6/p)/(1.+rp2/p); 356 To=LE+(pbt*lp2+38.2+52.7*rp2)/(1.+ 354 To=LE+(pbt*lp2+38.2+52.7*rp2)/(1.+2.72*rp2*rp2); 357 } 355 } 358 } 356 } 359 else if(I==1) // np 357 else if(I==1) // np/pn 360 { 358 { 361 if(p<pmi) 359 if(p<pmi) 362 { 360 { 363 G4double p2=p*p; 361 G4double p2=p*p; 364 El=1./(.00012+p2*(.051+.1*p2)); 362 El=1./(.00012+p2*(.051+.1*p2)); 365 To=El; 363 To=El; 366 } 364 } 367 else if(p>pma) 365 else if(p>pma) 368 { 366 { 369 G4double lp=G4Log(p)-lmi; 367 G4double lp=G4Log(p)-lmi; 370 G4double lp2=lp*lp; 368 G4double lp2=lp*lp; 371 El=pbe*lp2+6.72; 369 El=pbe*lp2+6.72; 372 To=pbt*lp2+38.2; 370 To=pbt*lp2+38.2; 373 } 371 } 374 else 372 else 375 { 373 { 376 G4double p2=p*p; 374 G4double p2=p*p; 377 G4double LE=1./(.00012+p2*(.051+.1 375 G4double LE=1./(.00012+p2*(.051+.1*p2)); 378 G4double lp=G4Log(p)-lmi; 376 G4double lp=G4Log(p)-lmi; 379 G4double lp2=lp*lp; 377 G4double lp2=lp*lp; 380 G4double rp2=1./p2; 378 G4double rp2=1./p2; 381 El=LE+(pbe*lp2+6.72+30./p)/(1.+.49 379 El=LE+(pbe*lp2+6.72+30./p)/(1.+.49*rp2/p); 382 To=LE+(pbt*lp2+38.2)/(1.+.54*rp2*r 380 To=LE+(pbt*lp2+38.2)/(1.+.54*rp2*rp2); 383 } 381 } 384 } 382 } 385 else if(I==2) // pi 383 else if(I==2) // pimp/pipn 386 { 384 { 387 G4double lp=G4Log(p); 385 G4double lp=G4Log(p); 388 if(p<pmi) 386 if(p<pmi) 389 { 387 { 390 G4double lr=lp+1.27; 388 G4double lr=lp+1.27; 391 El=1.53/(lr*lr+.0676); 389 El=1.53/(lr*lr+.0676); 392 To=El*3; 390 To=El*3; 393 } 391 } 394 else if(p>pma) 392 else if(p>pma) 395 { 393 { 396 G4double ld=lp-lmi; 394 G4double ld=lp-lmi; 397 G4double ld2=ld*ld; 395 G4double ld2=ld*ld; 398 G4double sp=std::sqrt(p); 396 G4double sp=std::sqrt(p); 399 El=pbe*ld2+2.4+7./sp; 397 El=pbe*ld2+2.4+7./sp; 400 To=pbt*ld2+22.3+12./sp; 398 To=pbt*ld2+22.3+12./sp; 401 } 399 } 402 else 400 else 403 { 401 { 404 G4double lr=lp+1.27; 402 G4double lr=lp+1.27; // p1 405 G4double LE=1.53/(lr*lr+.0676); 403 G4double LE=1.53/(lr*lr+.0676); // p2, p3 406 G4double ld=lp-lmi; 404 G4double ld=lp-lmi; // p4 (lmi=3.5) 407 G4double ld2=ld*ld; 405 G4double ld2=ld*ld; 408 G4double p2=p*p; 406 G4double p2=p*p; 409 G4double p4=p2*p2; 407 G4double p4=p2*p2; 410 G4double sp=std::sqrt(p); 408 G4double sp=std::sqrt(p); 411 G4double lm=lp+.36; 409 G4double lm=lp+.36; // p5 412 G4double md=lm*lm+.04; 410 G4double md=lm*lm+.04; // p6 413 G4double lh=lp-.017; 411 G4double lh=lp-.017; // p7 414 G4double hd=lh*lh+.0025; 412 G4double hd=lh*lh+.0025; // p8 415 El=LE+(pbe*ld2+2.4+7./sp)/(1.+.7/p 413 El=LE+(pbe*ld2+2.4+7./sp)/(1.+.7/p4)+.6/md+.05/hd;//p9(pbe=.0557),p10,p11,p12,p13,p14 416 To=LE*3+(pbt*ld2+22.3+12./sp)/(1.+ 414 To=LE*3+(pbt*ld2+22.3+12./sp)/(1.+.4/p4)+1./md+.06/hd; 417 } 415 } 418 } 416 } 419 else if(I==3) // pi 417 else if(I==3) // pipp/pimn 420 { 418 { 421 G4double lp=G4Log(p); 419 G4double lp=G4Log(p); 422 if(p<pmi) 420 if(p<pmi) 423 { 421 { 424 G4double lr=lp+1.27; 422 G4double lr=lp+1.27; 425 G4double lr2=lr*lr; 423 G4double lr2=lr*lr; 426 El=13./(lr2+lr2*lr2+.0676); 424 El=13./(lr2+lr2*lr2+.0676); 427 To=El; 425 To=El; 428 } 426 } 429 else if(p>pma) 427 else if(p>pma) 430 { 428 { 431 G4double ld=lp-lmi; 429 G4double ld=lp-lmi; 432 G4double ld2=ld*ld; 430 G4double ld2=ld*ld; 433 G4double sp=std::sqrt(p); 431 G4double sp=std::sqrt(p); 434 El=pbe*ld2+2.4+6./sp; 432 El=pbe*ld2+2.4+6./sp; 435 To=pbt*ld2+22.3+5./sp; 433 To=pbt*ld2+22.3+5./sp; 436 } 434 } 437 else 435 else 438 { 436 { 439 G4double lr=lp+1.27; 437 G4double lr=lp+1.27; // p1 440 G4double lr2=lr*lr; 438 G4double lr2=lr*lr; 441 G4double LE=13./(lr2+lr2*lr2+.0676 439 G4double LE=13./(lr2+lr2*lr2+.0676); // p2, p3 442 G4double ld=lp-lmi; 440 G4double ld=lp-lmi; // p4 (lmi=3.5) 443 G4double ld2=ld*ld; 441 G4double ld2=ld*ld; 444 G4double p2=p*p; 442 G4double p2=p*p; 445 G4double p4=p2*p2; 443 G4double p4=p2*p2; 446 G4double sp=std::sqrt(p); 444 G4double sp=std::sqrt(p); 447 G4double lm=lp-.32; 445 G4double lm=lp-.32; // p5 448 G4double md=lm*lm+.0576; 446 G4double md=lm*lm+.0576; // p6 449 El=LE+(pbe*ld2+2.4+6./sp)/(1.+3./p 447 El=LE+(pbe*ld2+2.4+6./sp)/(1.+3./p4)+.7/md; // p7(pbe=.0557), p8, p9, p10, p11 450 To=LE+(pbt*ld2+22.3+5./sp)/(1.+1./ 448 To=LE+(pbt*ld2+22.3+5./sp)/(1.+1./p4)+.8/md; 451 } 449 } 452 } 450 } 453 else if(I==4) // Km 451 else if(I==4) // Kmp/Kmn/K0p/K0n 454 { 452 { 455 453 456 if(p<pmi) 454 if(p<pmi) 457 { 455 { 458 G4double psp=p*std::sqrt(p); 456 G4double psp=p*std::sqrt(p); 459 El=5.2/psp; 457 El=5.2/psp; 460 To=14./psp; 458 To=14./psp; 461 } 459 } 462 else if(p>pma) 460 else if(p>pma) 463 { 461 { 464 G4double ld=G4Log(p)-lmi; 462 G4double ld=G4Log(p)-lmi; 465 G4double ld2=ld*ld; 463 G4double ld2=ld*ld; 466 El=pbe*ld2+2.23; 464 El=pbe*ld2+2.23; 467 To=pbt*ld2+19.5; 465 To=pbt*ld2+19.5; 468 } 466 } 469 else 467 else 470 { 468 { 471 G4double ld=G4Log(p)-lmi; 469 G4double ld=G4Log(p)-lmi; 472 G4double ld2=ld*ld; 470 G4double ld2=ld*ld; 473 G4double sp=std::sqrt(p); 471 G4double sp=std::sqrt(p); 474 G4double psp=p*sp; 472 G4double psp=p*sp; 475 G4double p2=p*p; 473 G4double p2=p*p; 476 G4double p4=p2*p2; 474 G4double p4=p2*p2; 477 G4double lm=p-.39; 475 G4double lm=p-.39; 478 G4double md=lm*lm+.000156; 476 G4double md=lm*lm+.000156; 479 G4double lh=p-1.; 477 G4double lh=p-1.; 480 G4double hd=lh*lh+.0156; 478 G4double hd=lh*lh+.0156; 481 El=5.2/psp+(pbe*ld2+2.23)/(1.-.7/s 479 El=5.2/psp+(pbe*ld2+2.23)/(1.-.7/sp+.075/p4)+.004/md+.15/hd; 482 To=14./psp+(pbt*ld2+19.5)/(1.-.21/ 480 To=14./psp+(pbt*ld2+19.5)/(1.-.21/sp+.52/p4)+.006/md+.30/hd; 483 } 481 } 484 } 482 } 485 else if(I==5) // Kp 483 else if(I==5) // Kpp/Kpn/aKp/aKn 486 { 484 { 487 if(p<pmi) 485 if(p<pmi) 488 { 486 { 489 G4double lr=p-.38; 487 G4double lr=p-.38; 490 G4double lm=p-1.; 488 G4double lm=p-1.; 491 G4double md=lm*lm+.372; 489 G4double md=lm*lm+.372; 492 El=.7/(lr*lr+.0676)+2./md; 490 El=.7/(lr*lr+.0676)+2./md; 493 To=El+.6/md; 491 To=El+.6/md; 494 } 492 } 495 else if(p>pma) 493 else if(p>pma) 496 { 494 { 497 G4double ld=G4Log(p)-lmi; 495 G4double ld=G4Log(p)-lmi; 498 G4double ld2=ld*ld; 496 G4double ld2=ld*ld; 499 El=pbe*ld2+2.23; 497 El=pbe*ld2+2.23; 500 To=pbt*ld2+19.5; 498 To=pbt*ld2+19.5; 501 } 499 } 502 else 500 else 503 { 501 { 504 G4double ld=G4Log(p)-lmi; 502 G4double ld=G4Log(p)-lmi; 505 G4double ld2=ld*ld; 503 G4double ld2=ld*ld; 506 G4double lr=p-.38; 504 G4double lr=p-.38; 507 G4double LE=.7/(lr*lr+.0676); 505 G4double LE=.7/(lr*lr+.0676); 508 G4double sp=std::sqrt(p); 506 G4double sp=std::sqrt(p); 509 G4double p2=p*p; 507 G4double p2=p*p; 510 G4double p4=p2*p2; 508 G4double p4=p2*p2; 511 G4double lm=p-1.; 509 G4double lm=p-1.; 512 G4double md=lm*lm+.372; 510 G4double md=lm*lm+.372; 513 El=LE+(pbe*ld2+2.23)/(1.-.7/sp+.1/ 511 El=LE+(pbe*ld2+2.23)/(1.-.7/sp+.1/p4)+2./md; 514 To=LE+(pbt*ld2+19.5)/(1.+.46/sp+1. 512 To=LE+(pbt*ld2+19.5)/(1.+.46/sp+1.6/p4)+2.6/md; 515 } 513 } 516 } 514 } 517 else if(I==6) // hy 515 else if(I==6) // hyperon-N 518 { 516 { 519 if(p<pmi) 517 if(p<pmi) 520 { 518 { 521 G4double p2=p*p; 519 G4double p2=p*p; 522 El=1./(.002+p2*(.12+p2)); 520 El=1./(.002+p2*(.12+p2)); 523 To=El; 521 To=El; 524 } 522 } 525 else if(p>pma) 523 else if(p>pma) 526 { 524 { 527 G4double lp=G4Log(p)-lmi; 525 G4double lp=G4Log(p)-lmi; 528 G4double lp2=lp*lp; 526 G4double lp2=lp*lp; 529 G4double sp=std::sqrt(p); 527 G4double sp=std::sqrt(p); 530 El=(pbe*lp2+6.72)/(1.+2./sp); 528 El=(pbe*lp2+6.72)/(1.+2./sp); 531 To=(pbt*lp2+38.2+900./sp)/(1.+27./ 529 To=(pbt*lp2+38.2+900./sp)/(1.+27./sp); 532 } 530 } 533 else 531 else 534 { 532 { 535 G4double p2=p*p; 533 G4double p2=p*p; 536 G4double LE=1./(.002+p2*(.12+p2)); 534 G4double LE=1./(.002+p2*(.12+p2)); 537 G4double lp=G4Log(p)-lmi; 535 G4double lp=G4Log(p)-lmi; 538 G4double lp2=lp*lp; 536 G4double lp2=lp*lp; 539 G4double p4=p2*p2; 537 G4double p4=p2*p2; 540 G4double sp=std::sqrt(p); 538 G4double sp=std::sqrt(p); 541 El=LE+(pbe*lp2+6.72+99./p2)/(1.+2. 539 El=LE+(pbe*lp2+6.72+99./p2)/(1.+2./sp+2./p4); 542 To=LE+(pbt*lp2+38.2+900./sp)/(1.+2 540 To=LE+(pbt*lp2+38.2+900./sp)/(1.+27./sp+3./p4); 543 } 541 } 544 } 542 } 545 else if(I==7) // an 543 else if(I==7) // antibaryon-N 546 { 544 { 547 if(p>pma) 545 if(p>pma) 548 { 546 { 549 G4double lp=G4Log(p)-lmi; 547 G4double lp=G4Log(p)-lmi; 550 G4double lp2=lp*lp; 548 G4double lp2=lp*lp; 551 El=pbe*lp2+6.72; 549 El=pbe*lp2+6.72; 552 To=pbt*lp2+38.2; 550 To=pbt*lp2+38.2; 553 } 551 } 554 else 552 else 555 { 553 { 556 G4double ye=G4Pow::GetInstance()-> 554 G4double ye=G4Pow::GetInstance()->powA(p,1.25); 557 G4double yt=G4Pow::GetInstance()-> 555 G4double yt=G4Pow::GetInstance()->powA(p,.35); 558 G4double lp=G4Log(p)-lmi; 556 G4double lp=G4Log(p)-lmi; 559 G4double lp2=lp*lp; 557 G4double lp2=lp*lp; 560 El=80./(ye+1.)+pbe*lp2+6.72; 558 El=80./(ye+1.)+pbe*lp2+6.72; 561 To=(80./yt+.3)/yt+pbt*lp2+38.2; 559 To=(80./yt+.3)/yt+pbt*lp2+38.2; 562 } 560 } 563 } 561 } 564 else 562 else 565 { 563 { 566 G4cout<<"*Error*G4QuasiElRatios::CalcE 564 G4cout<<"*Error*G4QuasiElRatios::CalcElTot:ind="<<I<<" is not defined (0-7)"<<G4endl; 567 G4Exception("G4QuasiElRatios::CalcElTo 565 G4Exception("G4QuasiElRatios::CalcElTot:","23",FatalException,"QEcrash"); 568 } 566 } 569 if(El>To) El=To; 567 if(El>To) El=To; 570 return std::make_pair(El,To); 568 return std::make_pair(El,To); 571 } // End of CalcElTot 569 } // End of CalcElTot 572 570 573 // For hadron PDG with momentum Mom (GeV/c) on 571 // For hadron PDG with momentum Mom (GeV/c) on N (p/n) calculate <sig_el,sig_tot> pair (mb) 574 std::pair<G4double,G4double> G4QuasiElRatios:: 572 std::pair<G4double,G4double> G4QuasiElRatios::GetElTotXS(G4double p, G4int PDG, G4bool F) 575 { 573 { 576 G4int ind=0; 574 G4int ind=0; // Prototype of the reaction index 577 G4bool kfl=true; 575 G4bool kfl=true; // Flag of K0/aK0 oscillation 578 G4bool kf=false; 576 G4bool kf=false; 579 if(PDG==130||PDG==310) 577 if(PDG==130||PDG==310) 580 { 578 { 581 kf=true; 579 kf=true; 582 if(G4UniformRand()>.5) kfl=false; 580 if(G4UniformRand()>.5) kfl=false; 583 } 581 } 584 if ( (PDG == 2212 && F) || (PDG == 21 582 if ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) ind=0; // pp/nn 585 else if ( (PDG == 2112 && F) || (PDG == 22 583 else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1; // np/pn 586 else if ( (PDG == -211 && F) || (PDG == 21 584 else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2; // pimp/pipn 587 else if ( (PDG == 211 && F) || (PDG == -21 585 else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3; // pipp/pimn 588 //AR-Jul2020: Extended to charmed and bott 586 //AR-Jul2020: Extended to charmed and bottom hadrons: 589 // - treat mesons with constituent c quark 587 // - treat mesons with constituent c quark or b quark as a meson with s quark (e.g. K-); 590 // - treat mesons with constituent cbar an 588 // - treat mesons with constituent cbar antiquark or bbar antiquark as a meson with sbar antiquark (e.g. K+); 591 // - treat all heavy baryons (i.e. hyperon 589 // - treat all heavy baryons (i.e. hyperons, charmed and bottom baryons) as lambda; 592 // - treat all heavy anti-baryons (i.e. an 590 // - treat all heavy anti-baryons (i.e. anti-hyperons, charmed and bottom anti-baryons) as anti-p/anti-n. 593 else if ( PDG == -321 || PDG == -311 || (k 591 else if ( PDG == -321 || PDG == -311 || (kf && !kfl) || // mesons with s quark 594 PDG == 411 || PDG == 421 || PD 592 PDG == 411 || PDG == 421 || PDG == 431 || // mesons with c quark 595 PDG == -521 || PDG == -511 || PDG == - 593 PDG == -521 || PDG == -511 || PDG == -531 || PDG == -541 ) ind=4; // mesons with b quark 596 else if ( PDG == 321 || PDG == 311 || (k 594 else if ( PDG == 321 || PDG == 311 || (kf && kfl) || // mesons with sbar antiquark 597 PDG == -411 || PDG == -421 || PD 595 PDG == -411 || PDG == -421 || PDG == -431 || // mesons with cbar antiquark 598 PDG == 521 || PDG == 511 || PDG == 596 PDG == 521 || PDG == 511 || PDG == 531 || PDG == 541 ) ind=5; // mesons with bbar antiquark 599 else if ( PDG > 3000 && PDG < 5333 ) ind 597 else if ( PDG > 3000 && PDG < 5333 ) ind=6; // @@ for all heavy baryons - take Lambda 600 else if ( PDG > -5333 && PDG < -2000 ) ind 598 else if ( PDG > -5333 && PDG < -2000 ) ind=7; // @@ for all anti-baryons (anti-p/anti-n) 601 else { 599 else { 602 G4cout<<"*Error*G4QuasiElRatios::CalcE 600 G4cout<<"*Error*G4QuasiElRatios::CalcElTotXS: PDG="<<PDG 603 <<", while it is defined only for p,n, 601 <<", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<G4endl; 604 G4Exception("G4QuasiElRatio::CalcElTot 602 G4Exception("G4QuasiElRatio::CalcElTotXS:","22",FatalException,"QEcrash"); 605 } 603 } 606 return CalcElTot(p,ind); 604 return CalcElTot(p,ind); 607 } 605 } 608 606 609 // Calculatio pair(hN_el,hN_tot)(mb): p in GeV 607 // Calculatio pair(hN_el,hN_tot)(mb): p in GeV/c, F=true -> N=proton, F=false -> N=neutron 610 std::pair<G4double,G4double> G4QuasiElRatios:: 608 std::pair<G4double,G4double> G4QuasiElRatios::FetchElTot(G4double p, G4int PDG, G4bool F) 611 { 609 { 612 // LogTable is created only if necessary. 610 // LogTable is created only if necessary. The ratio R(s>8100 mb) = 0 for any nuclei 613 G4int nDB=(G4int)vItot.size(); // A numbe << 611 G4int nDB=vItot.size(); // A number of hadrons already initialized in AMDB 614 if(nDB && lastHtot==PDG && lastFtot==F && 612 if(nDB && lastHtot==PDG && lastFtot==F && p>0. && p==lastPtot) return lastRtot;// VI don't use toler. 615 // if(nDB && lastH==PDG && lastF==F && p> 613 // if(nDB && lastH==PDG && lastF==F && p>0. && std::fabs(p-lastP)/p<toler) return lastR; 616 lastHtot=PDG; 614 lastHtot=PDG; 617 lastFtot=F; 615 lastFtot=F; 618 G4int ind=-1; // 616 G4int ind=-1; // Prototipe of the index of the PDG/F combination 619 // i=0: pp(nn), i=1: np(pn), i=2: pimp(pip 617 // i=0: pp(nn), i=1: np(pn), i=2: pimp(pipn), i=3: pipp(pimn), i=4: Kmp(Kmn,K0n,K0p), 620 // i=5: Kpp(Kpn,aK0n,aK0p), i=6: Hp(Hn), i 618 // i=5: Kpp(Kpn,aK0n,aK0p), i=6: Hp(Hn), i=7: app(apn,ann,anp) 621 G4bool kfl=true; 619 G4bool kfl=true; // Flag of K0/aK0 oscillation 622 G4bool kf=false; 620 G4bool kf=false; 623 if(PDG==130||PDG==310) 621 if(PDG==130||PDG==310) 624 { 622 { 625 kf=true; 623 kf=true; 626 if(G4UniformRand()>.5) kfl=false; 624 if(G4UniformRand()>.5) kfl=false; 627 } 625 } 628 if ( (PDG == 2212 && F) || (PDG == 21 626 if ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) ind=0; // pp/nn 629 else if ( (PDG == 2112 && F) || (PDG == 22 627 else if ( (PDG == 2112 && F) || (PDG == 2212 && !F) ) ind=1; // np/pn 630 else if ( (PDG == -211 && F) || (PDG == 21 628 else if ( (PDG == -211 && F) || (PDG == 211 && !F) ) ind=2; // pimp/pipn 631 else if ( (PDG == 211 && F) || (PDG == -21 629 else if ( (PDG == 211 && F) || (PDG == -211 && !F) ) ind=3; // pipp/pimn 632 //AR-Jul2020: Extended to charmed and bott 630 //AR-Jul2020: Extended to charmed and bottom hadrons: 633 // - treat mesons with constituent c quark 631 // - treat mesons with constituent c quark or b quark as a meson with s quark (e.g. K-); 634 // - treat mesons with constituent cbar an 632 // - treat mesons with constituent cbar antiquark or bbar antiquark as a meson with sbar antiquark (e.g. K+); 635 // - treat all heavy baryons (i.e. hyperon 633 // - treat all heavy baryons (i.e. hyperons, charmed and bottom baryons) as lambda; 636 // - treat all heavy anti-baryons (i.e. an 634 // - treat all heavy anti-baryons (i.e. anti-hyperons, charmed and bottom anti-baryons) as anti-p/anti-n. 637 else if ( PDG == -321 || PDG == -311 || (k 635 else if ( PDG == -321 || PDG == -311 || (kf && !kfl) || // mesons with s quark 638 PDG == 411 || PDG == 421 || PD 636 PDG == 411 || PDG == 421 || PDG == 431 || // mesons with c quark 639 PDG == -521 || PDG == -511 || PDG == - 637 PDG == -521 || PDG == -511 || PDG == -531 || PDG == -541 ) ind=4; // mesons with b quark 640 else if ( PDG == 321 || PDG == 311 || (k 638 else if ( PDG == 321 || PDG == 311 || (kf && kfl) || // mesons with sbar antiquark 641 PDG == -411 || PDG == -421 || PD 639 PDG == -411 || PDG == -421 || PDG == -431 || // mesons with cbar antiquark 642 PDG == 521 || PDG == 511 || PDG == 640 PDG == 521 || PDG == 511 || PDG == 531 || PDG == 541 ) ind=5; // mesons with bbar antiquark 643 else if ( PDG > 3000 && PDG < 5333 ) ind 641 else if ( PDG > 3000 && PDG < 5333 ) ind=6; // @@ for all heavy baryons - take Lambda 644 else if ( PDG > -5333 && PDG < -2000 ) ind 642 else if ( PDG > -5333 && PDG < -2000 ) ind=7; // @@ for all anti-baryons (anti-p/anti-n) 645 else { 643 else { 646 G4cout<<"*Error*G4QuasiElRatios::Fetch 644 G4cout<<"*Error*G4QuasiElRatios::FetchElTot: PDG="<<PDG 647 <<", while it is defined only for p,n, 645 <<", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<G4endl; 648 G4Exception("G4QuasiELRatio::FetchElTo 646 G4Exception("G4QuasiELRatio::FetchElTot:","22",FatalException,"QECrash"); 649 } 647 } 650 if(nDB && lastItot==ind && p>0. && p==last 648 if(nDB && lastItot==ind && p>0. && p==lastPtot) return lastRtot; // VI do not use toler 651 // if(nDB && lastI==ind && p>0. && std::f 649 // if(nDB && lastI==ind && p>0. && std::fabs(p-lastP)/p<toler) return lastR; 652 if(p<=mip || p>=map) return CalcElTot(p,in 650 if(p<=mip || p>=map) return CalcElTot(p,ind); // @@ Slow calculation ! (Warning?) 653 G4bool found=false; 651 G4bool found=false; 654 G4int i=-1; 652 G4int i=-1; 655 if(nDB) for (i=0; i<nDB; i++) if(ind==vIto 653 if(nDB) for (i=0; i<nDB; i++) if(ind==vItot[i]) // Sirch for this index in AMDB 656 { 654 { 657 found=true; 655 found=true; // The index is found 658 break; 656 break; 659 } 657 } 660 G4double lp=G4Log(p); 658 G4double lp=G4Log(p); 661 if(!nDB || !found) 659 if(!nDB || !found) // Create new line in the AMDB 662 { 660 { 663 lastXtot = new std::pair<G4double,G4do 661 lastXtot = new std::pair<G4double,G4double>[mlp]; // Create logarithmic Table for ElTot 664 lastItot = ind; 662 lastItot = ind; // Remember the initialized inex 665 lastKtot = static_cast<int>((lp-lpi)/d 663 lastKtot = static_cast<int>((lp-lpi)/dlp)+1; // MaxBin to be initialized in LogTaB 666 if(lastKtot>nlp) 664 if(lastKtot>nlp) 667 { 665 { 668 lastKtot=nlp; 666 lastKtot=nlp; 669 lastMtot=lpa-lpi; 667 lastMtot=lpa-lpi; 670 } 668 } 671 else lastMtot = lastKtot*dlp; 669 else lastMtot = lastKtot*dlp; // Calculate max initialized ln(p)-lpi for LogTab 672 G4double pv=mip; 670 G4double pv=mip; 673 for(G4int j=0; j<=lastKtot; j++) 671 for(G4int j=0; j<=lastKtot; j++) // Calculate LogTab values 674 { 672 { 675 lastXtot[j]=CalcElTot(pv,ind); 673 lastXtot[j]=CalcElTot(pv,ind); 676 if(j!=lastKtot) pv*=edlp; 674 if(j!=lastKtot) pv*=edlp; 677 } 675 } 678 i++; / 676 i++; // Make a new record to AMDB and position on it 679 vItot.push_back(lastItot); 677 vItot.push_back(lastItot); 680 vMtot.push_back(lastMtot); 678 vMtot.push_back(lastMtot); 681 vKtot.push_back(lastKtot); 679 vKtot.push_back(lastKtot); 682 vX->push_back(lastXtot); 680 vX->push_back(lastXtot); 683 } 681 } 684 else // 682 else // The A value was found in AMDB 685 { 683 { 686 lastItot=vItot[i]; 684 lastItot=vItot[i]; 687 lastMtot=vMtot[i]; 685 lastMtot=vMtot[i]; 688 lastKtot=vKtot[i]; 686 lastKtot=vKtot[i]; 689 lastXtot=(*vX)[i]; 687 lastXtot=(*vX)[i]; 690 G4int nextK=lastKtot+1; 688 G4int nextK=lastKtot+1; 691 G4double lpM=lastMtot+lpi; 689 G4double lpM=lastMtot+lpi; 692 if(lp>lpM && lastKtot<nlp) 690 if(lp>lpM && lastKtot<nlp) // LogTab must be updated 693 { 691 { 694 lastKtot = static_cast<int>((lp-lp 692 lastKtot = static_cast<int>((lp-lpi)/dlp)+1; // MaxBin to be initialized in LogTab 695 if(lastKtot>nlp) 693 if(lastKtot>nlp) 696 { 694 { 697 lastKtot=nlp; 695 lastKtot=nlp; 698 lastMtot=lpa-lpi; 696 lastMtot=lpa-lpi; 699 } 697 } 700 else lastMtot = lastKtot*dlp; 698 else lastMtot = lastKtot*dlp; // Calculate max initialized ln(p)-lpi for LogTab 701 G4double pv=G4Exp(lpM); // m 699 G4double pv=G4Exp(lpM); // momentum of the last calculated beam 702 for(G4int j=nextK; j<=lastKtot; j+ 700 for(G4int j=nextK; j<=lastKtot; j++)// Calculate LogTab values 703 { 701 { 704 pv*=edlp; 702 pv*=edlp; 705 lastXtot[j]=CalcElTot(pv,ind); 703 lastXtot[j]=CalcElTot(pv,ind); 706 } 704 } 707 } // End of LogTab update 705 } // End of LogTab update 708 if(lastKtot>=nextK) 706 if(lastKtot>=nextK) // The AMDB was apdated 709 { 707 { 710 vMtot[i]=lastMtot; 708 vMtot[i]=lastMtot; 711 vKtot[i]=lastKtot; 709 vKtot[i]=lastKtot; 712 } 710 } 713 } 711 } 714 // Now one can use tabeles to calculate th 712 // Now one can use tabeles to calculate the value 715 G4double dlpp=lp-lpi; 713 G4double dlpp=lp-lpi; // Shifted log(p) value 716 G4int n=static_cast<int>(dlpp/dlp); 714 G4int n=static_cast<int>(dlpp/dlp); // Low edge number of the bin 717 G4double d=dlpp-n*dlp; 715 G4double d=dlpp-n*dlp; // Log shift 718 G4double e=lastXtot[n].first; 716 G4double e=lastXtot[n].first; // E-Base 719 lastRtot.first=e+d*(lastXtot[n+1].first-e) 717 lastRtot.first=e+d*(lastXtot[n+1].first-e)/dlp; // E-Result 720 if(lastRtot.first<0.) lastRtot.first = 0. 718 if(lastRtot.first<0.) lastRtot.first = 0.; 721 G4double t=lastXtot[n].second; 719 G4double t=lastXtot[n].second; // T-Base 722 lastRtot.second=t+d*(lastXtot[n+1].second- 720 lastRtot.second=t+d*(lastXtot[n+1].second-t)/dlp; // T-Result 723 if(lastRtot.second<0.) lastRtot.second= 0. 721 if(lastRtot.second<0.) lastRtot.second= 0.; 724 if(lastRtot.first>lastRtot.second) lastRto 722 if(lastRtot.first>lastRtot.second) lastRtot.first = lastRtot.second; 725 return lastRtot; 723 return lastRtot; 726 } // End of FetchElTot 724 } // End of FetchElTot 727 725 728 // (Mean Elastic and Mean Total) Cross-Section 726 // (Mean Elastic and Mean Total) Cross-Sections (mb) for PDG+(Z,N) at P=p[GeV/c] 729 std::pair<G4double,G4double> G4QuasiElRatios:: 727 std::pair<G4double,G4double> G4QuasiElRatios::GetElTot(G4double pIU, G4int hPDG, 730 728 G4int Z, G4int N) 731 { 729 { 732 G4double pGeV=pIU/gigaelectronvolt; 730 G4double pGeV=pIU/gigaelectronvolt; 733 if(Z<1 && N<1) 731 if(Z<1 && N<1) 734 { 732 { 735 G4cout<<"-Warning-G4QuasiElRatio::GetE 733 G4cout<<"-Warning-G4QuasiElRatio::GetElTot:Z="<<Z<<",N="<<N<<", return zero"<<G4endl; 736 return std::make_pair(0.,0.); 734 return std::make_pair(0.,0.); 737 } 735 } 738 std::pair<G4double,G4double> hp=FetchElTot 736 std::pair<G4double,G4double> hp=FetchElTot(pGeV, hPDG, true); 739 std::pair<G4double,G4double> hn=FetchElTot 737 std::pair<G4double,G4double> hn=FetchElTot(pGeV, hPDG, false); 740 G4double A=(Z+N)/millibarn; 738 G4double A=(Z+N)/millibarn; // To make the result in independent units(IU) 741 return std::make_pair((Z*hp.first+N*hn.fir 739 return std::make_pair((Z*hp.first+N*hn.first)/A,(Z*hp.second+N*hn.second)/A); 742 } // End of GetElTot 740 } // End of GetElTot 743 741 744 // (Mean Elastic and Mean Total) Cross-Section 742 // (Mean Elastic and Mean Total) Cross-Sections (mb) for PDG+(Z,N) at P=p[GeV/c] 745 std::pair<G4double,G4double> G4QuasiElRatios:: 743 std::pair<G4double,G4double> G4QuasiElRatios::GetChExFactor(G4double pIU, G4int hPDG, 746 744 G4int Z, G4int N) 747 { 745 { 748 G4double pGeV=pIU/gigaelectronvolt; 746 G4double pGeV=pIU/gigaelectronvolt; 749 G4double resP=0.; 747 G4double resP=0.; 750 G4double resN=0.; 748 G4double resN=0.; 751 if(Z<1 && N<1) 749 if(Z<1 && N<1) 752 { 750 { 753 G4cout<<"-Warning-G4QuasiElRatio::GetC 751 G4cout<<"-Warning-G4QuasiElRatio::GetChExF:Z="<<Z<<",N="<<N<<", return zero"<<G4endl; 754 return std::make_pair(resP,resN); 752 return std::make_pair(resP,resN); 755 } 753 } 756 G4double A=Z+N; 754 G4double A=Z+N; 757 G4double pf=0.; 755 G4double pf=0.; // Possibility to interact with a proton 758 G4double nf=0.; 756 G4double nf=0.; // Possibility to interact with a neutron 759 if (hPDG==-211||hPDG==-321||hPDG==3112|| 757 if (hPDG==-211||hPDG==-321||hPDG==3112||hPDG==3212||hPDG==3312) pf=Z/(A+N); 760 else if(hPDG==211||hPDG==321||hPDG==3222|| 758 else if(hPDG==211||hPDG==321||hPDG==3222||hPDG==3212||hPDG==3322) nf=N/(A+Z); 761 else if(hPDG==-311||hPDG==311||hPDG==130|| 759 else if(hPDG==-311||hPDG==311||hPDG==130||hPDG==310) 762 { 760 { 763 G4double dA=A+A; 761 G4double dA=A+A; 764 pf=Z/(dA+N+N); 762 pf=Z/(dA+N+N); 765 nf=N/(dA+Z+Z); 763 nf=N/(dA+Z+Z); 766 } 764 } 767 G4double mult=1.; // Factor of increasing 765 G4double mult=1.; // Factor of increasing multiplicity ( ? @@) 768 if(pGeV>.5) 766 if(pGeV>.5) 769 { 767 { 770 mult=1./(1.+G4Log(pGeV+pGeV))/pGeV; 768 mult=1./(1.+G4Log(pGeV+pGeV))/pGeV; 771 if(mult>1.) mult=1.; 769 if(mult>1.) mult=1.; 772 } 770 } 773 if(pf) 771 if(pf) 774 { 772 { 775 std::pair<G4double,G4double> hp=FetchE 773 std::pair<G4double,G4double> hp=FetchElTot(pGeV, hPDG, true); 776 resP=pf*(hp.second/hp.first-1.)*mult; 774 resP=pf*(hp.second/hp.first-1.)*mult; 777 } 775 } 778 if(nf) 776 if(nf) 779 { 777 { 780 std::pair<G4double,G4double> hn=FetchE 778 std::pair<G4double,G4double> hn=FetchElTot(pGeV, hPDG, false); 781 resN=nf*(hn.second/hn.first-1.)*mult; 779 resN=nf*(hn.second/hn.first-1.)*mult; 782 } 780 } 783 return std::make_pair(resP,resN); 781 return std::make_pair(resP,resN); 784 } // End of GetChExFactor 782 } // End of GetChExFactor 785 783 786 // scatter (pPDG,p4M) on a virtual nucleon (NP 784 // scatter (pPDG,p4M) on a virtual nucleon (NPDG,N4M), result: final pair(newN4M,newp4M) 787 // if(newN4M.e()==0.) - below threshold, XS=0, 785 // if(newN4M.e()==0.) - below threshold, XS=0, no scattering of the progectile happened 788 std::pair<G4LorentzVector,G4LorentzVector> G4Q 786 std::pair<G4LorentzVector,G4LorentzVector> G4QuasiElRatios::Scatter(G4int NPDG, 789 787 G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M) 790 { 788 { 791 static const G4double mNeut= G4Neutron::Ne 789 static const G4double mNeut= G4Neutron::Neutron()->GetPDGMass(); 792 static const G4double mProt= G4Proton::Pro 790 static const G4double mProt= G4Proton::Proton()->GetPDGMass(); 793 static const G4double mDeut= G4Deuteron::D 791 static const G4double mDeut= G4Deuteron::Deuteron()->GetPDGMass(); 794 static const G4double mTrit= G4Triton::Tri 792 static const G4double mTrit= G4Triton::Triton()->GetPDGMass(); 795 static const G4double mHel3= G4He3::He3()- 793 static const G4double mHel3= G4He3::He3()->GetPDGMass(); 796 static const G4double mAlph= G4Alpha::Alph 794 static const G4double mAlph= G4Alpha::Alpha()->GetPDGMass(); 797 795 798 G4LorentzVector pr4M=p4M/megaelectronvolt; 796 G4LorentzVector pr4M=p4M/megaelectronvolt; // Convert 4-momenta in MeV (keep p4M) 799 N4M/=megaelectronvolt; 797 N4M/=megaelectronvolt; 800 G4LorentzVector tot4M=N4M+p4M; 798 G4LorentzVector tot4M=N4M+p4M; 801 G4double mT=mNeut; 799 G4double mT=mNeut; 802 G4int Z=0; 800 G4int Z=0; 803 G4int N=1; 801 G4int N=1; 804 if(NPDG==2212||NPDG==90001000) 802 if(NPDG==2212||NPDG==90001000) 805 { 803 { 806 mT=mProt; 804 mT=mProt; 807 Z=1; 805 Z=1; 808 N=0; 806 N=0; 809 } 807 } 810 else if(NPDG==90001001) 808 else if(NPDG==90001001) 811 { 809 { 812 mT=mDeut; 810 mT=mDeut; 813 Z=1; 811 Z=1; 814 N=1; 812 N=1; 815 } 813 } 816 else if(NPDG==90002001) 814 else if(NPDG==90002001) 817 { 815 { 818 mT=mHel3; 816 mT=mHel3; 819 Z=2; 817 Z=2; 820 N=1; 818 N=1; 821 } 819 } 822 else if(NPDG==90001002) 820 else if(NPDG==90001002) 823 { 821 { 824 mT=mTrit; 822 mT=mTrit; 825 Z=1; 823 Z=1; 826 N=2; 824 N=2; 827 } 825 } 828 else if(NPDG==90002002) 826 else if(NPDG==90002002) 829 { 827 { 830 mT=mAlph; 828 mT=mAlph; 831 Z=2; 829 Z=2; 832 N=2; 830 N=2; 833 } 831 } 834 else if(NPDG!=2112&&NPDG!=90000001) 832 else if(NPDG!=2112&&NPDG!=90000001) 835 { 833 { 836 G4cout<<"Error:G4QuasiElRatios::Scatte 834 G4cout<<"Error:G4QuasiElRatios::Scatter:NPDG="<<NPDG<<" is not 2212 or 2112"<<G4endl; 837 G4Exception("G4QuasiElRatios::Scatter: 835 G4Exception("G4QuasiElRatios::Scatter:","21",FatalException,"QEcomplain"); 838 //return std::make_pair(G4LorentzVecto 836 //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);// Use this if not exception 839 } 837 } 840 G4double mT2=mT*mT; 838 G4double mT2=mT*mT; 841 G4double mP2=pr4M.m2(); 839 G4double mP2=pr4M.m2(); 842 G4double E=(tot4M.m2()-mT2-mP2)/(mT+mT); 840 G4double E=(tot4M.m2()-mT2-mP2)/(mT+mT); 843 G4double E2=E*E; 841 G4double E2=E*E; 844 if(E<0. || E2<mP2) 842 if(E<0. || E2<mP2) 845 { 843 { 846 return std::make_pair(G4LorentzVector( 844 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action 847 } 845 } 848 G4double P=std::sqrt(E2-mP2); 846 G4double P=std::sqrt(E2-mP2); // Momentum in pseudo laboratory system 849 // @@ Temporary NN t-dependence for all ha 847 // @@ Temporary NN t-dependence for all hadrons 850 //if(pPDG>3400 || pPDG<-3400) G4cout<<"-Wa 848 //if(pPDG>3400 || pPDG<-3400) G4cout<<"-Warning-G4QE::Scatter: pPDG="<<pPDG<<G4endl; 851 G4int PDG=2212; 849 G4int PDG=2212; // *TMP* instead of pPDG 852 if(pPDG==2112||pPDG==-211||pPDG==-321) PDG 850 if(pPDG==2112||pPDG==-211||pPDG==-321) PDG=2112; // *TMP* instead of pPDG 853 if(!Z && N==1) // Change f 851 if(!Z && N==1) // Change for Quasi-Elastic on neutron 854 { 852 { 855 Z=1; 853 Z=1; 856 N=0; 854 N=0; 857 if (PDG==2212) PDG=2112; 855 if (PDG==2212) PDG=2112; 858 else if(PDG==2112) PDG=2212; 856 else if(PDG==2112) PDG=2212; 859 } 857 } 860 G4double xSec=0.; / 858 G4double xSec=0.; // Prototype of Recalculated Cross Section *TMP* 861 if(PDG==2212) xSec=PCSmanager->GetChipsCro 859 if(PDG==2212) xSec=PCSmanager->GetChipsCrossSection(P, Z, N, PDG); // P CrossSect *TMP* 862 else xSec=NCSmanager->GetChipsCro 860 else xSec=NCSmanager->GetChipsCrossSection(P, Z, N, PDG); // N CrossSect *TMP* 863 // @@ check a possibility to separate p, n 861 // @@ check a possibility to separate p, n, or alpha (!) 864 if(xSec <= 0.) 862 if(xSec <= 0.) // The cross-section iz 0 -> Do Nothing 865 { 863 { 866 return std::make_pair(G4LorentzVector( 864 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); //Do Nothing Action 867 } 865 } 868 G4double mint=0.; / 866 G4double mint=0.; // Prototype of functional rand -t (MeV^2) *TMP* 869 if(PDG==2212) mint=PCSmanager->GetExchange 867 if(PDG==2212) mint=PCSmanager->GetExchangeT(Z,N,PDG);// P functional rand -t(MeV^2) *TMP* 870 else mint=NCSmanager->GetExchange 868 else mint=NCSmanager->GetExchangeT(Z,N,PDG);// N functional rand -t(MeV^2) *TMP* 871 G4double maxt=0.; 869 G4double maxt=0.; // Prototype of max possible -t 872 if(PDG==2212) maxt=PCSmanager->GetHMaxT(); 870 if(PDG==2212) maxt=PCSmanager->GetHMaxT(); // max possible -t 873 else maxt=NCSmanager->GetHMaxT(); 871 else maxt=NCSmanager->GetHMaxT(); // max possible -t 874 G4double cost=1.-(mint+mint)/maxt; // cos( 872 G4double cost=1.-(mint+mint)/maxt; // cos(theta) in CMS 875 if(cost>1. || cost<-1. || !(cost>-1. || co 873 if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.)) 876 { 874 { 877 if (cost>1.) cost=1.; 875 if (cost>1.) cost=1.; 878 else if(cost<-1.) cost=-1.; 876 else if(cost<-1.) cost=-1.; 879 else 877 else 880 { 878 { 881 G4double tm=0.; 879 G4double tm=0.; 882 if(PDG==2212) tm=PCSmanager->GetHM 880 if(PDG==2212) tm=PCSmanager->GetHMaxT(); 883 else tm=NCSmanager->GetHM 881 else tm=NCSmanager->GetHMaxT(); 884 G4cerr<<"G4QuasiFreeRatio::Scat:*N 882 G4cerr<<"G4QuasiFreeRatio::Scat:*NAN* cost="<<cost<<",-t="<<mint<<",tm="<<tm<<G4endl; 885 return std::make_pair(G4LorentzVec 883 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action 886 } 884 } 887 } 885 } 888 G4LorentzVector reco4M=G4LorentzVector(0., 886 G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,mT); // 4mom of the recoil nucleon 889 G4LorentzVector dir4M=tot4M-G4LorentzVecto 887 G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-mT)*.01); 890 if(!RelDecayIn2(tot4M, pr4M, reco4M, dir4M 888 if(!RelDecayIn2(tot4M, pr4M, reco4M, dir4M, cost, cost)) 891 { 889 { 892 G4cerr<<"G4QFR::Scat:t="<<tot4M<<tot4M 890 G4cerr<<"G4QFR::Scat:t="<<tot4M<<tot4M.m()<<",mT="<<mT<<",mP="<<std::sqrt(mP2)<<G4endl; 893 //G4Exception("G4QFR::Scat:","009",Fat 891 //G4Exception("G4QFR::Scat:","009",FatalException,"Decay of ElasticComp"); 894 return std::make_pair(G4LorentzVector( 892 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action 895 } 893 } 896 return std::make_pair(reco4M*megaelectronv 894 return std::make_pair(reco4M*megaelectronvolt,pr4M*megaelectronvolt); // Result 897 } // End of Scatter 895 } // End of Scatter 898 896 899 // scatter (pPDG,p4M) on a virtual nucleon (NP 897 // scatter (pPDG,p4M) on a virtual nucleon (NPDG,N4M), result: final pair(newN4M,newp4M) 900 // if(newN4M.e()==0.) - below threshold, XS=0, 898 // if(newN4M.e()==0.) - below threshold, XS=0, no scattering of the progectile happened 901 // User should himself change the charge (PDG) 899 // User should himself change the charge (PDG) (e.g. pn->np, pi+n->pi0p, pi-p->pi0n etc.) 902 //AR-Jul2020: No need to change this method in 900 //AR-Jul2020: No need to change this method in order to extended quasi-elastic to 903 // charmed and bottom mesons, becau 901 // charmed and bottom mesons, because it is not used anywhere ! 904 std::pair<G4LorentzVector,G4LorentzVector> G4Q 902 std::pair<G4LorentzVector,G4LorentzVector> G4QuasiElRatios::ChExer(G4int NPDG, 905 903 G4LorentzVector N4M, G4int pPDG, G4LorentzVector p4M) 906 { 904 { 907 static const G4double mNeut= G4Neutron::Ne 905 static const G4double mNeut= G4Neutron::Neutron()->GetPDGMass(); 908 static const G4double mProt= G4Proton::Pro 906 static const G4double mProt= G4Proton::Proton()->GetPDGMass(); 909 G4LorentzVector pr4M=p4M/megaelectronvolt; 907 G4LorentzVector pr4M=p4M/megaelectronvolt; // Convert 4-momenta in MeV(keep p4M) 910 N4M/=megaelectronvolt; 908 N4M/=megaelectronvolt; 911 G4LorentzVector tot4M=N4M+p4M; 909 G4LorentzVector tot4M=N4M+p4M; 912 G4int Z=0; 910 G4int Z=0; 913 G4int N=1; 911 G4int N=1; 914 G4int sPDG=0; 912 G4int sPDG=0; // PDG code of the scattered hadron 915 G4double mS=0.; 913 G4double mS=0.; // proto of mass of scattered hadron 916 G4double mT=mProt; 914 G4double mT=mProt; // mass of the recoil nucleon 917 if(NPDG==2212) 915 if(NPDG==2212) 918 { 916 { 919 mT=mNeut; 917 mT=mNeut; 920 Z=1; 918 Z=1; 921 N=0; 919 N=0; 922 if(pPDG==-211) sPDG=111; 920 if(pPDG==-211) sPDG=111; // pi+ -> pi0 923 else if(pPDG==-321) 921 else if(pPDG==-321) 924 { 922 { 925 sPDG=310; 923 sPDG=310; // K+ -> K0S 926 if(G4UniformRand()>.5) sPDG=130; 924 if(G4UniformRand()>.5) sPDG=130; // K+ -> K0L 927 } 925 } 928 else if(pPDG==-311||pPDG==311||pPDG==1 926 else if(pPDG==-311||pPDG==311||pPDG==130||pPDG==310) sPDG=321; // K0 -> K+ (?) 929 else if(pPDG==3112) sPDG=3212; 927 else if(pPDG==3112) sPDG=3212; // Sigma- -> Sigma0 930 else if(pPDG==3212) sPDG=3222; 928 else if(pPDG==3212) sPDG=3222; // Sigma0 -> Sigma+ 931 else if(pPDG==3312) sPDG=3322; 929 else if(pPDG==3312) sPDG=3322; // Xi- -> Xi0 932 } 930 } 933 else if(NPDG==2112) // Default 931 else if(NPDG==2112) // Default 934 { 932 { 935 if(pPDG==211) sPDG=111; 933 if(pPDG==211) sPDG=111; // pi+ -> pi0 936 else if(pPDG==321) 934 else if(pPDG==321) 937 { 935 { 938 sPDG=310; 936 sPDG=310; // K+ -> K0S 939 if(G4UniformRand()>.5) sPDG=130; 937 if(G4UniformRand()>.5) sPDG=130; // K+ -> K0L 940 } 938 } 941 else if(pPDG==-311||pPDG==311||pPDG==1 939 else if(pPDG==-311||pPDG==311||pPDG==130||pPDG==310) sPDG=-321; // K0 -> K- (?) 942 else if(pPDG==3222) sPDG=3212; 940 else if(pPDG==3222) sPDG=3212; // Sigma+ -> Sigma0 943 else if(pPDG==3212) sPDG=3112; 941 else if(pPDG==3212) sPDG=3112; // Sigma0 -> Sigma- 944 else if(pPDG==3322) sPDG=3312; 942 else if(pPDG==3322) sPDG=3312; // Xi0 -> Xi- 945 } 943 } 946 else 944 else 947 { 945 { 948 G4cout<<"Error:G4QuasiElRatios::ChExer 946 G4cout<<"Error:G4QuasiElRatios::ChExer: NPDG="<<NPDG<<" is not 2212 or 2112"<<G4endl; 949 G4Exception("G4QuasiElRatios::ChExer:" 947 G4Exception("G4QuasiElRatios::ChExer:","21",FatalException,"QE complain"); 950 //return std::make_pair(G4LorentzVecto 948 //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);// Use this if not exception 951 } 949 } 952 if(sPDG) mS=mNeut; 950 if(sPDG) mS=mNeut; 953 else 951 else 954 { 952 { 955 G4cout<<"Error:G4QuasiElRatios::ChExer 953 G4cout<<"Error:G4QuasiElRatios::ChExer: BAD pPDG="<<pPDG<<", NPDG="<<NPDG<<G4endl; 956 G4Exception("G4QuasiElRatios::ChExer:" 954 G4Exception("G4QuasiElRatios::ChExer:","21",FatalException,"QE complain"); 957 //return std::make_pair(G4LorentzVecto 955 //return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M);// Use this if not exception 958 } 956 } 959 G4double mT2=mT*mT; 957 G4double mT2=mT*mT; 960 G4double mS2=mS*mS; 958 G4double mS2=mS*mS; 961 G4double E=(tot4M.m2()-mT2-mS2)/(mT+mT); 959 G4double E=(tot4M.m2()-mT2-mS2)/(mT+mT); 962 G4double E2=E*E; 960 G4double E2=E*E; 963 if(E<0. || E2<mS2) 961 if(E<0. || E2<mS2) 964 { 962 { 965 return std::make_pair(G4LorentzVector( 963 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action 966 } 964 } 967 G4double P=std::sqrt(E2-mS2); 965 G4double P=std::sqrt(E2-mS2); // Momentum in pseudo laboratory system 968 // @@ Temporary NN t-dependence for all ha 966 // @@ Temporary NN t-dependence for all hadrons 969 G4int PDG=2212; 967 G4int PDG=2212; // *TMP* instead of pPDG 970 if(pPDG==2112||pPDG==-211||pPDG==-321) PDG 968 if(pPDG==2112||pPDG==-211||pPDG==-321) PDG=2112; // *TMP* instead of pPDG 971 if(!Z && N==1) // Change f 969 if(!Z && N==1) // Change for Quasi-Elastic on neutron 972 { 970 { 973 Z=1; 971 Z=1; 974 N=0; 972 N=0; 975 if (PDG==2212) PDG=2112; 973 if (PDG==2212) PDG=2112; 976 else if(PDG==2112) PDG=2212; 974 else if(PDG==2112) PDG=2212; 977 } 975 } 978 G4double xSec=0.; / 976 G4double xSec=0.; // Prototype of Recalculated Cross Section *TMP* 979 if(PDG==2212) xSec=PCSmanager->GetChipsCro 977 if(PDG==2212) xSec=PCSmanager->GetChipsCrossSection(P, Z, N, PDG); // P CrossSect *TMP* 980 else xSec=NCSmanager->GetChipsCro 978 else xSec=NCSmanager->GetChipsCrossSection(P, Z, N, PDG); // N CrossSect *TMP* 981 // @@ check a possibility to separate p, n 979 // @@ check a possibility to separate p, n, or alpha (!) 982 if(xSec <= 0.) // The cross-section iz 0 - 980 if(xSec <= 0.) // The cross-section iz 0 -> Do Nothing 983 { 981 { 984 return std::make_pair(G4LorentzVector( 982 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); //Do Nothing Action 985 } 983 } 986 G4double mint=0.; / 984 G4double mint=0.; // Prototype of functional rand -t (MeV^2) *TMP* 987 if(PDG==2212) mint=PCSmanager->GetExchange 985 if(PDG==2212) mint=PCSmanager->GetExchangeT(Z,N,PDG);// P functional rand -t(MeV^2) *TMP* 988 else mint=NCSmanager->GetExchange 986 else mint=NCSmanager->GetExchangeT(Z,N,PDG);// N functional rand -t(MeV^2) *TMP* 989 G4double maxt=0.; 987 G4double maxt=0.; // Prototype of max possible -t 990 if(PDG==2212) maxt=PCSmanager->GetHMaxT(); 988 if(PDG==2212) maxt=PCSmanager->GetHMaxT(); // max possible -t 991 else maxt=NCSmanager->GetHMaxT(); 989 else maxt=NCSmanager->GetHMaxT(); // max possible -t 992 G4double cost=1.-mint/maxt; 990 G4double cost=1.-mint/maxt; // cos(theta) in CMS 993 if(cost>1. || cost<-1. || !(cost>-1. || co 991 if(cost>1. || cost<-1. || !(cost>-1. || cost<=1.)) 994 { 992 { 995 if (cost>1.) cost=1.; 993 if (cost>1.) cost=1.; 996 else if(cost<-1.) cost=-1.; 994 else if(cost<-1.) cost=-1.; 997 else 995 else 998 { 996 { 999 G4cerr<<"G4QuasiFreeRatio::ChExer: 997 G4cerr<<"G4QuasiFreeRatio::ChExer:*NAN* c="<<cost<<",t="<<mint<<",tm="<<maxt<<G4endl; 1000 return std::make_pair(G4LorentzVe 998 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action 1001 } 999 } 1002 } 1000 } 1003 G4LorentzVector reco4M=G4LorentzVector(0. 1001 G4LorentzVector reco4M=G4LorentzVector(0.,0.,0.,mT); // 4mom of the recoil nucleon 1004 pr4M=G4LorentzVector(0.,0.,0.,mS); 1002 pr4M=G4LorentzVector(0.,0.,0.,mS); // 4mom of the scattered hadron 1005 G4LorentzVector dir4M=tot4M-G4LorentzVect 1003 G4LorentzVector dir4M=tot4M-G4LorentzVector(0.,0.,0.,(tot4M.e()-mT)*.01); 1006 if(!RelDecayIn2(tot4M, pr4M, reco4M, dir4 1004 if(!RelDecayIn2(tot4M, pr4M, reco4M, dir4M, cost, cost)) 1007 { 1005 { 1008 G4cerr<<"G4QFR::ChEx:t="<<tot4M<<tot4 1006 G4cerr<<"G4QFR::ChEx:t="<<tot4M<<tot4M.m()<<",mT="<<mT<<",mP="<<mS<<G4endl; 1009 //G4Exception("G4QFR::ChExer:","009", 1007 //G4Exception("G4QFR::ChExer:","009",FatalException,"Decay of ElasticComp"); 1010 return std::make_pair(G4LorentzVector 1008 return std::make_pair(G4LorentzVector(0.,0.,0.,0.),p4M); // Do Nothing Action 1011 } 1009 } 1012 return std::make_pair(reco4M*megaelectron 1010 return std::make_pair(reco4M*megaelectronvolt,pr4M*megaelectronvolt); // Result 1013 } // End of ChExer 1011 } // End of ChExer 1014 1012 1015 // Calculate ChEx/El ratio (p is in independe 1013 // Calculate ChEx/El ratio (p is in independent units, (Z,N) is target, pPDG is projectile) 1016 G4double G4QuasiElRatios::ChExElCoef(G4double 1014 G4double G4QuasiElRatios::ChExElCoef(G4double p, G4int Z, G4int N, G4int pPDG) 1017 { 1015 { 1018 1016 1019 p/=MeV; // 1017 p/=MeV; // Converted from independent units 1020 G4double A=Z+N; 1018 G4double A=Z+N; 1021 if(A<1.5) return 0.; 1019 if(A<1.5) return 0.; 1022 G4double Cex=0.; 1020 G4double Cex=0.; 1023 if (pPDG==2212) Cex=N/(A+Z); 1021 if (pPDG==2212) Cex=N/(A+Z); 1024 else if(pPDG==2112) Cex=Z/(A+N); 1022 else if(pPDG==2112) Cex=Z/(A+N); 1025 else G4cout<<"*Warning*G4CohChrgExchange: 1023 else G4cout<<"*Warning*G4CohChrgExchange::ChExElCoef: wrong PDG="<<pPDG<<G4endl; 1026 Cex*=Cex; // Cohe 1024 Cex*=Cex; // Coherent processes squares the amplitude 1027 // @@ This is true only for nucleons: oth 1025 // @@ This is true only for nucleons: other projectiles must be treated differently 1028 G4double sp=std::sqrt(p); 1026 G4double sp=std::sqrt(p); 1029 G4double p2=p*p; 1027 G4double p2=p*p; 1030 G4double p4=p2*p2; 1028 G4double p4=p2*p2; 1031 G4double dl1=G4Log(p)-5.; 1029 G4double dl1=G4Log(p)-5.; 1032 G4double T=(6.75+.14*dl1*dl1+13./p)/(1.+. 1030 G4double T=(6.75+.14*dl1*dl1+13./p)/(1.+.14/p4)+.6/(p4+.00013); 1033 G4double U=(6.25+8.33e-5/p4/p)*(p*sp+.34) 1031 G4double U=(6.25+8.33e-5/p4/p)*(p*sp+.34)/p2/p; 1034 G4double R=U/T; 1032 G4double R=U/T; 1035 return Cex*R*R; 1033 return Cex*R*R; 1036 } 1034 } 1037 1035 1038 // Decay of Hadron In2Particles f&s, f is in 1036 // Decay of Hadron In2Particles f&s, f is in respect to the direction of HadronMomentumDir 1039 G4bool G4QuasiElRatios::RelDecayIn2(G4Lorentz 1037 G4bool G4QuasiElRatios::RelDecayIn2(G4LorentzVector& theMomentum, G4LorentzVector& f4Mom, G4LorentzVector& s4Mom, 1040 G4Lorent 1038 G4LorentzVector& dir, G4double maxCost, G4double minCost) 1041 { 1039 { 1042 G4double fM2 = f4Mom.m2(); 1040 G4double fM2 = f4Mom.m2(); 1043 G4double fM = std::sqrt(fM2); 1041 G4double fM = std::sqrt(fM2); // Mass of the 1st Hadron 1044 G4double sM2 = s4Mom.m2(); 1042 G4double sM2 = s4Mom.m2(); 1045 G4double sM = std::sqrt(sM2); 1043 G4double sM = std::sqrt(sM2); // Mass of the 2nd Hadron 1046 G4double iM2 = theMomentum.m2(); 1044 G4double iM2 = theMomentum.m2(); 1047 G4double iM = std::sqrt(iM2); 1045 G4double iM = std::sqrt(iM2); // Mass of the decaying hadron 1048 G4double vP = theMomentum.rho(); // 1046 G4double vP = theMomentum.rho(); // Momentum of the decaying hadron 1049 G4double dE = theMomentum.e(); // 1047 G4double dE = theMomentum.e(); // Energy of the decaying hadron 1050 if(dE<vP) 1048 if(dE<vP) 1051 { 1049 { 1052 G4cerr<<"***G4QHad::RelDecIn2: Tachio 1050 G4cerr<<"***G4QHad::RelDecIn2: Tachionic 4-mom="<<theMomentum<<", E-p="<<dE-vP<<G4endl; 1053 G4double accuracy=.000001*vP; 1051 G4double accuracy=.000001*vP; 1054 G4double emodif=std::fabs(dE-vP); 1052 G4double emodif=std::fabs(dE-vP); 1055 //if(emodif<accuracy) 1053 //if(emodif<accuracy) 1056 //{ 1054 //{ 1057 G4cerr<<"G4QHadron::RelDecIn2: *Boost 1055 G4cerr<<"G4QHadron::RelDecIn2: *Boost* E-p shift is corrected to "<<emodif<<G4endl; 1058 theMomentum.setE(vP+emodif+.01*accura 1056 theMomentum.setE(vP+emodif+.01*accuracy); 1059 //} 1057 //} 1060 } 1058 } 1061 G4ThreeVector ltb = theMomentum.boostVect 1059 G4ThreeVector ltb = theMomentum.boostVector();// Boost vector for backward Lorentz Trans. 1062 G4ThreeVector ltf = -ltb; // 1060 G4ThreeVector ltf = -ltb; // Boost vector for forward Lorentz Trans. 1063 G4LorentzVector cdir = dir; // 1061 G4LorentzVector cdir = dir; // A copy to make a transformation to CMS 1064 cdir.boost(ltf); // 1062 cdir.boost(ltf); // Direction transpormed to CMS of the Momentum 1065 G4ThreeVector vdir = cdir.vect(); // 1063 G4ThreeVector vdir = cdir.vect(); // 3-Vector of the direction-particle 1066 G4ThreeVector vx(0.,0.,1.); // 1064 G4ThreeVector vx(0.,0.,1.); // Ort in the direction of the reference particle 1067 G4ThreeVector vy(0.,1.,0.); // 1065 G4ThreeVector vy(0.,1.,0.); // First ort orthogonal to the direction 1068 G4ThreeVector vz(1.,0.,0.); // 1066 G4ThreeVector vz(1.,0.,0.); // Second ort orthoganal to the direction 1069 if(vdir.mag2() > 0.) // 1067 if(vdir.mag2() > 0.) // the refference particle isn't at rest in CMS 1070 { 1068 { 1071 vx = vdir.unit(); 1069 vx = vdir.unit(); // Ort in the direction of the reference particle 1072 G4ThreeVector vv= vx.orthogonal(); 1070 G4ThreeVector vv= vx.orthogonal(); // Not normed orthogonal vector (!) 1073 vy = vv.unit(); 1071 vy = vv.unit(); // First ort orthogonal to the direction 1074 vz = vx.cross(vy); 1072 vz = vx.cross(vy); // Second ort orthoganal to the direction 1075 } 1073 } 1076 if(maxCost> 1.) maxCost= 1.; 1074 if(maxCost> 1.) maxCost= 1.; 1077 if(minCost<-1.) minCost=-1.; 1075 if(minCost<-1.) minCost=-1.; 1078 if(maxCost<-1.) maxCost=-1.; 1076 if(maxCost<-1.) maxCost=-1.; 1079 if(minCost> 1.) minCost= 1.; 1077 if(minCost> 1.) minCost= 1.; 1080 if(minCost> maxCost) minCost=maxCost; 1078 if(minCost> maxCost) minCost=maxCost; 1081 if(std::fabs(iM-fM-sM)<.00000001) 1079 if(std::fabs(iM-fM-sM)<.00000001) 1082 { 1080 { 1083 G4double fR=fM/iM; 1081 G4double fR=fM/iM; 1084 G4double sR=sM/iM; 1082 G4double sR=sM/iM; 1085 f4Mom=fR*theMomentum; 1083 f4Mom=fR*theMomentum; 1086 s4Mom=sR*theMomentum; 1084 s4Mom=sR*theMomentum; 1087 return true; 1085 return true; 1088 } 1086 } 1089 else if (iM+.001<fM+sM || iM==0.) 1087 else if (iM+.001<fM+sM || iM==0.) 1090 {//@@ Later on make a quark content check 1088 {//@@ Later on make a quark content check for the decay 1091 G4cerr<<"***G4QH::RelDecIn2: fM="<<fM 1089 G4cerr<<"***G4QH::RelDecIn2: fM="<<fM<<"+sM="<<sM<<">iM="<<iM<<",d="<<iM-fM-sM<<G4endl; 1092 return false; 1090 return false; 1093 } 1091 } 1094 G4double d2 = iM2-fM2-sM2; 1092 G4double d2 = iM2-fM2-sM2; 1095 G4double p2 = (d2*d2/4.-fM2*sM2)/iM2; 1093 G4double p2 = (d2*d2/4.-fM2*sM2)/iM2; // Decay momentum(^2) in CMS of Quasmon 1096 if(p2<0.) 1094 if(p2<0.) 1097 { 1095 { 1098 p2=0.; 1096 p2=0.; 1099 } 1097 } 1100 G4double p = std::sqrt(p2); 1098 G4double p = std::sqrt(p2); 1101 G4double ct = maxCost; 1099 G4double ct = maxCost; 1102 if(maxCost>minCost) 1100 if(maxCost>minCost) 1103 { 1101 { 1104 G4double dcost=maxCost-minCost; 1102 G4double dcost=maxCost-minCost; 1105 ct = minCost+dcost*G4UniformRand(); 1103 ct = minCost+dcost*G4UniformRand(); 1106 } 1104 } 1107 G4double phi= twopi*G4UniformRand(); // 1105 G4double phi= twopi*G4UniformRand(); // @@ Change 360.*deg to M_TWOPI (?) 1108 G4double pss=0.; 1106 G4double pss=0.; 1109 if(std::fabs(ct)<1.) pss = p * std::sqrt( 1107 if(std::fabs(ct)<1.) pss = p * std::sqrt(1.-ct*ct); 1110 else 1108 else 1111 { 1109 { 1112 if(ct>1.) ct=1.; 1110 if(ct>1.) ct=1.; 1113 if(ct<-1.) ct=-1.; 1111 if(ct<-1.) ct=-1.; 1114 } 1112 } 1115 G4ThreeVector pVect=(pss*std::sin(phi))*v 1113 G4ThreeVector pVect=(pss*std::sin(phi))*vz+(pss*std::cos(phi))*vy+p*ct*vx; 1116 1114 1117 f4Mom.setVect(pVect); 1115 f4Mom.setVect(pVect); 1118 f4Mom.setE(std::sqrt(fM2+p2)); 1116 f4Mom.setE(std::sqrt(fM2+p2)); 1119 s4Mom.setVect((-1)*pVect); 1117 s4Mom.setVect((-1)*pVect); 1120 s4Mom.setE(std::sqrt(sM2+p2)); 1118 s4Mom.setE(std::sqrt(sM2+p2)); 1121 1119 1122 if(f4Mom.e()+.001<f4Mom.rho())G4cerr<<"*G 1120 if(f4Mom.e()+.001<f4Mom.rho())G4cerr<<"*G4QH::RDIn2:*Boost* f4M="<<f4Mom<<",e-p=" 1123 <<f4Mom.e()-f4Mom.rho()<<G4endl; 1121 <<f4Mom.e()-f4Mom.rho()<<G4endl; 1124 f4Mom.boost(ltb); 1122 f4Mom.boost(ltb); // Lor.Trans. of 1st hadron back to LS 1125 if(s4Mom.e()+.001<s4Mom.rho())G4cerr<<"*G 1123 if(s4Mom.e()+.001<s4Mom.rho())G4cerr<<"*G4QH::RDIn2:*Boost* s4M="<<s4Mom<<",e-p=" 1126 <<s4Mom.e()-s4Mom.rho()<<G4endl; 1124 <<s4Mom.e()-s4Mom.rho()<<G4endl; 1127 s4Mom.boost(ltb); 1125 s4Mom.boost(ltb); // Lor.Trans. of 2nd hadron back to LS 1128 return true; 1126 return true; 1129 } // End of "RelDecayIn2" 1127 } // End of "RelDecayIn2" 1130 1128 1131 1129 1132 1130 1133 1131 1134 1132 1135 1133 1136 1134