source: trunk/examples/extended/field/field04/src/F04PhysicsList.cc @ 1027

Last change on this file since 1027 was 807, checked in by garnier, 16 years ago

update
File size: 21.0 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//
28
29#include "F04PhysicsList.hh"
30#include "F04PhysicsListMessenger.hh"
31
32#include "F04ExtraPhysics.hh"
33//#include "F04OpticalPhysics.hh"
34
35#include "G4DecayPhysics.hh"
36
37#include "G4EmStandardPhysics_option1.hh"
38#include "G4EmStandardPhysics_option2.hh"
39#include "G4EmStandardPhysics.hh"
40
41#include "G4EmExtraPhysics.hh"
42#include "G4EmProcessOptions.hh"
43
44#include "G4HadronElasticPhysics.hh"
45#include "G4HadronQElasticPhysics.hh"
46#include "G4HadronHElasticPhysics.hh"
47#include "G4NeutronTrackingCut.hh"
48#include "G4QStoppingPhysics.hh"
49#include "G4LHEPStoppingPhysics.hh"
50#include "G4IonBinaryCascadePhysics.hh"
51#include "G4IonPhysics.hh"
52
53#include "HadronPhysicsFTFP.hh"
54#include "HadronPhysicsFTFC.hh"
55
56#include "HadronPhysicsLHEP.hh"
57#include "HadronPhysicsLHEP_BERT.hh"
58#include "HadronPhysicsLHEP_EMV.hh"
59#include "HadronPhysicsLHEP_PRECO_HP.hh"
60
61#include "HadronPhysicsQGSC.hh"
62#include "HadronPhysicsQGSC_EFLOW.hh"
63
64#include "HadronPhysicsQGSP.hh"
65#include "HadronPhysicsQGSP_BERT.hh"
66#include "HadronPhysicsQGSP_BERT_HP.hh"
67#include "HadronPhysicsQGSP_BERT_TRV.hh"
68#include "HadronPhysicsQGSP_BIC.hh"
69#include "HadronPhysicsQGSP_BIC_HP.hh"
70
71#include "G4HadronInelasticQBBC.hh"
72#include "G4HadronInelasticQLHEP.hh"
73
74#include "G4IonPhysics.hh"
75
76#include "G4HadronProcessStore.hh"
77
78#include "G4LossTableManager.hh"
79
80#include "G4ProcessManager.hh"
81#include "G4ParticleTypes.hh"
82#include "G4ParticleTable.hh"
83
84#include "G4Gamma.hh"
85#include "G4Electron.hh"
86#include "G4Positron.hh"
87
88#include "F04StepMax.hh"
89
90#include "G4ProcessTable.hh"
91
92#include "G4PionDecayMakeSpin.hh"
93#include "G4DecayWithSpin.hh"
94
95#include "G4DecayTable.hh"
96#include "G4MuonDecayChannelWithSpin.hh"
97#include "G4MuonRadiativeDecayChannelWithSpin.hh"
98
99F04PhysicsList::F04PhysicsList(G4String physicsList) : G4VModularPhysicsList() 
100{
101    G4LossTableManager::Instance();
102
103    defaultCutValue  = 1.*mm;
104    fCutForGamma     = defaultCutValue;
105    fCutForElectron  = defaultCutValue;
106    fCutForPositron  = defaultCutValue;
107    fDump            = false;
108
109    fMessenger = new F04PhysicsListMessenger(this);
110
111    fEMPhysics     = new PhysicsListVector();
112    fHadronPhysics = new PhysicsListVector();
113
114    // Particles
115    fParticleList = new G4DecayPhysics("decays");
116
117    // EM physics
118    AddPhysicsList("emstandard_opt1");
119
120    // Set the default hadronic physics.
121    AddPhysicsList(physicsList);
122
123    stepMaxProcess = new F04StepMax();
124}
125
126F04PhysicsList::~F04PhysicsList()
127{
128    delete fMessenger;
129    delete fParticleList;
130
131    ClearEMPhysics();
132    ClearHadronPhysics();
133
134    delete fEMPhysics;
135    delete fHadronPhysics;
136
137    delete stepMaxProcess;
138}
139
140void F04PhysicsList::ClearEMPhysics()
141{
142    for (PhysicsListVector::iterator p  = fEMPhysics->begin();
143                                     p != fEMPhysics->end(); ++p) {
144        delete (*p);
145    }
146    fEMPhysics->clear();
147}
148
149void F04PhysicsList::ClearHadronPhysics()
150{
151    for (PhysicsListVector::iterator p  = fHadronPhysics->begin();
152                                     p != fHadronPhysics->end(); ++p) {
153        delete (*p);
154    }
155    fHadronPhysics->clear();
156}
157
158void F04PhysicsList::ConstructParticle()
159{
160    fParticleList->ConstructParticle();
161
162    G4DecayTable* MuonPlusDecayTable = new G4DecayTable();
163    MuonPlusDecayTable -> Insert(new
164                           G4MuonDecayChannelWithSpin("mu+",0.986));
165    MuonPlusDecayTable -> Insert(new
166                           G4MuonRadiativeDecayChannelWithSpin("mu+",0.014));
167    G4MuonPlus::MuonPlusDefinition() -> SetDecayTable(MuonPlusDecayTable);
168
169    G4DecayTable* MuonMinusDecayTable = new G4DecayTable();
170    MuonMinusDecayTable -> Insert(new
171                            G4MuonDecayChannelWithSpin("mu-",0.986));
172    MuonMinusDecayTable -> Insert(new
173                            G4MuonRadiativeDecayChannelWithSpin("mu-",0.014));
174    G4MuonMinus::MuonMinusDefinition() -> SetDecayTable(MuonMinusDecayTable);
175
176}
177
178void F04PhysicsList::ConstructProcess()
179{
180    AddTransportation();
181
182    for (PhysicsListVector::iterator p  = fEMPhysics->begin();
183                                     p != fEMPhysics->end(); ++p) {
184        (*p)->ConstructProcess();
185    }
186
187    fParticleList->ConstructProcess();
188
189    G4DecayWithSpin* decayWithSpin = new G4DecayWithSpin();
190
191    G4ProcessTable* processTable = G4ProcessTable::GetProcessTable();
192
193    G4VProcess* decay;
194    decay = processTable->FindProcess("Decay",G4MuonPlus::MuonPlus());
195
196    G4ProcessManager* fManager;
197    fManager = G4MuonPlus::MuonPlus()->GetProcessManager();
198
199    if (fManager) {
200      if (decay) fManager->RemoveProcess(decay);
201      fManager->AddProcess(decayWithSpin);
202      // set ordering for PostStepDoIt and AtRestDoIt
203      fManager ->SetProcessOrdering(decayWithSpin, idxPostStep);
204      fManager ->SetProcessOrdering(decayWithSpin, idxAtRest);
205    }
206
207    decay = processTable->FindProcess("Decay",G4MuonMinus::MuonMinus());
208
209    fManager = G4MuonMinus::MuonMinus()->GetProcessManager();
210
211    if (fManager) {
212      if (decay) fManager->RemoveProcess(decay);
213      fManager->AddProcess(decayWithSpin);
214      // set ordering for PostStepDoIt and AtRestDoIt
215      fManager ->SetProcessOrdering(decayWithSpin, idxPostStep);
216      fManager ->SetProcessOrdering(decayWithSpin, idxAtRest);
217    }
218
219    G4PionDecayMakeSpin* poldecay = new G4PionDecayMakeSpin();
220
221    decay = processTable->FindProcess("Decay",G4PionPlus::PionPlus());
222
223    fManager = G4PionPlus::PionPlus()->GetProcessManager();
224
225    if (fManager) {
226      if (decay) fManager->RemoveProcess(decay);
227      fManager->AddProcess(poldecay);
228      // set ordering for PostStepDoIt and AtRestDoIt
229      fManager ->SetProcessOrdering(poldecay, idxPostStep);
230      fManager ->SetProcessOrdering(poldecay, idxAtRest);
231    }
232
233    decay = processTable->FindProcess("Decay",G4PionMinus::PionMinus());
234
235    fManager = G4PionMinus::PionMinus()->GetProcessManager();
236
237    if (fManager) {
238      if (decay) fManager->RemoveProcess(decay);
239      fManager->AddProcess(poldecay);
240      // set ordering for PostStepDoIt and AtRestDoIt
241      fManager ->SetProcessOrdering(poldecay, idxPostStep);
242      fManager ->SetProcessOrdering(poldecay, idxAtRest);
243    }
244
245    for (PhysicsListVector::iterator p  = fHadronPhysics->begin();
246                                     p != fHadronPhysics->end();++p) {
247        (*p)->ConstructProcess();
248    }
249
250    AddStepMax();
251
252    if (fDump) G4HadronProcessStore::Instance()->Dump(1);
253}
254
255void F04PhysicsList::RemoveFromEMPhysicsList(const G4String& name)
256{
257    G4bool success = false;
258    for (PhysicsListVector::iterator p  = fEMPhysics->begin();
259                                     p != fEMPhysics->end(); ++p) {
260        G4VPhysicsConstructor* e = (*p);
261        if (e->GetPhysicsName() == name) {
262           fEMPhysics->erase(p);
263           success = true;
264           break;
265        }
266    }
267    if (!success) {
268       std::ostringstream message;
269       message << "PhysicsList::RemoveFromEMPhysicsList "<< name << "not found";
270       G4Exception(message.str().c_str());
271    }
272}
273
274void F04PhysicsList::RemoveFromHadronPhysicsList(const G4String& name)
275{
276    G4bool success = false;
277    for (PhysicsListVector::iterator p  = fHadronPhysics->begin();
278                                     p != fHadronPhysics->end(); ++p) {
279        G4VPhysicsConstructor* e = (*p);
280        if (e->GetPhysicsName() == name) {
281           fHadronPhysics->erase(p);
282           success = true;
283           break;
284        }
285    }
286    if (!success) {
287       G4String message("F04PhysicsList::RemoveFromHadronPhysicsList \"");
288       message += name;
289       message += " not found \"";
290       G4Exception(message);
291    }
292}
293
294void F04PhysicsList::AddPhysicsList(const G4String& name)
295{
296    if (verboseLevel>0) {
297        G4cout << "F04PhysicsList::AddPhysicsList: <" << name
298               << ">" << G4endl;
299    }
300
301    if (name == "emstandard_opt1") {
302
303       ClearEMPhysics();
304       fEMPhysics->push_back(new G4EmStandardPhysics_option1());
305       fEMPhysics->push_back(new F04ExtraPhysics());
306//       fEMPhysics->push_back(new F04OpticalPhysics());
307
308    } else if (name == "emstandard_opt2") {
309
310       ClearEMPhysics();
311       fEMPhysics->push_back(new G4EmStandardPhysics_option2());
312       fEMPhysics->push_back(new F04ExtraPhysics());
313//       fEMPhysics->push_back(new F04OpticalPhysics());
314
315   } else if (name == "FTFC") {
316
317       SetStandardList(false, false);
318       fHadronPhysics->push_back(new HadronPhysicsFTFC("hadron",true));
319
320    } else if (name == "FTFP") {
321
322       SetStandardList(false, false);
323       fHadronPhysics->push_back(new HadronPhysicsFTFP("hadron",true));
324
325    } else if (name == "FTFP_EMV") {
326
327       AddPhysicsList("emstandard_opt1");
328       AddPhysicsList("FTFP");
329
330    } else if (name == "LHEP") {
331
332       ClearHadronPhysics();
333       fHadronPhysics->push_back(new G4EmExtraPhysics("gamma_nuc"));
334       fHadronPhysics->push_back(new G4HadronElasticPhysics("LElastic",
335                                                           verboseLevel,
336                                                           false));
337       fHadronPhysics->push_back(new HadronPhysicsLHEP("hadron"));
338       fHadronPhysics->push_back(new G4IonPhysics("ion"));
339       fDump = true;
340 
341    } else if (name == "LHEP_BERT") {
342
343       ClearHadronPhysics();
344       fHadronPhysics->push_back(new G4EmExtraPhysics("gamma_nuc"));
345       fHadronPhysics->push_back(new G4HadronElasticPhysics("LElastic",
346                                                           verboseLevel,
347                                                           false));
348       fHadronPhysics->push_back(new HadronPhysicsLHEP_BERT("hadron"));
349       fHadronPhysics->push_back(new G4QStoppingPhysics("stopping",
350                                                       verboseLevel));
351       fHadronPhysics->push_back(new G4IonPhysics("ion"));
352       fDump = true;
353
354    } else if (name == "LHEP_EMV") {
355
356       ClearHadronPhysics();
357       AddPhysicsList("emstandard_opt1");
358       fHadronPhysics->push_back(new G4EmExtraPhysics("gamma_nuc"));
359       fHadronPhysics->push_back(new G4HadronElasticPhysics("LElastic",
360                                                           verboseLevel,
361                                                           false));
362       fHadronPhysics->push_back(new HadronPhysicsLHEP_EMV("hadron"));
363       fHadronPhysics->push_back(new G4QStoppingPhysics("stopping",
364                                                        verboseLevel));
365       fHadronPhysics->push_back(new G4IonPhysics("ion"));
366       fDump = true;
367
368    } else if (name == "LHEP_PRECO_HP") {
369
370       ClearHadronPhysics();
371       fHadronPhysics->push_back(new G4EmExtraPhysics("gamma_nuc"));
372       fHadronPhysics->push_back(new G4HadronElasticPhysics("LElastic",
373                                                           verboseLevel,
374                                                           false));
375       fHadronPhysics->push_back(new HadronPhysicsLHEP_PRECO_HP("hadron"));
376       fHadronPhysics->push_back(new G4QStoppingPhysics("stopping",
377                                                       verboseLevel));
378       fHadronPhysics->push_back(new G4IonPhysics("ion"));
379       fDump = true;
380
381    } else if (name == "QBBC") {
382
383       SetStandardList(false,true);
384       fHadronPhysics->push_back(new G4HadronInelasticQBBC("QBBC",
385                                                          verboseLevel,
386                                                          false,false,
387                                                          false,false,true));
388
389    } else if (name == "QBBCG") {
390
391       SetStandardList(false, false);
392       fHadronPhysics->push_back(new G4HadronInelasticQBBC("QBBC",
393                                                          verboseLevel,
394                                                          false,false,
395                                                          false,false,false));
396
397    } else if (name == "QBEC") {
398
399       SetStandardList(false,true);
400       fHadronPhysics->push_back(new G4HadronInelasticQBBC("QBBC",
401                                                          verboseLevel,
402                                                          false,true,
403                                                          false,false,true));
404
405    } else if (name == "QBBC_HP") {
406
407       SetStandardList(true,true);
408       fHadronPhysics->push_back(new G4HadronInelasticQBBC("QBBC",
409                                                          verboseLevel,
410                                                          false,false,
411                                                          false,true,true));
412
413    } else if (name == "QBEC_HP") {
414
415       SetStandardList(true,true);
416       fHadronPhysics->push_back(new G4HadronInelasticQBBC("QBBC",
417                                                          verboseLevel,
418                                                          false,true,
419                                                          false,true,true));
420
421    } else if (name == "QGSC") {
422
423       ClearHadronPhysics();
424       fHadronPhysics->push_back(new G4EmExtraPhysics("gamma_nuc"));
425       fHadronPhysics->push_back(new G4HadronQElasticPhysics("QElastic",
426                                                            verboseLevel));
427       fHadronPhysics->push_back(new HadronPhysicsQGSC("hadron",true));
428       fHadronPhysics->push_back(new G4QStoppingPhysics("stopping",
429                                                       verboseLevel,false));
430       fHadronPhysics->push_back(new G4IonPhysics("ion"));
431       fHadronPhysics->push_back(new G4NeutronTrackingCut());
432       fDump = true;
433
434    } else if (name == "QGSC_EFLOW") {
435
436       ClearHadronPhysics();
437       fHadronPhysics->push_back(new G4EmExtraPhysics("gamma_nuc"));
438       fHadronPhysics->push_back(new G4HadronQElasticPhysics("QElastic",
439                                                            verboseLevel));
440       fHadronPhysics->push_back(new HadronPhysicsQGSC_EFLOW("hadron",true));
441       fHadronPhysics->push_back(new G4QStoppingPhysics("stopping",
442                                                       verboseLevel,false));
443       fHadronPhysics->push_back(new G4IonPhysics("ion"));
444       fHadronPhysics->push_back(new G4NeutronTrackingCut());
445       fDump = true;
446
447    } else if (name == "QGSC_EMV") {
448
449       AddPhysicsList("emstandard_opt1");
450       AddPhysicsList("QGSC");
451
452    } else if (name == "QGSP") {
453
454       SetStandardList(false, false);
455       fHadronPhysics->push_back(new HadronPhysicsQGSP("hadron",true));
456
457    } else if (name == "QGSP_BERT") {
458
459       SetStandardList(false, false);
460       fHadronPhysics->push_back(new HadronPhysicsQGSP_BERT("hadron",true));
461
462    } else if (name == "QGSP_BERT_EMV") {
463
464       AddPhysicsList("emstandard_opt1");
465       AddPhysicsList("QGSP_BERT");
466
467    } else if (name == "QGSP_BERT_HP") {
468
469       SetStandardList(true, false);
470       fHadronPhysics->push_back(new HadronPhysicsQGSP_BERT_HP("hadron",true));
471
472    } else if (name == "QGSP_BERT_NQE") {
473
474       SetStandardList(false, false);
475       fHadronPhysics->push_back(new HadronPhysicsQGSP_BERT("hadron",false));
476
477    } else if (name == "QGSP_BERT_TRV") {
478
479       SetStandardList(false, false);
480       fHadronPhysics->push_back(new HadronPhysicsQGSP_BERT_TRV("hadron",true));
481
482    } else if (name == "QGSP_BIC") {
483
484       SetStandardList(false, false);
485       fHadronPhysics->push_back(new HadronPhysicsQGSP_BIC("hadron",true));
486
487    } else if (name == "QGSP_BIC_HP") {
488
489       SetStandardList(true, false);
490       fHadronPhysics->push_back(new HadronPhysicsQGSP_BIC_HP("hadron",true));
491
492    } else if (name == "QGSP_EMV") {
493
494       AddPhysicsList("emstandard_opt1");
495       AddPhysicsList("QGSP");
496
497    } else if (name == "QGSP_EMV_NQE") {
498
499       AddPhysicsList("emstandard_opt1");
500       AddPhysicsList("QGSP_NQE");
501
502    } else if (name == "QGSP_EMX") {
503
504       AddPhysicsList("emstandard_opt2");
505       AddPhysicsList("QGSP");
506
507    } else if (name == "QGSP_NQE") {
508
509       SetStandardList(false, false);
510       fHadronPhysics->push_back(new HadronPhysicsQGSP("hadron",false));
511
512    } else if (name == "QGSP_QEL") {
513
514       ClearHadronPhysics();
515       fHadronPhysics->push_back(new G4EmExtraPhysics("gamma_nuc"));
516       fHadronPhysics->push_back(new G4HadronQElasticPhysics("QElastic",
517                                                            verboseLevel));
518       fHadronPhysics->push_back(new HadronPhysicsQGSP("hadron",true));
519       fHadronPhysics->push_back(new G4QStoppingPhysics("stopping",
520                                                       verboseLevel));
521       fHadronPhysics->push_back(new G4IonPhysics("ion"));
522       fHadronPhysics->push_back(new G4NeutronTrackingCut());
523       fDump = true;
524
525    } else {
526
527       G4cout << "F04PhysicsList::AddPhysicsList: <" << name << ">"
528              << " is not defined" << G4endl;
529    }
530}
531
532void F04PhysicsList::SetStandardList(G4bool flagHP, G4bool glauber)
533{
534    ClearHadronPhysics();
535
536    fHadronPhysics->push_back(new G4EmExtraPhysics("gamma_nuc"));
537    fHadronPhysics->push_back(new G4HadronElasticPhysics("elastic",
538                                                        verboseLevel,
539                                                        flagHP,glauber));
540    fHadronPhysics->push_back(new  G4QStoppingPhysics("stopping",
541                                                     verboseLevel));
542    fHadronPhysics->push_back(new G4IonBinaryCascadePhysics("binary_ion"));
543    fHadronPhysics->push_back(new G4NeutronTrackingCut());
544
545    fDump = true;
546}
547
548void F04PhysicsList::SetCuts()
549{
550    if (verboseLevel >0) {
551        G4cout << "F04PhysicsList::SetCuts:";
552        G4cout << "CutLength : " << G4BestUnit(defaultCutValue,"Length")
553               << G4endl;
554    }
555
556    // set cut values for gamma at first and for e- second and next for e+,
557    // because some processes for e+/e- need cut values for gamma
558    SetCutValue(fCutForGamma, "gamma");
559    SetCutValue(fCutForElectron, "e-");
560    SetCutValue(fCutForPositron, "e+");
561
562    if (verboseLevel>0) DumpCutValuesTable();
563}
564
565void F04PhysicsList::SetCutForGamma(G4double cut)
566{
567    fCutForGamma = cut;
568    SetParticleCuts(fCutForGamma, G4Gamma::Gamma());
569}
570
571void F04PhysicsList::SetCutForElectron(G4double cut)
572{
573    fCutForElectron = cut;
574    SetParticleCuts(fCutForElectron, G4Electron::Electron());
575}
576
577void F04PhysicsList::SetCutForPositron(G4double cut)
578{
579    fCutForPositron = cut;
580    SetParticleCuts(fCutForPositron, G4Positron::Positron());
581}
582
583void F04PhysicsList::List()
584{
585    G4cout << "### PhysicsLists available: FTFC FTFP FTFP_EMV"
586           << " LHEP LHEP_BERT LHEP_EMV"
587           << G4endl;
588    G4cout << "                            LHEP_PRECO_HP QBBC QBBCG"
589           << " QBEC QBBC_HP QBEC_HP"
590           << G4endl;
591    G4cout << "                            QGSC QGSC_EFLOW QGSC_EMV"
592           << " QGSP QGSP_BERT QGSP_BER_EMV"
593           << G4endl;
594    G4cout << "                            QGSP_BERT_HP QGSP_BERT_NQE"
595           << " QGSP_BERT_TRV"
596           << G4endl;
597    G4cout << "                            QGSP_BIC QGSP_BIC_HP QGSP_EMV"
598           << " QGSP_EMV_NQE"
599           << G4endl;
600    G4cout << "                            QGSC_EMX QGSP_NQE QGSP_QEL"
601           << G4endl;
602}
603
604void F04PhysicsList::SetStepMax(G4double step)
605{
606  MaxChargedStep = step ;
607  stepMaxProcess->SetStepMax(MaxChargedStep);
608}
609
610F04StepMax* F04PhysicsList::GetStepMaxProcess()
611{
612  return stepMaxProcess;
613}
614
615void F04PhysicsList::AddStepMax()
616{
617  // Step limitation seen as a process
618
619  theParticleIterator->reset();
620  while ((*theParticleIterator)()){
621      G4ParticleDefinition* particle = theParticleIterator->value();
622      G4ProcessManager* pmanager = particle->GetProcessManager();
623
624      if (stepMaxProcess->IsApplicable(*particle) && !particle->IsShortLived())
625      {
626         if (pmanager) pmanager ->AddDiscreteProcess(stepMaxProcess);
627      }
628  }
629}
Note: See TracBrowser for help on using the repository browser.