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

Last change on this file since 819 was 819, checked in by garnier, 16 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.