source: trunk/source/processes/hadronic/management/src/G4HadronicProcessStore.cc @ 1034

Last change on this file since 1034 was 1007, checked in by garnier, 15 years ago

update to geant4.9.2
File size: 18.3 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// $Id: G4HadronicProcessStore.cc,v 1.7 2008/10/22 07:58:20 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-02 $
28//
29// -------------------------------------------------------------------
30//
31// GEANT4 Class file
32//
33//
34// File name:     G4HadronicProcessStore
35//
36// Author:        Vladimir Ivanchenko
37//
38// Creation date: 09.05.2008
39//
40// Modifications:
41//
42//
43// Class Description:
44//
45// -------------------------------------------------------------------
46//
47//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
48//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
49
50#include "G4HadronicProcessStore.hh"
51#include "G4Element.hh"
52#include "G4ProcessManager.hh"
53#include "G4Electron.hh"
54#include "G4Proton.hh"
55
56G4HadronicProcessStore* G4HadronicProcessStore::theInstance = 0;
57
58//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
59
60G4HadronicProcessStore* G4HadronicProcessStore::Instance()
61{
62  if(0 == theInstance) {
63    static G4HadronicProcessStore manager;
64    theInstance = &manager;
65  }
66  return theInstance;
67}
68
69//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
70
71G4HadronicProcessStore::~G4HadronicProcessStore()
72{
73  /*
74  for (G4int i=0; i<n_proc; i++) {
75    if( process[i] ) delete process[i];
76  }
77  */
78}
79
80//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
81
82G4HadronicProcessStore::G4HadronicProcessStore()
83{
84  n_proc = 0;
85  n_part = 0;
86  n_model= 0;
87  n_extra= 0;
88  currentProcess  = 0;
89  currentParticle = 0;
90  verbose = 1;
91  buildTableStart = true;
92}
93
94//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
95
96G4double G4HadronicProcessStore::GetElasticCrossSectionPerVolume(
97    const G4ParticleDefinition *aParticle,
98    G4double kineticEnergy,
99    const G4Material *material)
100{
101  G4double cross = 0.0;
102  const G4ElementVector* theElementVector = material->GetElementVector();
103  const G4double* theAtomNumDensityVector = material->GetVecNbOfAtomsPerVolume();
104  size_t nelm = material->GetNumberOfElements();
105  for (size_t i=0; i<nelm; i++) {
106    const G4Element* elm = (*theElementVector)[i];
107    cross += theAtomNumDensityVector[i]*
108      GetElasticCrossSectionPerAtom(aParticle,kineticEnergy,elm);
109  }
110  return cross;
111}
112
113//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
114
115G4double G4HadronicProcessStore::GetElasticCrossSectionPerAtom(
116    const G4ParticleDefinition *aParticle,
117    G4double kineticEnergy,
118    const G4Element *anElement)
119{
120  G4HadronicProcess* hp = FindProcess(aParticle, fHadronElastic);
121  localDP.SetKineticEnergy(kineticEnergy);
122  G4double cross = 0.0;
123  if(hp) cross = hp->GetMicroscopicCrossSection(&localDP,
124                                                anElement,
125                                                STP_Temperature);
126  return cross;
127}
128
129//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
130
131G4double G4HadronicProcessStore::GetElasticCrossSectionPerIsotope(
132    const G4ParticleDefinition*,
133    G4double,
134    G4int, G4int)
135{
136  return 0.0;
137}
138
139//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
140
141G4double G4HadronicProcessStore::GetInelasticCrossSectionPerVolume(
142    const G4ParticleDefinition *aParticle,
143    G4double kineticEnergy,
144    const G4Material *material)
145{
146  G4double cross = 0.0;
147  const G4ElementVector* theElementVector = material->GetElementVector();
148  const G4double* theAtomNumDensityVector = material->GetVecNbOfAtomsPerVolume();
149  size_t nelm = material->GetNumberOfElements();
150  for (size_t i=0; i<nelm; i++) {
151    const G4Element* elm = (*theElementVector)[i];
152    cross += theAtomNumDensityVector[i]*
153      GetInelasticCrossSectionPerAtom(aParticle,kineticEnergy,elm);
154  }
155  return cross;
156}
157
158//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
159
160G4double G4HadronicProcessStore::GetInelasticCrossSectionPerAtom(
161    const G4ParticleDefinition *aParticle,
162    G4double kineticEnergy,
163    const G4Element *anElement)
164{
165  G4HadronicProcess* hp = FindProcess(aParticle, fHadronInelastic);
166  localDP.SetDefinition(const_cast<G4ParticleDefinition*>(aParticle));
167  localDP.SetKineticEnergy(kineticEnergy);
168  G4double cross = 0.0;
169  if(hp) cross = hp->GetMicroscopicCrossSection(&localDP,
170                                                anElement,
171                                                STP_Temperature);
172  return cross;
173}
174
175//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
176
177G4double G4HadronicProcessStore::GetInelasticCrossSectionPerIsotope(
178    const G4ParticleDefinition *,
179    G4double,
180    G4int, G4int)
181{
182  return 0.0;
183}
184
185//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
186
187G4double G4HadronicProcessStore::GetCaptureCrossSectionPerVolume(
188    const G4ParticleDefinition *aParticle,
189    G4double kineticEnergy,
190    const G4Material *material)
191{
192  G4double cross = 0.0;
193  const G4ElementVector* theElementVector = material->GetElementVector();
194  const G4double* theAtomNumDensityVector = material->GetVecNbOfAtomsPerVolume();
195  size_t nelm = material->GetNumberOfElements();
196  for (size_t i=0; i<nelm; i++) {
197    const G4Element* elm = (*theElementVector)[i];
198    cross += theAtomNumDensityVector[i]*
199      GetCaptureCrossSectionPerAtom(aParticle,kineticEnergy,elm);
200  }
201  return cross;
202}
203
204//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
205
206G4double G4HadronicProcessStore::GetCaptureCrossSectionPerAtom(
207    const G4ParticleDefinition *aParticle,
208    G4double kineticEnergy,
209    const G4Element *anElement)
210{
211  G4HadronicProcess* hp = FindProcess(aParticle, fCapture);
212  localDP.SetDefinition(const_cast<G4ParticleDefinition*>(aParticle));
213  localDP.SetKineticEnergy(kineticEnergy);
214  G4double cross = 0.0;
215  if(hp) cross = hp->GetMicroscopicCrossSection(&localDP,
216                                                anElement,
217                                                STP_Temperature);
218  return cross;
219}
220
221//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
222
223G4double G4HadronicProcessStore::GetCaptureCrossSectionPerIsotope(
224    const G4ParticleDefinition *,
225    G4double,
226    G4int, G4int)
227{
228  return 0.0;
229}
230
231//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
232
233G4double G4HadronicProcessStore::GetFissionCrossSectionPerVolume(
234    const G4ParticleDefinition *aParticle,
235    G4double kineticEnergy,
236    const G4Material *material)
237{
238  G4double cross = 0.0;
239  const G4ElementVector* theElementVector = material->GetElementVector();
240  const G4double* theAtomNumDensityVector = material->GetVecNbOfAtomsPerVolume();
241  size_t nelm = material->GetNumberOfElements();
242  for (size_t i=0; i<nelm; i++) {
243    const G4Element* elm = (*theElementVector)[i];
244    cross += theAtomNumDensityVector[i]*
245      GetFissionCrossSectionPerAtom(aParticle,kineticEnergy,elm);
246  }
247  return cross;
248}
249
250//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
251
252G4double G4HadronicProcessStore::GetFissionCrossSectionPerAtom(
253    const G4ParticleDefinition *aParticle,
254    G4double kineticEnergy,
255    const G4Element *anElement)
256{
257  G4HadronicProcess* hp = FindProcess(aParticle, fFission);
258  localDP.SetDefinition(const_cast<G4ParticleDefinition*>(aParticle));
259  localDP.SetKineticEnergy(kineticEnergy);
260  G4double cross = 0.0;
261  if(hp) cross = hp->GetMicroscopicCrossSection(&localDP,
262                                                anElement,
263                                                STP_Temperature);
264  return cross;
265}
266
267//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
268
269G4double G4HadronicProcessStore::GetFissionCrossSectionPerIsotope(
270    const G4ParticleDefinition *,
271    G4double,
272    G4int, G4int)
273{
274  return 0.0;
275}
276
277//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
278
279G4double G4HadronicProcessStore::GetChargeExchangeCrossSectionPerVolume(
280    const G4ParticleDefinition *aParticle,
281    G4double kineticEnergy,
282    const G4Material *material)
283{
284  G4double cross = 0.0;
285  const G4ElementVector* theElementVector = material->GetElementVector();
286  const G4double* theAtomNumDensityVector = material->GetVecNbOfAtomsPerVolume();
287  size_t nelm = material->GetNumberOfElements();
288  for (size_t i=0; i<nelm; i++) {
289    const G4Element* elm = (*theElementVector)[i];
290    cross += theAtomNumDensityVector[i]*
291      GetChargeExchangeCrossSectionPerAtom(aParticle,kineticEnergy,elm);
292  }
293  return cross;
294}
295
296//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
297
298G4double G4HadronicProcessStore::GetChargeExchangeCrossSectionPerAtom(
299    const G4ParticleDefinition *aParticle,
300    G4double kineticEnergy,
301    const G4Element *anElement)
302{
303  G4HadronicProcess* hp = FindProcess(aParticle, fChargeExchange);
304  localDP.SetDefinition(const_cast<G4ParticleDefinition*>(aParticle));
305  localDP.SetKineticEnergy(kineticEnergy);
306  G4double cross = 0.0;
307  if(hp) cross = hp->GetMicroscopicCrossSection(&localDP,
308                                                anElement,
309                                                STP_Temperature);
310  return cross;
311}
312
313//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
314
315G4double G4HadronicProcessStore::GetChargeExchangeCrossSectionPerIsotope(
316    const G4ParticleDefinition *,
317    G4double,
318    G4int, G4int)
319{
320  return 0.0;
321}
322
323//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
324
325void G4HadronicProcessStore::Register(G4HadronicProcess* proc) 
326{ 
327  for(G4int i=0; i<n_proc; i++) {if(process[i] == proc) return;}
328   
329  n_proc++;
330  process.push_back(proc);
331}
332
333//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
334
335void G4HadronicProcessStore::RegisterParticle(G4HadronicProcess* proc, 
336                                              const G4ParticleDefinition* part) 
337{ 
338  G4int i=0;
339  for(; i<n_proc; i++) {if(process[i] == proc) break;}
340  G4int j=0;
341  for(; j<n_part; j++) {if(particle[j] == part) break;}
342
343  if(j == n_part) {
344    n_part++;
345    particle.push_back(part);
346    wasPrinted.push_back(0);
347  }
348 
349  // the pair should be added?
350  if(i < n_proc) {
351    std::multimap<PD,HP,std::less<PD> >::iterator it;
352    for(it=p_map.lower_bound(part); it!=p_map.upper_bound(part); ++it) {
353      if(it->first == part) {
354        HP process = (it->second);
355        if(proc == process) return;
356      }
357    }
358  }
359 
360  p_map.insert(std::multimap<PD,HP>::value_type(part,proc));
361}
362
363//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
364
365void G4HadronicProcessStore::RegisterInteraction(G4HadronicProcess* proc,
366                                                 G4HadronicInteraction* mod)
367{
368  G4int i=0;
369  for(; i<n_proc; i++) {if(process[i] == proc) break;}
370  G4int k=0;
371  for(; k<n_model; k++) {if(model[k] == mod) break;}
372   
373  m_map.insert(std::multimap<HP,HI>::value_type(proc,mod));
374   
375  if(k == n_model) {
376    n_model++;
377    model.push_back(mod);
378    modelName.push_back(mod->GetModelName());
379  }
380}
381
382//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
383
384void G4HadronicProcessStore::DeRegister(G4HadronicProcess* proc)
385{
386  for(G4int i=0; i<n_proc; i++) {
387    if(process[i] == proc) {
388      process[i] = 0;
389      break;
390    }
391  }
392} 
393
394//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
395
396void G4HadronicProcessStore::RegisterExtraProcess(G4VProcess* proc)
397{
398  for(G4int i=0; i<n_extra; i++) {if(extraProcess[i] == proc) return;}
399   
400  n_extra++;
401  extraProcess.push_back(proc);
402}
403
404//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
405
406void G4HadronicProcessStore::RegisterParticleForExtraProcess(
407                             G4VProcess* proc,
408                             const G4ParticleDefinition* part)
409{
410  G4int i=0;
411  for(; i<n_extra; i++) {if(extraProcess[i] == proc) break;}
412  G4int j=0;
413  for(; j<n_part; j++) {if(particle[j] == part) break;}
414
415  if(j == n_part) {
416    n_part++;
417    particle.push_back(part);
418    wasPrinted.push_back(0);
419  }
420 
421  // the pair should be added?
422  if(i < n_extra) {
423    std::multimap<PD,G4VProcess*,std::less<PD> >::iterator it;
424    for(it=ep_map.lower_bound(part); it!=ep_map.upper_bound(part); ++it) {
425      if(it->first == part) {
426        G4VProcess* process = (it->second);
427        if(proc == process) return;
428      }
429    }
430  }
431 
432  ep_map.insert(std::multimap<PD,G4VProcess*>::value_type(part,proc));
433} 
434
435//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
436
437void G4HadronicProcessStore::DeRegisterExtraProcess(G4VProcess* proc)
438{
439  for(G4int i=0; i<n_extra; i++) {
440    if(extraProcess[i] == proc) {
441      extraProcess[i] = 0;
442      break;
443    }
444  }
445}
446
447//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
448
449void G4HadronicProcessStore::PrintInfo(const G4ParticleDefinition* part) 
450{
451  if(buildTableStart && part == particle[n_part - 1]) {
452    buildTableStart = false;
453    Dump(verbose);
454  }
455}
456
457//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
458
459void G4HadronicProcessStore::Dump(G4int level)
460{
461  if(level > 0) {
462    G4cout << "=============================================================="
463           << "=============================="
464           << G4endl;
465      G4cout << "             HADRONIC PROCESSES SUMMARY (verbose level " << level
466             << ")" << G4endl;
467  }
468  for(G4int i=0; i<n_part; i++) {
469    PD part = particle[i];
470    G4String pname = part->GetParticleName();
471    G4bool yes = false;
472    if(level >= 2) yes = true;
473    else if(level == 1 && (pname == "proton" || 
474                           pname == "neutron" ||
475                           pname == "pi+" ||
476                           pname == "pi-" ||
477                           pname == "gamma" ||
478                           pname == "e-" ||
479                           pname == "mu-" ||
480                           pname == "kaon+" ||
481                           pname == "kaon-" ||
482                           pname == "lambda" ||
483                           pname == "anti_neutron" ||
484                           pname == "anti_proton")) yes = true;
485    if(yes) {
486      // main processes
487      std::multimap<PD,HP,std::less<PD> >::iterator it;
488      for(it=p_map.lower_bound(part); it!=p_map.upper_bound(part); ++it) {
489        if(it->first == part) {
490          HP proc = (it->second);
491          G4int j=0;
492          for(; j<n_proc; j++) {
493            if(process[j] == proc) {
494              Print(j, i);
495            }
496          }
497        }
498      }
499      // extra processes
500      std::multimap<PD,G4VProcess*,std::less<PD> >::iterator itp;
501      for(itp=ep_map.lower_bound(part); itp!=ep_map.upper_bound(part); ++itp) {
502        if(itp->first == part) {
503          G4VProcess* proc = (itp->second);
504          if(wasPrinted[i] == 0) {
505            wasPrinted[i] = 1;
506            G4cout<<G4endl;
507            G4cout << "                     Hadronic Processes for <" 
508                   <<part->GetParticleName() << ">" << G4endl; 
509          }
510          G4cout << "          " << proc->GetProcessName() << G4endl;
511        }
512      }
513    }
514  }
515  if(level > 0) {
516    G4cout << "=============================================================="
517           << "=============================="
518           << G4endl;
519  }
520}
521
522//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
523
524void G4HadronicProcessStore::Print(G4int idxProc, G4int idxPart)
525{
526  G4HadronicProcess* proc = process[idxProc];
527  const G4ParticleDefinition* part = particle[idxPart];
528  if(wasPrinted[idxPart] == 0) {
529    wasPrinted[idxPart] = 1;
530    G4cout<<G4endl;
531    G4cout << "                     Hadronic Processes for <" 
532           <<part->GetParticleName() << ">" << G4endl; 
533  }
534  HI hi = 0;
535  G4bool first;
536  std::multimap<HP,HI,std::less<HP> >::iterator ih;
537  G4cout << std::setw(20) << proc->GetProcessName() 
538         << "  Models: ";
539  first = true;
540  for(ih=m_map.lower_bound(proc); ih!=m_map.upper_bound(proc); ++ih) {
541    if(ih->first == proc) {
542      hi = ih->second;
543      G4int i=0;
544      for(; i<n_model; i++) {
545        if(model[i] == hi) break;
546      }
547      if(!first) G4cout << "                              ";
548      first = false;
549      G4cout << std::setw(25) << modelName[i] 
550             << ": Emin(GeV)= " 
551             << std::setw(5) << hi->GetMinEnergy()/GeV
552             << "  Emax(GeV)= " 
553             << hi->GetMaxEnergy()/GeV
554             << G4endl;
555    }
556  }
557}
558
559//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
560
561void G4HadronicProcessStore::SetVerbose(G4int val)
562{
563  verbose = val;
564  G4int i;
565  for(i=0; i<n_proc; i++) {
566    if(process[i]) process[i]->SetVerboseLevel(val);
567  }
568  for(i=0; i<n_model; i++) {
569    if(model[i]) model[i]->SetVerboseLevel(val);
570  }
571}
572
573//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
574
575G4int G4HadronicProcessStore::GetVerbose()
576{
577  return verbose;
578}
579
580//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
581
582G4HadronicProcess* G4HadronicProcessStore::FindProcess(
583   const G4ParticleDefinition* part, G4HadronicProcessType subType)
584{
585  bool isNew = false;
586  G4HadronicProcess* hp = 0;
587
588  if(part != currentParticle) {
589    isNew = true;
590    currentParticle = part;
591    localDP.SetDefinition(const_cast<G4ParticleDefinition*>(part));
592  } else if(!currentProcess) {
593    isNew = true;
594  } else if(subType == currentProcess->GetProcessSubType()) {
595    hp = currentProcess;
596  } else {
597    isNew = true;
598  }
599
600  if(isNew) {
601    std::multimap<PD,HP,std::less<PD> >::iterator it;
602    for(it=p_map.lower_bound(part); it!=p_map.upper_bound(part); ++it) {
603      if(it->first == part && subType == (it->second)->GetProcessSubType()) {
604        hp = it->second;
605        break;
606      }
607    } 
608    currentProcess = hp;
609  }
610
611  return hp;
612}
613
614//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
Note: See TracBrowser for help on using the repository browser.