source: trunk/source/processes/hadronic/models/leading_particle/src/G4Mars5GeV.cc@ 906

Last change on this file since 906 was 819, checked in by garnier, 17 years ago

import all except CVS
File size: 26.5 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// $Id: G4Mars5GeV.cc,v 1.14 2006/06/29 20:43:24 gunter Exp $
28// GEANT4 tag $Name: $
29//
30//
31// ------------------------------------------------------------
32// GEANT 4 class header file
33//
34//
35// ------------------------------------------------------------
36// First Implemention 17 Nov. 1998 M.Asai, H.Kurahige
37//
38// History:
39// modified as hadronic model 28 Oct 2001 N.Kanaya
40// ------------------------------------------------------------
41// This is a Event Biasing mechanism based on MARS code
42// This model is applicable to
43// proton/neutron/pi+-/K+-/gamma/anti_proton
44// with energy < 5.0GeV
45//
46// Original code is MARS13 written by Nikolai Mokhov (FNAL)
47//**************************************************************
48//* MARS13: 9. hA EVENT GENERATOR:
49//* Copyright Nikolai Mokhov (Fermilab)
50//*
51//* LAST CHANGE: 14-NOV-1998
52//**************************************************************
53//* Copyright Nikolai Mokhov (Fermilab)
54//*
55//* MARS13(98)
56//*
57//* INCLUSIVE HADRON(photon)-NUCLEUS VERTEX AT E < 5 GEV !!!
58//* THREE WEIGHTED HADRONS IN FINAL STATE: !!!
59//* IP+A -> N/P(CASC)+ PI+/PI-(K+/K-) + PI0 *//
60
61#include "globals.hh"
62#include "G4Mars5GeV.hh"
63#include <iostream>
64
65G4Mars5GeV::G4Mars5GeV() : G4InelasticInteraction(),
66 maxWeight(1000.0),
67 minWeight(perMillion)
68{
69 std::cout << " MARS13(98)"<<std::endl;
70 std::cout << std::endl;
71 std::cout << " INCLUSIVE HADRON(photon)-NUCLEUS VERTEX AT E < 5 GEV !!! "<<std::endl;
72 std::cout << " THREE WEIGHTED HADRONS IN FINAL STATE: !!!"<<std::endl;
73 std::cout << " IP+A -> N/P(CASC)+ PI+/PI-(K+/K-) + PI0 "<<std::endl;
74 std::cout << " *********************************************************"<<std::endl;
75 std::cout << " Important notice! "<< std::endl;
76 std::cout << " Since 1998 MARS codes used CEM (Cascade-Exciton Model) " << std::endl;
77 std::cout << " for nuclear interactions below 5 GeV " << std::endl;
78 std::cout << " and do NOT use this inclusive model. " << std::endl;
79
80 std::cout << std::endl;
81
82 SetMinEnergy( 1.0*MeV );
83 SetMaxEnergy( 5.0*GeV );
84
85 theParticleTable = G4ParticleTable::GetParticleTable();
86 G4ParticleDefinition* pProton =
87 theParticleTable->FindParticle("proton");
88 if(pProton) ProtonMass = pProton->GetPDGMass();
89
90 // set some constants
91 selec3.Eth = 0.001*MeV;
92}
93
94G4HadFinalState* G4Mars5GeV::ApplyYourself(const G4HadProjectile& aTrack,
95 G4Nucleus& aTarget
96 )
97{
98 theParticleChange.Clear();
99#ifdef G4VERBOSE
100 if (GetVerboseLevel() > 2) {
101 G4cout << " G4Mars5GeV:ApplyYourself" << G4endl;
102 }
103#endif
104
105 // get the incident particle type
106 incidentParticle = &aTrack;
107 // get the incident particle energy/momentum
108 incidentMarsEncoding = GetMarsEncoding(incidentParticle->GetDefinition());
109
110 // Atomic and charge number
111 //GetTargetNuclei( aTrack.GetMaterial() );
112 fANucl = aTarget.GetN();
113 fZNucl = aTarget.GetZ();
114
115 // initialize secondary information
116 numberOfSecondaries = 0;
117 secondaries.Initialize(FastVectorSize);
118
119 G4int idx;
120
121 // invoke MARS
122 Treem5();
123
124 // create secondaries
125 // set max. number of secondaries
126
127 for (idx=0; idx<numberOfSecondaries; idx+=1){
128
129 // check secondary weight
130 G4double fweight = weightOfSecondaries[idx];
131 if (fweight > maxWeight) {
132#ifdef G4VERBOSE
133 if (GetVerboseLevel() > 0) {
134 G4cout << "G4Mars5GeV::ApplyYourself : too big secondary weight ";
135 G4cout << " Weight = " << fweight << G4endl;
136 secondaries[idx]->DumpInfo();
137 }
138#endif
139 } else if (fweight < minWeight) {
140 // track with small weight is neglected
141#ifdef G4VERBOSE
142 if (GetVerboseLevel() > 2) {
143 G4cout << "G4Mars5GeV::ApplyYourself : too small secondary weight ";
144 G4cout << " Weight = " << fweight << G4endl;
145 secondaries[idx]->DumpInfo();
146 }
147#endif
148 } else {
149 G4HadSecondary *track = new G4HadSecondary(secondaries[idx], fweight);
150// Remain unchanged, because this is a member function of G4HadFinalSate
151 theParticleChange.AddSecondary(track);
152 }
153 }
154 // kill incident particle
155 theParticleChange.SetStatusChange(stopAndKill);
156
157 return &theParticleChange;
158
159}
160
161
162void G4Mars5GeV::GetTargetNuclei(const G4Material* material)
163{
164 // get elements in the actual material,
165 const G4ElementVector* theElementVector = material->GetElementVector();
166 const G4double* theAtomicNumDensityVector = material->GetAtomicNumDensityVector();
167 const G4int numberOfElements = material->GetNumberOfElements() ;
168#ifdef G4VERBOSE
169 if (GetVerboseLevel() > 2) {
170 G4cout << " G4Mars5GeV::GetTargetNuclei" << G4endl;
171 }
172#endif
173
174 fANucl = 0.0;
175 fZNucl = 0.0;
176 G4double totNumAtoms = 0.0;
177 for (G4int iel=0; iel < numberOfElements; iel +=1) {
178 totNumAtoms += theAtomicNumDensityVector[iel];
179
180 fZNucl +=
181 theAtomicNumDensityVector[iel]*((*theElementVector)[iel]->GetZ());
182 fANucl +=
183 theAtomicNumDensityVector[iel]*((*theElementVector)[iel]->GetN());
184
185//#ifdef G4VERBOSE
186// if (GetVerboseLevel() > 2) {
187 G4cout << iel << ": " << theAtomicNumDensityVector[iel];
188 G4cout << " Z=" << (*theElementVector)[iel]->GetZ() << " A=" <<
189 (*theElementVector)[iel]->GetN();
190 G4cout << G4endl;
191// }
192//#endif
193 }
194 fANucl /= totNumAtoms;
195 fZNucl /= totNumAtoms;
196#ifdef G4VERBOSE
197 if (GetVerboseLevel() > 2) {
198 G4cout << "<Z>=" << fZNucl;
199 G4cout << "<A>=" << fANucl;
200 G4cout << G4endl;
201 }
202#endif
203}
204
205
206void G4Mars5GeV::Treem5()
207{
208
209 // G4double pMass = incidentParticle->GetDefinition()->GetPDGMass();
210 G4double pE = incidentParticle->GetKineticEnergy();
211 G4int pType = incidentMarsEncoding;
212
213#ifdef G4VERBOSE
214 if (GetVerboseLevel() > 2) {
215 G4cout << " G4Mars5GeV::Treem5() ";
216 G4cout << " Incident Particle: " << incidentParticle->GetDefinition()->GetParticleName();
217 G4cout << " : energy = " << pE/GeV << "[GeV]" << G4endl;
218 }
219#endif
220
221 // CoulombBarrier
222 if (CoulombBarrier(pType, pE)) return;
223
224 G4int ib;
225 if (pType==MarsAP) {
226 ib = MarsP;
227 } else if (pType==MarsGAM){
228 if ( G4UniformRand() >0.5) {
229 ib = MarsPIplus;
230 } else {
231 ib = MarsPIminus;
232 }
233 } else {
234 ib = pType;
235 }
236
237 selec1.Einc = pE;
238 if (pE < 0.5*MeV) pE = 0.5*MeV;
239 selec3.Emax = pE;
240 selec3.X = 0.0;
241 selec3.Pt = 0.0;
242 selec3.P = 0.0;
243
244 // Nucleons at E < 5GeV
245 CreateNucleon(ib, pType, pE);
246
247 // Pion+- or Kaon+- at E < 5GeV
248 CreatePion(ib, pType, pE);
249
250 // Pi0 at E < 5GeV
251 CreatePionZero(ib, pType, pE);
252}
253
254G4bool G4Mars5GeV::CoulombBarrier(G4int pType, G4double pE){
255 static const G4double EthCoulombBarrier = 20.0* MeV;
256 static const G4double AvCoulomb = 1.11*MeV;
257 static const G4double RCoulombTh = 1.0e-5;
258
259 // CoulombBarrier
260 if ( ( pType == MarsP) || ( pType ==MarsPIplus) || ( pType ==MarsKplus) ) {
261 if ( ( pE < EthCoulombBarrier ) && (fANucl >=1.5) ) {
262 G4double pMass = GetParticleDefinition(pType)->GetPDGMass();
263 G4double vCoulomb = AvCoulomb*std::pow(fZNucl/fANucl, 1./3.);
264 G4double tc = pE*(fANucl*ProtonMass)/(pMass+(fANucl*ProtonMass));
265 G4double rCoulomb = 1.0-vCoulomb/tc;
266 if ( rCoulomb < RCoulombTh ) {
267#ifdef G4VERBOSE
268 if (GetVerboseLevel() > 2) {
269 G4cout << " Can not interact because of Coulomb Barrier " << G4endl;
270 }
271#endif
272 return true;
273 }
274 }
275 }
276 return false;
277}
278
279void G4Mars5GeV::CreateNucleon(G4int ib, G4int pType, G4double )
280{
281#ifdef G4VERBOSE
282 if (GetVerboseLevel() > 2) {
283 G4cout << " G4Mars5GeV::CreateNucleon()" << G4endl;
284 }
285#endif
286 if ( pType == MarsGAM) {
287 selec1.Treac = MarsPIplus;
288 selec1.Tprod = MarsN;
289 selec1.V10 = 2.5;
290 } else {
291 if ( ib == MarsP ) {
292 selec1.Treac = MarsP;
293 } else if ( ib == MarsN ) {
294 selec1.Treac = MarsN;
295 } else if ( ib == MarsPIplus ) {
296 selec1.Treac = MarsPIplus;
297 } else if ( ib == MarsPIminus ) {
298 selec1.Treac = MarsPIminus;
299 } else if ( ib == MarsKplus ) {
300 selec1.Treac = MarsPIplus;
301 } else if ( ib == MarsKminus ) {
302 selec1.Treac = MarsPIminus;
303 } else {
304 selec1.Treac = MarsPIminus;
305 }
306 if (G4UniformRand()<0.5) {
307 selec1.Tprod = MarsN;
308 } else {
309 selec1.Tprod = MarsP;
310 }
311 selec1.V10 = 2.0;
312 }
313
314 //if ( SelBS(pType, fANucl, fZNucl) >0.0 ) AddSecondary();
315 if ( SelBS(pType, fANucl, fZNucl) >0.0 ) AddSecondaryToMarsList();
316
317}
318
319void G4Mars5GeV::CreatePion(G4int ib, G4int pType, G4double pE)
320{
321#ifdef G4VERBOSE
322 if (GetVerboseLevel() > 2) {
323 G4cout << " G4Mars5GeV::CreatePion()" << G4endl;
324 }
325#endif
326 static const G4double PionProductionEth = 0.28*GeV;
327 static const G4double KaonProductionEth = 2.0*GeV;
328
329 if ( pE<PionProductionEth ) {
330 if ((ib==MarsP)||(ib==MarsN)) return;
331 pE += GetParticleDefinition(MarsPIminus)->GetPDGMass();
332 }
333 selec1.Einc = pE;
334
335 if ( ib == MarsP ) {
336 selec1.Treac = MarsP;
337 } else if ( ib == MarsN ) {
338 selec1.Treac = MarsN;
339 } else if ( ib == MarsPIplus ) {
340 selec1.Treac = MarsPIplus;
341 } else if ( ib == MarsPIminus ) {
342 selec1.Treac = MarsPIminus;
343 } else if ( ib == MarsKplus ) {
344 selec1.Treac = MarsPIplus;
345 } else if ( ib == MarsKminus ) {
346 selec1.Treac = MarsPIminus;
347 } else {
348 selec1.Treac = MarsPIminus;
349 }
350 if (G4UniformRand()<0.5) {
351 selec1.Tprod = MarsPIplus;
352 } else {
353 selec1.Tprod = MarsPIminus;
354 }
355 selec1.V10 = 2.1;
356 if ( SelBS(pType, fANucl, fZNucl) >0.0 ){
357 // change secondary into Kaon
358 if ( pE > KaonProductionEth ) {
359 if ( Rkaon(ib,selec1.Tprod,pE) > G4UniformRand()) {
360 if (selec1.Tprod==MarsPIminus) {
361 selec1.Tprod=MarsKminus;
362 } else {
363 selec1.Tprod=MarsKplus;
364 }
365 }
366 }
367 //AddSecondary();
368 AddSecondaryToMarsList();
369 }
370}
371
372void G4Mars5GeV::CreatePionZero(G4int ib, G4int pType, G4double pE)
373{
374#ifdef G4VERBOSE
375 if (GetVerboseLevel() > 2) {
376 G4cout << " G4Mars5GeV::CreatePionZero()" << G4endl;
377 }
378#endif
379 static const G4double PionProductionEth = 0.28*GeV;
380
381 if ( pE<PionProductionEth ) {
382 if ((ib==MarsP)||(ib==MarsN)) return;
383 }
384
385 if ( ib == MarsP ) {
386 selec1.Treac = MarsP;
387 } else if ( ib == MarsN ) {
388 selec1.Treac = MarsN;
389 } else if ( ib == MarsPIplus ) {
390 selec1.Treac = MarsPIplus;
391 } else if ( ib == MarsPIminus ) {
392 selec1.Treac = MarsPIminus;
393 } else if ( ib == MarsKplus ) {
394 selec1.Treac = MarsPIplus;
395 } else if ( ib == MarsKminus ) {
396 selec1.Treac = MarsPIminus;
397 } else {
398 selec1.Treac = MarsPIminus;
399 }
400 selec1.Tprod = MarsKplus;
401 selec1.V10 = 1.0;
402 if ( SelBS(pType, fANucl, fZNucl) >0.0 ) {
403 selec1.Tprod = MarsPI0;
404 //AddSecondary();
405 AddSecondaryToMarsList();
406 }
407}
408
409//void G4Mars5GeV::AddSecondary()
410void G4Mars5GeV::AddSecondaryToMarsList()
411{
412#ifdef G4VERBOSE
413 if (GetVerboseLevel() > 2) {
414 G4cout << " G4Mars5GeV::AddSecondaryToMarsList()" << G4endl;
415 G4cout << " Particle :" << selec1.Tprod;
416 G4cout << ":" << GetParticleName(selec1.Tprod) <<G4endl;
417 G4cout << " Energy :" << selec1.EN <<G4endl;
418 G4cout << "Weight :" << selec1.V << G4endl;
419 }
420#endif
421 // determine direction cosine
422 G4double g = 1.0;
423 while (g>=1.0) {
424 G4double g1 = G4UniformRand();
425 G4double g2 = G4UniformRand();
426 G4double gg = 2.0*g1 - 1.0;
427 g = gg*gg + g2*g2;
428 selec2.Ch = (gg*gg - g2*g2)/g;
429 selec2.Sh = 2.0*gg*g2/g;
430 }
431 G4ThreeVector pin = incidentParticle->Get4Momentum().vect().unit();
432 G4ThreeVector pout;
433
434 Trans(&pin, &pout);
435 if (numberOfSecondaries>=FastVectorSize) {
436 G4String text = " G4Mars5GeV::AddSecondaryToMarsList() too many secondaries";
437 throw G4HadronicException(__FILE__, __LINE__, text);
438 }
439
440 // create seconday Dynamic Particle
441 G4DynamicParticle* secondary =
442 new G4DynamicParticle(GetParticleDefinition(selec1.Tprod),
443 pout.unit(),
444 selec1.EN);
445 // add secondary into list
446
447 secondaries.SetElement(numberOfSecondaries, secondary);
448 weightOfSecondaries[numberOfSecondaries] = selec1.V;
449 numberOfSecondaries +=1;
450}
451
452G4double G4Mars5GeV::SelBS(G4int pType, G4double aNucl, G4double zNucl)
453{
454 static const G4double Atau= 0.2;
455 static const G4double Btau= 0.5*GeV;
456
457 G4int nc = 0;
458 G4int ip = selec1.Treac; // reaction particle type
459 G4int jp = selec1.Tprod; // procduction particle type
460 G4int jj = pType; // incident particle type
461 G4double e0 = selec1.Einc;
462 G4double en;
463 G4double v2 = 0.0;
464
465#ifdef G4VERBOSE
466 if (GetVerboseLevel() > 2) {
467 G4cout << " G4Mars5GeV::SelBS" << G4endl;
468 G4cout << " pType = " << pType << " e0 = " << e0 << G4endl;
469 G4cout << " aNucl = " << aNucl << " zNucl = " << zNucl << G4endl;
470 G4cout << " Treac = " <<selec1.Treac;
471 G4cout << " Tprod = " <<selec1.Tprod << G4endl;
472 }
473#endif
474
475 while(1){
476
477 G4double g1 = G4UniformRand();
478 G4double g2 = G4UniformRand();
479
480 // calculate energy
481 G4double dw = 0.0;
482 if (ip==jp) {
483 G4double ea = e0 * 0.01;
484 if (ea < selec3.Eth) {
485 dw = selec3.Emax-selec3.Eth;
486 en = selec3.Eth + g1*dw;
487 } else {
488 G4double cb = std::log(ea/selec3.Eth);
489 G4double ca = cb + 99.0;
490 if (g1<cb/ca) {
491 en = selec3.Eth*std::exp(g1*ca);
492 dw = en*ca;
493 } else {
494 en = ea*(g1*ca + 1.0 - cb);
495 dw = ea*ca;
496 }
497 }
498 } else {
499 en = selec3.Eth*std::pow(selec3.Emax/selec3.Eth, g1);
500 dw = en*std::log(selec3.Emax/selec3.Eth);
501 }
502
503 selec1.EN = en;
504
505#ifdef G4VERBOSE
506 if (GetVerboseLevel() > 2) {
507 G4cout << "selec1.EN = " << en << G4endl;
508 }
509#endif
510
511 if (en<0.5*MeV) {
512 selec1.V = 0.0;
513 return selec1.V;
514 }
515
516 // calculate direction cosine
517 G4double tau = en/Atau/e0*(Btau+e0)/GeV;
518 G4double c5 = 1.0-std::exp(-pi*tau);
519 G4double c4 = 1.0-g2*c5;
520 G4double t1 = -std::log(c4)/tau;
521 G4double rcs = std::cos(t1);
522 G4double rss = std::sqrt(1.0-rcs*rcs);
523 G4double da = 2.0*pi*rss*c5/(tau*c4);
524 selec2.Cs = rcs;
525 selec2.Ss = rss;
526
527 // select particle type
528 G4int ib = ip;
529 if (ip == MarsP) {
530 ib = MarsN;
531 } else if (ip == MarsN) {
532 ib = MarsP;
533 }
534 G4int jb = jp;
535 if ( ( jj==MarsGAM ) && ((jp!=MarsP)||(jp!=MarsN)) ){
536 jb = MarsKplus;
537 } else if (jp == MarsP) {
538 jb = MarsN;
539 } else if (jp == MarsN) {
540 jb = MarsP;
541 }
542
543 // calculate V
544 nc +=1;
545 v2 = dw*D2N2(jj, e0, en, t1, ib, jb, aNucl, zNucl)*da*(selec1.V10);
546#ifdef G4VERBOSE
547 if (GetVerboseLevel() > 2) {
548 G4cout << " D2N2 = " << v2/(dw*da*(selec1.V10));
549 G4cout << " v2 = " << v2 << G4endl;
550 }
551#endif
552
553 if (v2>0.0) break;
554
555 if (nc >=3) {
556 selec1.V = 0.0;
557#ifdef G4VERBOSE
558 if (GetVerboseLevel() > 2) {
559 G4cout << "exceed retry limit !!" << G4endl;
560 }
561#endif
562 return selec1.V;
563 }
564 }
565 selec1.V = v2;
566 return v2;
567}
568
569
570G4double G4Mars5GeV::D2N2(G4int pType, G4double incidentE,
571 G4double prodE, G4double tin,
572 G4int reacType, G4int proType,
573 G4double ai, G4double z)
574{
575 // Hadron inclusive yield at E0 < 5 GeV
576 // All parametrizations are based on
577 // the energy unit of MeV
578 //
579 // Original code is written by Nikolai Mokhov (Fermilab)
580 //C Copyright Nikolai Mokhov (Fermilab)
581 //C
582 //C MARS13(98)
583 //C
584 //C HADRON INCLUSIVE YIELD AT E0 < 5 GEV
585 //C-----
586 //C CREATED: 1979 BY B.SYCHEV
587 //C MODIFIED: 1979-1998 BY NVM
588 //C LAST CHANGE: 16-JUL-1998 BY NVM
589
590 static const G4double o2pi = 1./twopi;
591 static const G4double ospi = 1./std::sqrt(pi);
592
593 static G4double abu = 0.0;
594 static G4double alga = 0.0;
595 static G4double a13 = 1.0;
596 static G4double a23 = 1.0;
597 static G4double a125 = 1.0;
598 static G4double am25 = 0.0;
599 static G4double sqa = 1.0;
600 static G4double sqa1 = 0.0;
601 static G4double bm = 2.0;
602 static G4double sl;
603 static G4double sa;
604
605 // input of this method
606 G4double e0 = incidentE/MeV; // SHOULD BE GIVEN BY MEV !!!
607 G4double e = prodE/MeV; // SHOULD BE GIVEN BY MEV !!!
608 G4double t = tin;
609 G4int i = reacType;
610 G4int j = proType;
611 G4int jj = pType;
612
613 // output of this method
614 G4double d2n = 0.0;
615 G4double dnde = 0.0; // this value is not used anywhere else
616
617 if(ai<1.0) return 0.0;
618
619 G4double a = ai;
620 if(abu!=a)
621 {
622 abu = a;
623 if(a<=2.0)
624 {
625 alga = 0.0;
626 a13 = 1.0;
627 a23 = 1.0;
628 a125 = 1.0;
629 am25 = 0.0;
630 sqa = 1.0;
631 sqa1 = 0.0;
632 bm = 2.0;
633 }
634 else
635 {
636 alga = std::log(a);
637 a13 = std::pow(a,1./3.);
638 a23 = a13*a13;
639 a125 = std::pow(a,-1.25);
640 am25 = std::pow(a-1.0,0.25);
641 sqa = std::sqrt(a);
642 sqa1 = std::sqrt(a-1.);
643 bm = 1.0 + sqa;
644 }
645 sl = 0.72/std::pow(1.+alga,0.4);
646 sa = 0.087*a23 + 4.15;
647 }
648
649 G4double bn;
650 if(a<=2.0)
651 {
652 if(i*j<9 && t>halfpi) return 0.0;
653 bn = 2.0;
654 if(i>=3) bn = 1.0;
655 }
656 else
657 { bn = bm*std::exp(-sa*std::pow(3.68/e0,sl)); }
658
659 G4double emm = e0;
660 G4double e1ge = 0.001*e0;
661 G4double e2ge = e1ge*e1ge;
662 G4double f21 = 0.04/(e2ge*e2ge);
663 G4double f31 = 0.38*std::pow(e1ge,-0.65);
664 G4double f22 = 0.25/e2ge;
665 G4double f32 = am25*0.7/(e1ge+1.);
666 G4double ei2 = 0.0;
667 G4double x1 = f21 + f31;
668 if(x1<60.) ei2 = 0.8*std::exp(-x1);
669 G4double ei1 = ei2;
670 if(a>2.0)
671 {
672 ei1 = 0.;
673 x1 = f22 + f32;
674 if(x1<60.) ei1 = std::exp(-x1);
675 }
676
677 G4double ew1 = ei2;
678 G4double dnl = 0.0;
679 G4double dnl1 = 0.0;
680 G4double eli = 0.0;
681 if(ew1>=1.e-19)
682 {
683 G4double dli = 35.0*ew1/(a+69.0);
684 eli = 0.5*dli*e0;
685 if(i==j)
686 { dnl1 = dli*(2./3.)/e0; }
687 else
688 { dnl1 = dli*(1.-e/e0)/e0; }
689 dnl = dli*(5./3.-e/e0)/e0;
690 }
691
692 G4double qel = 1.0 - ei1;
693 if(a>2.0)
694 {
695 G4double e02 = std::pow(e0/350.,1.5);
696 G4double ex2 = 1.0;
697 if(e02<60.) ex2 = 1.0-std::exp(-e02);
698 qel = 0.0;
699 if(t<halfpi) qel = 1.17*ex2*std::exp(-0.08*sqa1)*(1.-ew1);
700 }
701
702 G4int in = i; // save i
703 G4double sw2 = (0.5+10.*e2ge/(2.+e1ge))*(4.+e0/470.);
704 G4double sql = sw2/(2.+e0/940.)-1.0;
705 G4double eql = 0.0;
706 G4double dnq = 0.0;
707
708 if(jj!=MarsGAM || j!=5)
709 {
710 if(qel>1.e-25)
711 {
712 eql = qel*e0*(sql+1.)/(sql+2.);
713 G4double bp1x = -60./std::log(e/e0);
714 if(sql<=bp1x)
715 {
716 bp1x = std::pow(e/e0,sql);
717 dnq = qel/e0*(sql+1.)*bp1x;
718 }
719 }
720 }
721
722 G4double bp1 = e0;
723 if(e0<1.e9) bp1 = std::sqrt(e0*e0+1880.*e0);
724 G4double pul = 1.e-3*bp1;
725 bp1 = 3.*std::pow(pul,0.25) - 2.0;
726 if(bp1<1.) bp1 = 1.;
727 G4double bpi = 0.0;
728 if(ei1>0.) bpi = bp1*std::exp(0.075*sqa1)*ei1;
729 G4double ec = 0.0;
730 if(a>2.0) ec = 10.5 - 0.02*a;
731 G4double g = 0.1*alga + 0.2;
732 G4double eog = std::pow(e0,g);
733 G4double f1 = 1./3.*ec*a/(1.8*eog);
734 x1 = 1.0;
735 if(f1<60.) x1 = 1.0 - std::exp(-f1);
736 G4double fm = 1.8*eog*x1;
737 G4double ez = ec + fm;
738 if(fm>=e0) ez = ec + e0;
739 G4double d = 1.0;
740 if(i>=3) d = 0.0;
741 x1 = 1.0;
742 if(a<=44.) x1 = std::exp(-std::exp(4.-a));
743 G4double ez2 = 33.5*a125*x1*(ez-ec)*(1.-ez/(ec-e0));
744 G4double epw = e0 - ez - (bn-d)*ec - 140.*bpi - ez2;
745 G4double e2 = epw - eli - eql;
746
747 G4double ak1 = 3.0;
748 G4double ak20 = 5.e-4*(1.+a13)*e0;
749 x1 = 1.0;
750 if(ak20<60.) x1 = 1.0 - std::exp(-ak20);
751 G4double ga = std::pow(e1ge,0.06)*ak1*x1;
752 G4double egr = e0/(ga+1.);
753 G4double d2 = 250.*(1.+2.5*e0*std::exp(-0.02*a)/(e0+1.e3))/sqa;
754 G4double aea = e2/(e0*(1./(1.+ga)-d2*std::log(1.+egr/d2)/e0));
755 aea *= 1./(1.+d2*(ga+1.)/(3.*e0*(d2/e0+1.75)));
756 if(i<=2 && j>=3)
757 {
758 emm = e0 - 140.;
759 if(j!=5 && a!=1. && (i+j)!=5) emm = e0 - 280;
760 }
761 G4double dn = 0.0;
762 if(e<=emm) dn = aea*(e0/emm)*std::pow(1.-e/emm,ga)/(e+d2);
763 if(i>=3) bpi += 1;
764 dnde = dn + dnq + dnl;
765 // In original code, check nupr. But in this code, nupr is aliways set to 0
766 // if(nupr==1) return;
767
768 G4double pna = bpi/bp1;
769 G4double pns = bn+bpi;
770
771 // Angular distribution
772 G4double qe = 0.0;
773 if((jj!=MarsGAM || j!=5)
774 &&(t<halfpi)
775 &&(i<3||i==j||j>4)
776 &&(i>2||j<3)
777 &&(a>2.0||i!=2||j!=1)
778 &&(qel>=1.e-26))
779 {
780 G4double d1 = 25.*(1.+0.008*e0*t);
781 G4double dp;
782 if(i==j)
783 { dp = 0.8; }
784 else
785 { dp = 0.2; }
786 if(a<=2.&&i==1&&j<=2) dp = 0.5;
787 if(a<=2.&&i==2&&j==2) dp = 1.0;
788 G4double eq = e0*sqr(std::cos(t))/(1.+e0*sqr(std::sin(t))/1880.) - 25.0;
789 G4double exq = sqr((e-eq)/d1) + 0.5*sw2*t*t;
790 if(exq<60.) qe = qel*sw2*std::exp(-exq)*dp*ospi/d1;
791 }
792
793 G4int iold = i;
794 if(i==3) i = 2;
795 if(i==4) i = 1;
796 G4bool condA = a<2. && i==2;
797 G4bool condB = a<2.;
798 G4double az(0);
799 G4double pn;
800 G4double pr(0);
801 if(!condB)
802 {
803 if(i==2)
804 { az = (z+1.)/(a-z); }
805 else
806 { az = (a+1.-z)/z; }
807 x1 = 0.5*e1ge;
808 pn = az;
809 if(x1<60.) pn *= 1. + std::exp(-x1);
810 if(i==j)
811 { pr = bn*pn/(1.+pn); }
812 else
813 { pr = bn/(1.+pn); }
814 }
815 if(condA || !condB)
816 {
817 if(i==1) az = z/(a-z);
818 if(i==2) az = (a-z)/z;
819 G4double bp = 1.0;
820 G4double e0g = e1ge*e2ge;
821 if(e0g<60.) bp -= 0.5*std::exp(-e0g);
822 G4double ap = az * bp;
823 bp = 6.*(1.+ap);
824 if((i==1&&j==3)||(i==2&&j==4)) pr = pna*(bp1/3.-(2.+ap)/bp);
825 if((i==2&&j==3)||(i==1&&j==4)) pr = pna*(bp1/3.+(3.-ap)/bp);
826 if(j==5) pr = pna*(bp1/3.+(2.*ap-1.)/bp);
827 }
828 if(condB)
829 {
830 switch(i)
831 {
832 case 1:
833 switch(j)
834 {
835 case 1:
836 case 2:
837 pr = bn/2.; break;
838 case 3:
839 case 4:
840 pr = pna*(bp1/3.-1./6.); break;
841 case 5:
842 pr = pna*(bp1+1.)/3.; break;
843 }
844 break;
845 case 2:
846 switch(j)
847 {
848 case 1:
849 pr = 0.33*ew1/bn; break;
850 case 2:
851 pr = (1.-0.33*ew1)/bn; break;
852 }
853 break;
854 }
855 }
856
857 G4double ek3 = 0.01*std::sqrt(e1ge)*(1.+alga/4.);
858 G4double tay = 200.*e0/(e0+560.);
859 G4double w = e/tay;
860 G4double ek4 = 1.21*e0*w/(std::sqrt(1.+alga)*(e0+2000.));
861 if(j>=3) ek4 = 0.3*w*(e0-1000.)/(e0+1000.);
862 G4double wpic = w*pi;
863 G4double w2 = w*w;
864 G4double ex8 = 1.0;
865 if(wpic<60.) ex8 /= 1.0 + std::exp(-wpic);
866 G4double ek = (1.+w2)*(1.+5.2*ek4/(2.+w2))*ex8;
867 G4double wtw = 2.*(std::sqrt(1.+ek3*e*t*1.e-3)-1.)/(tay*ek3*1.e-3)+ek4*t*t;
868 G4double sm = 0.0;
869 if(dn>=1.e-26 && wtw<60.) sm = pr*dn*ek*std::exp(-wtw)/pns;
870
871 G4double dl = 0.0;
872 i = iold;
873 if((dnl1>=1.e-20)
874 &&(i<3||i==j||j>4)
875 &&(i>2||j<3))
876 {
877 tay = 200.*e0/(e0+2600.);
878 w = e/tay;
879 ek4 = 1.21*e0*w/(std::sqrt(1.+alga)*(e0+2000.));
880 i = in;
881 wpic = w*pi;
882 w2 = w*w;
883 ex8 = 1.0;
884 if(wpic<60.) ex8 /= 1.0 + std::exp(-wpic);
885 wtw = w*t + ek4*t*t;
886 if(wtw<60.)
887 {
888 G4double ft = (1.+w2)*(1.+5.2*ek4/(2.+w2))*std::exp(-wtw)*ex8;
889 if(ft>=1.-16)
890 {
891 G4double dp;
892 if(i==j)
893 { dp = 2./3.; }
894 else
895 { dp = 1./3.; }
896 if(a<=2.&&i==1&&j==2) dp = 0.5;
897 dl = dp*dnl1*ft;
898 }
899 }
900 }
901
902 if(jj==MarsGAM && j==5)
903 {
904 sm *= 0.6;
905 dl *= 2.0;
906 }
907 d2n = o2pi*(qe+sm+dl);
908 // d2n value is calculated in unit of [1/MeV]
909
910 d2n *= (1./MeV);
911 return d2n;
912}
913
914
915G4double G4Mars5GeV::Rkaon(G4int ib, G4int jp, G4double eRaw)
916{
917 // Energy dependent K/pi ratio
918 // Parametrizations are valid for energy range of
919 // incident particle as 2.0 GeV to 100 GeV
920 // All parametrizations in this method are based on
921 // the energy unit of GeV.
922 //
923 // Original code is written by Nikolai Mokhov (Fermilab)
924 //C Copyright Nikolai Mokhov (Fermilab)
925 //C
926 //C MARS13(98)
927 //C ENERGY DEPENDENT K/PI RATIO
928 //C FOR GIVEN TREEM AND SELMO PARAMETERS
929 //C-----
930 //C CREATED: 1996 BY N.MOKHOV (NVM)
931 //C LAST CHANGE: 12-FEB-1996 BY NVM
932
933 static const G4double rkp = 0.071;
934 static const G4double rkm = 0.083;
935 static const G4double al2 = 0.69314718;
936 static const G4double al100 = 4.6051702;
937 static const G4double al21 = 3.0445224;
938 static const G4double al51 = 3.9318256;
939
940 G4double eGeV = eRaw / GeV;
941 G4double rK = 0.;
942 if(eGeV < 2.1) return rK;
943 G4double ale = std::log(eGeV);
944
945 // No.1
946 rK = rkp;
947 if(jp == MarsPIminus) rK = rkm;
948 if(ib == MarsPIplus || ib == MarsPIminus) rK *= 1.3;
949 else if(ib == MarsKplus || ib == MarsKminus) rK *= 2.0;
950 G4double rK1 = rK;
951 if(eGeV<100.)
952 {
953 G4double rmi = 0.03;
954 if(ib >= MarsPIplus) rmi = 0.08;
955 rK1 = rmi + (rK-rmi)*(ale-al2)/(al100-al2);
956 }
957
958 // No.2
959 if(eGeV<=5.2 || eGeV>=51.0) {
960 rK = 1.3*rK1;
961 } else if(eGeV<7.2) {
962 rK = rK1*(1.3+0.15*(eGeV-5.2));
963 } else if(eGeV<21.) {
964 rK = 1.6*rK1;
965 } else {
966 rK = rK1*(1.3+0.3*(al51-ale)/(al51-al21));
967 }
968
969 return rK;
970}
971
972void G4Mars5GeV::Trans(G4ThreeVector* d1,G4ThreeVector* d2)
973{
974
975#ifdef G4VERBOSE
976 if (GetVerboseLevel() > 2) {
977 G4cout << " G4Mars5GeV::Trans() " << G4endl;
978 }
979#endif
980 // Direction cosine transformation
981 // using (cs,ss,ch,sh)
982
983 // inputs
984 G4double cs = selec2.Cs;
985 G4double ss = selec2.Ss;
986 G4double ch = selec2.Ch;
987 G4double sh = selec2.Sh;
988
989 G4double sss, ttt, uuu;
990 G4double dx1 = d1->x();
991 G4double dy1 = d1->y();
992 G4double dz1 = d1->z();
993 G4double sz = dx1*dx1 + dy1*dy1;
994 if(sz > 1.e-50)
995 {
996 sz = std::sqrt(sz);
997 sss = ss*(ch*dz1*dx1-sh*dy1)/sz + cs*dx1;
998 ttt = ss*(ch*dz1*dy1+sh*dx1)/sz + cs*dy1;
999 uuu = - ss*ch*sz + cs*dz1;
1000 }
1001 else
1002 {
1003 sss = ss*ch + dx1;
1004 ttt = ss*sh + dy1;
1005 uuu = cs*dz1;
1006 }
1007 G4double den = std::sqrt(sss*sss+uuu*uuu+ttt*ttt);
1008 d2->setX(sss/den);
1009 d2->setY(ttt/den);
1010 d2->setZ(uuu/den);
1011 return;
1012}
Note: See TracBrowser for help on using the repository browser.