From ad528b7143b1d1f6716db7eabd17a5975b0a86ac Mon Sep 17 00:00:00 2001
From: adrien matta <matta@lpccaen.in2p3.fr>
Date: Tue, 1 Oct 2019 17:30:52 +0200
Subject: [PATCH] * Adding custom level decay option - User can add a
level decay file in geant4 format - this is done in the event
generator file - the level scheme file should be saved in the tree ->
to be done
---
NPSimulation/Process/PhysicsList.cc | 778 ++++++++++++++------------
NPSimulation/Process/PhysicsList.hh | 2 +-
Projects/e748/Analysis.cxx | 388 ++++++-------
Projects/e748/Analysis.h | 8 +-
Projects/e748/configs/ConfigMust2.dat | 0
Projects/e748/e748.detector | 2 +-
6 files changed, 607 insertions(+), 571 deletions(-)
mode change 100644 => 100755 Projects/e748/configs/ConfigMust2.dat
mode change 100644 => 100755 Projects/e748/e748.detector
diff --git a/NPSimulation/Process/PhysicsList.cc b/NPSimulation/Process/PhysicsList.cc
index ee7b76530..2b58820af 100644
--- a/NPSimulation/Process/PhysicsList.cc
+++ b/NPSimulation/Process/PhysicsList.cc
@@ -23,6 +23,10 @@
//NPS
#include "PhysicsList.hh"
+//NPL
+#include "NPInputParser.h"
+#include "NPOptionManager.h"
+#include "NPNucleus.h"
// Local physic directly implemented from the Hadronthrapy example
// Physic dedicated to the ion-ion inelastic processes
@@ -44,7 +48,8 @@
#include "G4IonParametrisedLossModel.hh"
#include "G4UniversalFluctuation.hh"
#include "G4UImanager.hh"
-
+#include "G4NuclearLevelData.hh"
+#include "G4LevelManager.hh"
#include "G4PAIModel.hh"
#include "G4PAIPhotModel.hh"
#include "menate_R.hh"
@@ -55,227 +60,227 @@ PhysicsList::PhysicsList() : G4VUserPhysicsList(){
// UI commands to activate one or more High Energy Processes
// G4UImanager* UI = G4UImanager::GetUIpointer();
// UI->ApplyCommand("/physics_list/em/MuonNuclear true");
-
- m_EmList = "Option4";
- defaultCutValue = 1*mm;//0.2*mm;
- opticalPhysicsList = NULL;
- driftElectronPhysicsList = NULL;
- ReadConfiguration("PhysicsListOption.txt");
- G4LossTableManager::Instance();
- SetVerboseLevel(0);
-
- // ****** Definition of defaults for the physics processes *****
- // ****** in case no physics is called by the macro file *****
- //
- // The default physics corresponds to the actual QGSP_BIC_HP list
- // but with the following differences:
- // --> G4EmStandardPhysics_option4 for the electromagnetic processes
- // is used n place of the less accurate G4EmStandardPhysics
- // --> The G4RadioactiveDecayPhysics is add
- // --> G4HadronPhysicsQGSP_BIC is used in place of G4HadronPhysicsQGSP_BIC_HP
- // --> G4HadronElasticPhysics is used in place of G4HadronElasticPhysics_HP
-
- // Elecromagnetic physics
- // Using the more accurate option4
-
- emPhysicsList=NULL;
- if(m_EmList == "Option4"){
- emPhysicsList = new G4EmStandardPhysics_option4();
- cout << "//// Using G4EmStandardPhysics_option4 Physics List ////" << endl;
- }
- else if(m_EmList== "Option1"){
- emPhysicsList = new G4EmStandardPhysics_option1();
- cout << "//// Using G4EmStandardPhysics_option1 Physics List ////" << endl;
- }
- else if(m_EmList== "Option2"){
- emPhysicsList = new G4EmStandardPhysics_option2();
- cout << "//// Using G4EmStandardPhysics_option2 Physics List ////" << endl;
- }
- else if(m_EmList== "Option3"){
- emPhysicsList = new G4EmStandardPhysics_option3();
- cout << "//// Using G4EmStandardPhysics_option3 Physics List ////" << endl;
- }
- else if(m_EmList== "Standard"){
- emPhysicsList = new G4EmStandardPhysics();
- cout << "//// Using G4EmStandardPhysics default EM constructor Physics List ////" << endl;
- }
- else if(m_EmList== "Livermore"){
- emPhysicsList = new G4EmLivermorePhysics();
- cout << "//// Using G4EmLivermorePhysics Physics List ////" << endl;
- }
-
- else if( m_EmList == "Penelope"){
- emPhysicsList = new G4EmPenelopePhysics();
- cout << "//// Using G4EmPenelopePhysics Physics List ////" << endl;
- }
-
- else{
- std::cout << "\r\033[1;31mERROR: User given physics list " << m_EmList << " is not supported, option are Option4 Livermore Penelope\033[0m" <<std::endl;
- exit(1);
- }
- emConfig = G4LossTableManager::Instance()->EmConfigurator();
-
- // Hadronic physics
- if(m_IonBinaryCascadePhysics){
- m_PhysList["IonBinaryCascadePhysics"] = new G4IonBinaryCascadePhysics();
- cout << "//// Using G4IonBinaryCascadePhysics Physics List ////" << endl;
- }
- if(m_EmExtraPhysics)
- m_PhysList["EmExtraPhysics"]=new G4EmExtraPhysics();
-
- if(m_HadronElasticPhysics){
- m_PhysList["G4HadronElasticPhysics"]=new G4HadronElasticPhysics();
- m_PhysList["G4IonElasticPhysics"]=new G4IonElasticPhysics();
- cout << "//// Using G4HadronElasticPhysics Physics List ////" << endl;
- cout << "//// Using G4IonElasticPhysics Physics List ////" << endl;
- }
-
- if(m_NPIonInelasticPhysics){
- m_PhysList["NPIonInelasticPhysics"]=new NPIonIonInelasticPhysic();
- cout << "//// Using NPIonIonInelasticPhysic Physics List ////" << endl;
- }
-
- if(m_StoppingPhysics){
- m_PhysList["StoppingPhysics"]=new G4StoppingPhysics();
- cout << "//// Using G4StoppingPhysics Physics List ////" << endl;
- }
-
- if(m_HadronPhysicsINCLXX){
- m_PhysList["HadronPhysicsINCLXX"] = new G4HadronPhysicsINCLXX();
- m_PhysList["IonPhysicsINCLXX"] = new G4IonINCLXXPhysics();
- cout << "//// Using INCLXX Physics List ////" << endl;
- }
-
- if(m_HadronPhysicsQGSP_BIC_HP){
+
+ m_EmList = "Option4";
+ defaultCutValue = 1*mm;//0.2*mm;
+ opticalPhysicsList = NULL;
+ driftElectronPhysicsList = NULL;
+ ReadConfiguration("PhysicsListOption.txt");
+ G4LossTableManager::Instance();
+ SetVerboseLevel(0);
+
+ // ****** Definition of defaults for the physics processes *****
+ // ****** in case no physics is called by the macro file *****
+ //
+ // The default physics corresponds to the actual QGSP_BIC_HP list
+ // but with the following differences:
+ // --> G4EmStandardPhysics_option4 for the electromagnetic processes
+ // is used n place of the less accurate G4EmStandardPhysics
+ // --> The G4RadioactiveDecayPhysics is add
+ // --> G4HadronPhysicsQGSP_BIC is used in place of G4HadronPhysicsQGSP_BIC_HP
+ // --> G4HadronElasticPhysics is used in place of G4HadronElasticPhysics_HP
+
+ // Elecromagnetic physics
+ // Using the more accurate option4
+
+ emPhysicsList=NULL;
+ if(m_EmList == "Option4"){
+ emPhysicsList = new G4EmStandardPhysics_option4();
+ cout << "//// Using G4EmStandardPhysics_option4 Physics List ////" << endl;
+ }
+ else if(m_EmList== "Option1"){
+ emPhysicsList = new G4EmStandardPhysics_option1();
+ cout << "//// Using G4EmStandardPhysics_option1 Physics List ////" << endl;
+ }
+ else if(m_EmList== "Option2"){
+ emPhysicsList = new G4EmStandardPhysics_option2();
+ cout << "//// Using G4EmStandardPhysics_option2 Physics List ////" << endl;
+ }
+ else if(m_EmList== "Option3"){
+ emPhysicsList = new G4EmStandardPhysics_option3();
+ cout << "//// Using G4EmStandardPhysics_option3 Physics List ////" << endl;
+ }
+ else if(m_EmList== "Standard"){
+ emPhysicsList = new G4EmStandardPhysics();
+ cout << "//// Using G4EmStandardPhysics default EM constructor Physics List ////" << endl;
+ }
+ else if(m_EmList== "Livermore"){
+ emPhysicsList = new G4EmLivermorePhysics();
+ cout << "//// Using G4EmLivermorePhysics Physics List ////" << endl;
+ }
+
+ else if( m_EmList == "Penelope"){
+ emPhysicsList = new G4EmPenelopePhysics();
+ cout << "//// Using G4EmPenelopePhysics Physics List ////" << endl;
+ }
+
+ else{
+ std::cout << "\r\033[1;31mERROR: User given physics list " << m_EmList << " is not supported, option are Option4 Livermore Penelope\033[0m" <<std::endl;
+ exit(1);
+ }
+ emConfig = G4LossTableManager::Instance()->EmConfigurator();
+
+ // Hadronic physics
+ if(m_IonBinaryCascadePhysics){
+ m_PhysList["IonBinaryCascadePhysics"] = new G4IonBinaryCascadePhysics();
+ cout << "//// Using G4IonBinaryCascadePhysics Physics List ////" << endl;
+ }
+ if(m_EmExtraPhysics)
+ m_PhysList["EmExtraPhysics"]=new G4EmExtraPhysics();
+
+ if(m_HadronElasticPhysics){
+ m_PhysList["G4HadronElasticPhysics"]=new G4HadronElasticPhysics();
+ m_PhysList["G4IonElasticPhysics"]=new G4IonElasticPhysics();
+ cout << "//// Using G4HadronElasticPhysics Physics List ////" << endl;
+ cout << "//// Using G4IonElasticPhysics Physics List ////" << endl;
+ }
+
+ if(m_NPIonInelasticPhysics){
+ m_PhysList["NPIonInelasticPhysics"]=new NPIonIonInelasticPhysic();
+ cout << "//// Using NPIonIonInelasticPhysic Physics List ////" << endl;
+ }
+
+ if(m_StoppingPhysics){
+ m_PhysList["StoppingPhysics"]=new G4StoppingPhysics();
+ cout << "//// Using G4StoppingPhysics Physics List ////" << endl;
+ }
+
+ if(m_HadronPhysicsINCLXX){
+ m_PhysList["HadronPhysicsINCLXX"] = new G4HadronPhysicsINCLXX();
+ m_PhysList["IonPhysicsINCLXX"] = new G4IonINCLXXPhysics();
+ cout << "//// Using INCLXX Physics List ////" << endl;
+ }
+
+ if(m_HadronPhysicsQGSP_BIC_HP){
#if NPS_GEANT4_VERSION_MAJOR > 9
- m_PhysList["HadronPhysicsQGSP_BIC_HP"] = new G4HadronPhysicsQGSP_BIC_HP();
- cout << "//// Using QGSP_BIC_HP Physics List ////" << endl;
-
+ m_PhysList["HadronPhysicsQGSP_BIC_HP"] = new G4HadronPhysicsQGSP_BIC_HP();
+ cout << "//// Using QGSP_BIC_HP Physics List ////" << endl;
+
#else
- std::cout << "\r\032[1;31m Warning: physics list HadronPhysicsQGSP_BIC_HP require Geant4 10, process not activated \033[0m" <<std::endl;
+ std::cout << "\r\032[1;31m Warning: physics list HadronPhysicsQGSP_BIC_HP require Geant4 10, process not activated \033[0m" <<std::endl;
#endif
- }
- // Optical Photon for scintillator simulation
- if(m_OpticalPhysics){
- opticalPhysicsList = new G4OpticalPhysics(0);
- opticalPhysicsList->SetMaxNumPhotonsPerStep(100);
- opticalPhysicsList->SetScintillationYieldFactor(0.1);
- opticalPhysicsList->SetTrackSecondariesFirst(kScintillation,true);
- opticalPhysicsList->SetTrackSecondariesFirst(kCerenkov,true);
-
- cout << "//// Using Optical Photon Physics List ////" << endl;
- }
-
- // Drift electron for gazeous detector simulation
- if(m_DriftElectronPhysics){
- driftElectronPhysicsList = new G4DriftElectronPhysics(0);
- driftElectronPhysicsList->SetMaxNumDriftElectronPerStep(1e6);
- cout << "//// Using Drift Electron Physics List ////" << endl;
- }
- if(m_IonGasModels){
- AddIonGasModels();
- cout << "//// Using Ion Gas Model Physics List ////" << endl;
- }
-
- if(m_pai){
- AddPAIModel("pai");
- cout << "//// Using PAI Model Physics List ////" << endl;
- }
- if(m_pai_photon){
- AddPAIModel("pai_photon");
- cout << "//// Using PAI PHOTON Model Physics List ////" << endl;
- }
-
- // Decay physics
- // Add Radioactive decay
- // Gamma decay of known states
- if(m_Decay){
- decay_List = new G4DecayPhysics();
- radioactiveDecay_List = new G4RadioactiveDecayPhysics() ;
- m_PhysList["decay_list"]= decay_List;
- m_PhysList["radioactiveDecay"]=radioactiveDecay_List;
- }
-
- else{
- decay_List = 0;
- radioactiveDecay_List = 0;
- }
-
+ }
+ // Optical Photon for scintillator simulation
+ if(m_OpticalPhysics){
+ opticalPhysicsList = new G4OpticalPhysics(0);
+ opticalPhysicsList->SetMaxNumPhotonsPerStep(100);
+ opticalPhysicsList->SetScintillationYieldFactor(0.1);
+ opticalPhysicsList->SetTrackSecondariesFirst(kScintillation,true);
+ opticalPhysicsList->SetTrackSecondariesFirst(kCerenkov,true);
+
+ cout << "//// Using Optical Photon Physics List ////" << endl;
+ }
+
+ // Drift electron for gazeous detector simulation
+ if(m_DriftElectronPhysics){
+ driftElectronPhysicsList = new G4DriftElectronPhysics(0);
+ driftElectronPhysicsList->SetMaxNumDriftElectronPerStep(1e6);
+ cout << "//// Using Drift Electron Physics List ////" << endl;
+ }
+ if(m_IonGasModels){
+ AddIonGasModels();
+ cout << "//// Using Ion Gas Model Physics List ////" << endl;
+ }
+
+ if(m_pai){
+ AddPAIModel("pai");
+ cout << "//// Using PAI Model Physics List ////" << endl;
+ }
+ if(m_pai_photon){
+ AddPAIModel("pai_photon");
+ cout << "//// Using PAI PHOTON Model Physics List ////" << endl;
+ }
+
+ // Decay physics
+ // Add Radioactive decay
+ // Gamma decay of known states
+ if(m_Decay){
+ decay_List = new G4DecayPhysics();
+ radioactiveDecay_List = new G4RadioactiveDecayPhysics() ;
+ m_PhysList["decay_list"]= decay_List;
+ m_PhysList["radioactiveDecay"]=radioactiveDecay_List;
+ }
+
+ else{
+ decay_List = 0;
+ radioactiveDecay_List = 0;
+ }
+
}
////////////////////////////////////////////////////////////////////////////////
void PhysicsList::ReadConfiguration(std::string filename){
- m_IonBinaryCascadePhysics = 0;
- m_EmExtraPhysics = 0;
- m_HadronElasticPhysics = 0;
- m_NPIonInelasticPhysics = 0;
- m_StoppingPhysics = 0;
- m_OpticalPhysics = 0;
- m_DriftElectronPhysics = 0;
- m_HadronPhysicsQGSP_BIC_HP = 0;
- m_HadronPhysicsINCLXX = 0;
- m_Decay = 0;
- m_IonGasModels = 0;
- m_pai= 0;
- m_pai_photon= 0;
- m_Menate_R = 0;
-
- std::ifstream file(filename.c_str());
- if(!file.is_open()){
- std::cout << "WARNING: Could not find Physics List option file " << filename << " | Using default Physics List" << std::endl;
- return;
- }
-
- std::cout << "Reading Physics list option file " << filename << std::endl;
-
- std::string name;
- std::string str_value;
- double value;
- while(file >> name >> str_value){
- value = std::atof(str_value.c_str());
- if(name == "EmPhysicsList")
- m_EmList = str_value;
- else if(name == "DefaultCutOff")
- defaultCutValue = value*mm;
- else if(name == "IonBinaryCascadePhysics")
- m_IonBinaryCascadePhysics = value;
- else if(name == "NPIonInelasticPhysics")
- m_NPIonInelasticPhysics = value;
- else if (name == "EmExtraPhysics")
- m_EmExtraPhysics= value;
- else if (name == "HadronElasticPhysics")
- m_HadronElasticPhysics= value;
- else if (name == "StoppingPhysics")
- m_StoppingPhysics= value;
- else if (name == "OpticalPhysics")
- m_OpticalPhysics= value;
- else if (name == "DriftElectronPhysics")
- m_DriftElectronPhysics= value;
- else if (name == "HadronPhysicsQGSP_BIC_HP")
- m_HadronPhysicsQGSP_BIC_HP= value;
- else if (name == "HadronPhysicsINCLXX")
- m_HadronPhysicsINCLXX= value;
- else if (name == "Decay")
- m_Decay = value;
- else if (name == "IonGasModels")
- m_IonGasModels= value;
- else if (name == "pai")
- m_pai = value;
- else if (name == "pai_photon")
- m_pai_photon = value;
- else if (name == "Menate_R")
- m_Menate_R = value;
- else
- std::cout <<"WARNING: Physics List Token '" << name << "' unknown. Token is ignored." << std::endl;
- }
-
- // Most special process need decay to be activated
- if( (m_IonBinaryCascadePhysics || m_EmExtraPhysics || m_HadronElasticPhysics || m_NPIonInelasticPhysics
- || m_StoppingPhysics || m_HadronPhysicsQGSP_BIC_HP || m_HadronPhysicsINCLXX) && !m_Decay){
- m_Decay = true;
- std::cout << "Information: Selected process require Decay to be activated." << std::endl;
- }
-
+ m_IonBinaryCascadePhysics = 0;
+ m_EmExtraPhysics = 0;
+ m_HadronElasticPhysics = 0;
+ m_NPIonInelasticPhysics = 0;
+ m_StoppingPhysics = 0;
+ m_OpticalPhysics = 0;
+ m_DriftElectronPhysics = 0;
+ m_HadronPhysicsQGSP_BIC_HP = 0;
+ m_HadronPhysicsINCLXX = 0;
+ m_Decay = 0;
+ m_IonGasModels = 0;
+ m_pai= 0;
+ m_pai_photon= 0;
+ m_Menate_R = 0;
+
+ std::ifstream file(filename.c_str());
+ if(!file.is_open()){
+ std::cout << "WARNING: Could not find Physics List option file " << filename << " | Using default Physics List" << std::endl;
+ return;
+ }
+
+ std::cout << "Reading Physics list option file " << filename << std::endl;
+
+ std::string name;
+ std::string str_value;
+ double value;
+ while(file >> name >> str_value){
+ value = std::atof(str_value.c_str());
+ if(name == "EmPhysicsList")
+ m_EmList = str_value;
+ else if(name == "DefaultCutOff")
+ defaultCutValue = value*mm;
+ else if(name == "IonBinaryCascadePhysics")
+ m_IonBinaryCascadePhysics = value;
+ else if(name == "NPIonInelasticPhysics")
+ m_NPIonInelasticPhysics = value;
+ else if (name == "EmExtraPhysics")
+ m_EmExtraPhysics= value;
+ else if (name == "HadronElasticPhysics")
+ m_HadronElasticPhysics= value;
+ else if (name == "StoppingPhysics")
+ m_StoppingPhysics= value;
+ else if (name == "OpticalPhysics")
+ m_OpticalPhysics= value;
+ else if (name == "DriftElectronPhysics")
+ m_DriftElectronPhysics= value;
+ else if (name == "HadronPhysicsQGSP_BIC_HP")
+ m_HadronPhysicsQGSP_BIC_HP= value;
+ else if (name == "HadronPhysicsINCLXX")
+ m_HadronPhysicsINCLXX= value;
+ else if (name == "Decay")
+ m_Decay = value;
+ else if (name == "IonGasModels")
+ m_IonGasModels= value;
+ else if (name == "pai")
+ m_pai = value;
+ else if (name == "pai_photon")
+ m_pai_photon = value;
+ else if (name == "Menate_R")
+ m_Menate_R = value;
+ else
+ std::cout <<"WARNING: Physics List Token '" << name << "' unknown. Token is ignored." << std::endl;
+ }
+
+ // Most special process need decay to be activated
+ if( (m_IonBinaryCascadePhysics || m_EmExtraPhysics || m_HadronElasticPhysics || m_NPIonInelasticPhysics
+ || m_StoppingPhysics || m_HadronPhysicsQGSP_BIC_HP || m_HadronPhysicsINCLXX) && !m_Decay){
+ m_Decay = true;
+ std::cout << "Information: Selected process require Decay to be activated." << std::endl;
+ }
+
}
/////////////////////////////////////////////////////////////////////////////
PhysicsList::~PhysicsList(){
@@ -284,199 +289,230 @@ PhysicsList::~PhysicsList(){
/////////////////////////////////////////////////////////////////////////////
void PhysicsList::AddIonGasModels()
{
- G4ParticleTable::G4PTblDicIterator* particleIterator = G4ParticleTable::GetParticleTable()->GetIterator();
- particleIterator->reset();
- while ((*particleIterator)())
- {
- G4ParticleDefinition* particle = particleIterator->value();
- G4String partname = particle->GetParticleName();
- if(partname == "p" || partname == "proton" || partname == "alpha" || partname == "He3" || partname == "GenericIon") {
- G4BraggIonGasModel* mod1 = new G4BraggIonGasModel();
- G4BetheBlochIonGasModel* mod2 = new G4BetheBlochIonGasModel();
- G4double eth = 2.*MeV*particle->GetPDGMass()/proton_mass_c2;
- emConfig->SetExtraEmModel(partname,"ionIoni",mod1,"",0.0,eth,new G4IonFluctuations());
- emConfig->SetExtraEmModel(partname,"ionIoni",mod2,"",eth,100*TeV,new G4UniversalFluctuation());
- }
- }
+ G4ParticleTable::G4PTblDicIterator* particleIterator = G4ParticleTable::GetParticleTable()->GetIterator();
+ particleIterator->reset();
+ while ((*particleIterator)())
+ {
+ G4ParticleDefinition* particle = particleIterator->value();
+ G4String partname = particle->GetParticleName();
+ if(partname == "p" || partname == "proton" || partname == "alpha" || partname == "He3" || partname == "GenericIon") {
+ G4BraggIonGasModel* mod1 = new G4BraggIonGasModel();
+ G4BetheBlochIonGasModel* mod2 = new G4BetheBlochIonGasModel();
+ G4double eth = 2.*MeV*particle->GetPDGMass()/proton_mass_c2;
+ emConfig->SetExtraEmModel(partname,"ionIoni",mod1,"",0.0,eth,new G4IonFluctuations());
+ emConfig->SetExtraEmModel(partname,"ionIoni",mod2,"",eth,100*TeV,new G4UniversalFluctuation());
+ }
+ }
}
/////////////////////////////////////////////////////////////////////////////
void PhysicsList::AddPAIModel(const G4String& modname)
{
- G4ParticleTable::G4PTblDicIterator* particleIterator = G4ParticleTable::GetParticleTable()->GetIterator();
- particleIterator->reset();
- while ((*particleIterator)())
- {
- G4ParticleDefinition* particle = particleIterator->value();
- G4String partname = particle->GetParticleName();
- if(partname == "e-" || partname == "e+") {
- NewPAIModel(particle, modname, "eIoni");
-
- } else if(partname == "mu-" || partname == "mu+") {
- NewPAIModel(particle, modname, "muIoni");
-
- } else if(partname == "proton" ||
- partname == "pi+" ||
- partname == "pi-"
- ) {
- NewPAIModel(particle, modname, "hIoni");
- }
+ G4ParticleTable::G4PTblDicIterator* particleIterator = G4ParticleTable::GetParticleTable()->GetIterator();
+ particleIterator->reset();
+ while ((*particleIterator)())
+ {
+ G4ParticleDefinition* particle = particleIterator->value();
+ G4String partname = particle->GetParticleName();
+ if(partname == "e-" || partname == "e+") {
+ NewPAIModel(particle, modname, "eIoni");
+
+ } else if(partname == "mu-" || partname == "mu+") {
+ NewPAIModel(particle, modname, "muIoni");
+
+ } else if(partname == "proton" ||
+ partname == "pi+" ||
+ partname == "pi-"
+ ) {
+ NewPAIModel(particle, modname, "hIoni");
}
+ }
}
/////////////////////////////////////////////////////////////////////////////
void PhysicsList::NewPAIModel(const G4ParticleDefinition* part,
- const G4String& modname,
- const G4String& procname)
+ const G4String& modname,
+ const G4String& procname)
{
- G4String partname = part->GetParticleName();
- if(modname == "pai"){
- G4PAIModel* pai = new G4PAIModel(part,"PAIModel");
- emConfig->SetExtraEmModel(partname,procname,pai,"GasDetector",0.0,100.*TeV,pai);
- }
-
- else if(modname == "pai_photon") {
- G4PAIPhotModel* pai = new G4PAIPhotModel(part,"PAIPhotModel");
- emConfig->SetExtraEmModel(partname,procname,pai,"GasDetector",0.0,100.*TeV,pai);
- }
+ G4String partname = part->GetParticleName();
+ if(modname == "pai"){
+ G4PAIModel* pai = new G4PAIModel(part,"PAIModel");
+ emConfig->SetExtraEmModel(partname,procname,pai,"GasDetector",0.0,100.*TeV,pai);
+ }
+
+ else if(modname == "pai_photon") {
+ G4PAIPhotModel* pai = new G4PAIPhotModel(part,"PAIPhotModel");
+ emConfig->SetExtraEmModel(partname,procname,pai,"GasDetector",0.0,100.*TeV,pai);
+ }
}
/////////////////////////////////////////////////////////////////////////////
void PhysicsList::ConstructParticle(){
- if(m_OpticalPhysics){
- //G4OpticalPhoton::OpticalPhotonDefinition();
- ((G4VPhysicsConstructor*) opticalPhysicsList)->ConstructParticle();
-
- }
-
- if(m_DriftElectronPhysics){
- ((G4VPhysicsConstructor*) driftElectronPhysicsList)->ConstructParticle();
- }
-
- if(decay_List){
- decay_List -> ConstructParticle();
- radioactiveDecay_List->ConstructParticle();
- }
-
- else{
- // If decay is not activated we have to declare the particle ourself
- G4He3::He3Definition();
- G4Deuteron::DeuteronDefinition();
- G4Triton::TritonDefinition();
- G4Alpha::AlphaDefinition();
- G4Proton::ProtonDefinition();
- G4AntiProton::AntiProtonDefinition();
- G4Neutron::NeutronDefinition();
- G4Electron::ElectronDefinition();
- G4Positron::PositronDefinition();
- G4NeutrinoE::NeutrinoEDefinition();
- G4MuonPlus::MuonPlusDefinition();
- G4MuonMinus::MuonMinusDefinition();
- G4AntiNeutrinoE::AntiNeutrinoEDefinition();
- G4Geantino::GeantinoDefinition();
- G4ChargedGeantino::ChargedGeantinoDefinition();
- G4Gamma::GammaDefinition();
- // mesons
- G4PionPlus ::PionPlusDefinition();
- G4PionMinus ::PionMinusDefinition();
- G4PionZero ::PionZeroDefinition();
- G4Eta ::EtaDefinition();
- G4EtaPrime ::EtaPrimeDefinition();
- // G4RhoZero ::RhoZeroDefinition();
- G4KaonPlus ::KaonPlusDefinition();
- G4KaonMinus ::KaonMinusDefinition();
- G4KaonZero ::KaonZeroDefinition();
- G4AntiKaonZero ::AntiKaonZeroDefinition();
- G4KaonZeroLong ::KaonZeroLongDefinition();
- G4KaonZeroShort::KaonZeroShortDefinition();
- // Ion
- G4IonConstructor ionConstructor ;
- ionConstructor.ConstructParticle() ;
- }
+ if(m_OpticalPhysics){
+ //G4OpticalPhoton::OpticalPhotonDefinition();
+ ((G4VPhysicsConstructor*) opticalPhysicsList)->ConstructParticle();
+
+ }
+
+ if(m_DriftElectronPhysics){
+ ((G4VPhysicsConstructor*) driftElectronPhysicsList)->ConstructParticle();
+ }
+
+ if(decay_List){
+ decay_List -> ConstructParticle();
+ radioactiveDecay_List->ConstructParticle();
+ }
+
+ else{
+ // If decay is not activated we have to declare the particle ourself
+ G4He3::He3Definition();
+ G4Deuteron::DeuteronDefinition();
+ G4Triton::TritonDefinition();
+ G4Alpha::AlphaDefinition();
+ G4Proton::ProtonDefinition();
+ G4AntiProton::AntiProtonDefinition();
+ G4Neutron::NeutronDefinition();
+ G4Electron::ElectronDefinition();
+ G4Positron::PositronDefinition();
+ G4NeutrinoE::NeutrinoEDefinition();
+ G4MuonPlus::MuonPlusDefinition();
+ G4MuonMinus::MuonMinusDefinition();
+ G4AntiNeutrinoE::AntiNeutrinoEDefinition();
+ G4Geantino::GeantinoDefinition();
+ G4ChargedGeantino::ChargedGeantinoDefinition();
+ G4Gamma::GammaDefinition();
+ // mesons
+ G4PionPlus ::PionPlusDefinition();
+ G4PionMinus ::PionMinusDefinition();
+ G4PionZero ::PionZeroDefinition();
+ G4Eta ::EtaDefinition();
+ G4EtaPrime ::EtaPrimeDefinition();
+ // G4RhoZero ::RhoZeroDefinition();
+ G4KaonPlus ::KaonPlusDefinition();
+ G4KaonMinus ::KaonMinusDefinition();
+ G4KaonZero ::KaonZeroDefinition();
+ G4AntiKaonZero ::AntiKaonZeroDefinition();
+ G4KaonZeroLong ::KaonZeroLongDefinition();
+ G4KaonZeroShort::KaonZeroShortDefinition();
+ // Ion
+ G4IonConstructor ionConstructor ;
+ ionConstructor.ConstructParticle() ;
+ }
}
/////////////////////////////////////////////////////////////////////////////
void PhysicsList::ConstructProcess(){
- // Transportation
- AddTransportation();
-
- // Electromagnetic physics
- emPhysicsList -> ConstructProcess();
- if(opticalPhysicsList){
- ((G4VPhysicsConstructor*) opticalPhysicsList)->ConstructProcess();
- }
- if(driftElectronPhysicsList){
- ((G4VPhysicsConstructor*) driftElectronPhysicsList)->ConstructProcess();
- }
-
- // Hadronic physics
- std::map<std::string,G4VPhysicsConstructor*>::iterator it;
- for(it = m_PhysList.begin(); it!= m_PhysList.end(); it++){
- it->second -> ConstructProcess();
- }
- BiasCrossSectionByFactor(m_IonBinaryCascadePhysics);
-
- em_option.SetBuildCSDARange(true);
- em_option.SetDEDXBinningForCSDARange(10*10);
- em_option.SetFluo(true);
- em_option.SetAuger(true);
-
- AddParametrisation();
-
- return;
+ // Transportation
+ AddTransportation();
+
+ // Electromagnetic physics
+ emPhysicsList -> ConstructProcess();
+ if(opticalPhysicsList){
+ ((G4VPhysicsConstructor*) opticalPhysicsList)->ConstructProcess();
+ }
+ if(driftElectronPhysicsList){
+ ((G4VPhysicsConstructor*) driftElectronPhysicsList)->ConstructProcess();
+ }
+
+ // Hadronic physics
+ std::map<std::string,G4VPhysicsConstructor*>::iterator it;
+ for(it = m_PhysList.begin(); it!= m_PhysList.end(); it++){
+ it->second -> ConstructProcess();
+ }
+ BiasCrossSectionByFactor(m_IonBinaryCascadePhysics);
+
+ em_option.SetBuildCSDARange(true);
+ em_option.SetDEDXBinningForCSDARange(10*10);
+ em_option.SetFluo(true);
+ em_option.SetAuger(true);
+
+ AddParametrisation();
+ AddLevelData();
+ return;
}
/////////////////////////////////////////////////////////////////////////////
void PhysicsList::AddStepMax(){
}
/////////////////////////////////////////////////////////////////////////////
void PhysicsList::AddParametrisation(){
- G4FastSimulationManagerProcess* BeamReaction =
+ G4FastSimulationManagerProcess* BeamReaction =
new G4FastSimulationManagerProcess("NPSimulationProcess");
-
- // For 10.3 and higher
+
+ // For 10.3 and higher
#ifndef theParticleIterator
- G4ParticleTable::G4PTblDicIterator* theParticleIterator = GetParticleIterator();
+ G4ParticleTable::G4PTblDicIterator* theParticleIterator = GetParticleIterator();
#endif
-
- theParticleIterator->reset();
- while ((*theParticleIterator)()){
- G4ParticleDefinition* particle = theParticleIterator->value();
- G4ProcessManager* pmanager = particle->GetProcessManager();
- std::string name = particle->GetParticleName();
- pmanager->AddDiscreteProcess(BeamReaction);
- // Add a Step limiter to the beam particle.
- // This will be used to limit the step of the beam in the target
- pmanager->AddProcess(new G4StepLimiter,-1,-1,5);
-
- // Add menate R to the process manager for neutrons
- // (if selected as a option)
- if(m_Menate_R > 0 && name == "neutron") {
- menate_R* theMenate = new menate_R("menateR_neutron");
- theMenate->SetMeanFreePathCalcMethod("ORIGINAL");
- pmanager->AddDiscreteProcess(theMenate);
- std::cout <<"||--------------------------------------------------||" << std::endl;
- std::cout <<" MENATE_R Added to Process Manager " << std::endl;
- std::cout <<"||--------------------------------------------------||" << std::endl;
- }
- }
+
+ theParticleIterator->reset();
+ while ((*theParticleIterator)()){
+ G4ParticleDefinition* particle = theParticleIterator->value();
+ G4ProcessManager* pmanager = particle->GetProcessManager();
+ std::string name = particle->GetParticleName();
+ pmanager->AddDiscreteProcess(BeamReaction);
+ // Add a Step limiter to the beam particle.
+ // This will be used to limit the step of the beam in the target
+ pmanager->AddProcess(new G4StepLimiter,-1,-1,5);
+
+ // Add menate R to the process manager for neutrons
+ // (if selected as a option)
+ if(m_Menate_R > 0 && name == "neutron") {
+ menate_R* theMenate = new menate_R("menateR_neutron");
+ theMenate->SetMeanFreePathCalcMethod("ORIGINAL");
+ pmanager->AddDiscreteProcess(theMenate);
+ std::cout <<"||--------------------------------------------------||" << std::endl;
+ std::cout <<" MENATE_R Added to Process Manager " << std::endl;
+ std::cout <<"||--------------------------------------------------||" << std::endl;
+ }
+ }
}
+////////////////////////////////////////////////////////////////////////////////
+void PhysicsList::AddLevelData(){
+ NPL::InputParser parser(NPOptionManager::getInstance()->GetReactionFile());
+ vector<NPL::InputBlock*> blocks;
+ blocks = parser.GetAllBlocksWithToken("LevelData");
+ // Loop over all blocks
+ if (blocks.size()>0 && NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << endl << "//// Reading user level data " << endl;
+
+ for (unsigned int i = 0; i < blocks.size(); i++) {
+ NPL::Nucleus N(blocks[i]->GetMainValue());
+ int Z =N.GetZ() ; int A=N.GetA();
+ cout << " Nuclei: " << N.GetName() << endl;
+ string level_path = blocks[i]->GetString("Path");
+ G4IonTable* ionTable = G4IonTable::GetIonTable();
+ G4NuclearLevelData* levelData = G4NuclearLevelData::GetInstance();
+ levelData->AddPrivateData(Z, A, level_path);
+ const G4LevelManager* levelManager = levelData->GetLevelManager(Z, A);
+ G4int Nentries = levelManager->NumberOfTransitions()+1;
+ for(G4int j = 1; j < Nentries; j++){ // Excited states
+ cout << " - Level " << j
+ << " energy = " << levelManager->LevelEnergy(j)
+ << endl;
+ G4ParticleDefinition* excitedState
+ = ionTable->GetIon(Z,A,levelManager->LevelEnergy(j));
+ excitedState->SetPDGStable(false);
+ excitedState->SetPDGLifeTime(levelManager->LifeTime(j));
+ }
+ }
+
+}
/////////////////////////////////////////////////////////////////////////////
void PhysicsList::SetCuts(){
-
- if (verboseLevel >0){
- G4cout << "PhysicsList::SetCuts:";
- G4cout << "CutLength : " << G4BestUnit(defaultCutValue,"Length") << G4endl;
- }
-
- // Special Cut for optical photon to be emmitted
- SetCutsWithDefault();
+
+ if (verboseLevel >0){
+ G4cout << "PhysicsList::SetCuts:";
+ G4cout << "CutLength : " << G4BestUnit(defaultCutValue,"Length") << G4endl;
+ }
+
+ // Special Cut for optical photon to be emmitted
+ SetCutsWithDefault();
}
////////////////////////////////////////////////////////////////////////////////
/////// Friend Method for CS biasing
////////////////////////////////////////////////////////////////////////////////
void PhysicsList::BiasCrossSectionByFactor(double factor){
- factor++;
+ factor++;
}
diff --git a/NPSimulation/Process/PhysicsList.hh b/NPSimulation/Process/PhysicsList.hh
index 33f7eb807..5a7b488b3 100644
--- a/NPSimulation/Process/PhysicsList.hh
+++ b/NPSimulation/Process/PhysicsList.hh
@@ -89,7 +89,7 @@ class PhysicsList: public G4VUserPhysicsList{
void AddIonGasModels();
void AddPAIModel(const G4String& modname);
void NewPAIModel(const G4ParticleDefinition* part, const G4String& modname,const G4String& procname);
-
+ void AddLevelData();
private:
std::map<std::string,G4VPhysicsConstructor*> m_PhysList;
G4EmConfigurator em_config;
diff --git a/Projects/e748/Analysis.cxx b/Projects/e748/Analysis.cxx
index e837582a9..376c1dd90 100644
--- a/Projects/e748/Analysis.cxx
+++ b/Projects/e748/Analysis.cxx
@@ -22,7 +22,6 @@
#include<iostream>
using namespace std;
#include"Analysis.h"
-#include"NPFunction.h"
#include"NPAnalysisFactory.h"
#include"NPDetectorManager.h"
#include"NPOptionManager.h"
@@ -60,16 +59,14 @@ void Analysis::Init(){
ZTarget = 0;
TimeCorr=0;
TargetThickness = m_DetectorManager->GetTargetThickness();
- string TargetMaterial = m_DetectorManager->GetTargetMaterial();
- // energy losses
- string light=NPL::ChangeNameToG4Standard(myReaction->GetNucleus3()->GetName());
- string beam=NPL::ChangeNameToG4Standard(myReaction->GetNucleus1()->GetName());
- LightTarget = NPL::EnergyLoss(light+"_"+TargetMaterial+".G4table","G4Table",10 );
- LightAl = NPL::EnergyLoss(light+"_Al.G4table","G4Table",10);
- LightSi = NPL::EnergyLoss(light+"_Si.G4table","G4Table",10);
- BeamTarget = NPL::EnergyLoss(beam+"_"+TargetMaterial+".G4table","G4Table",10);
- BeamMylar= NPL::EnergyLoss(beam+"_Mylar.G4table","G4Table",10);
- BeamIsobutane = EnergyLoss(beam+"_iC4H10.G4table","G4Table",10);
+ // Energy loss table: the G4Table are generated by the simulation
+ He3CD2 = EnergyLoss("Example/He3_CD2.G4table","G4Table",10 );
+ He3Al = EnergyLoss("Example/He3_Al.G4table","G4Table",10);
+ He3Si = EnergyLoss("Example/He3_Si.G4table","G4Table",10);
+ BeamCD2 = EnergyLoss("Be12_CD2.G4table","G4Table",10);
+ BeamMylar = EnergyLoss("Be12_Mylar.G4table","G4Table",10);
+ BeamIsobutane = EnergyLoss("Be12_iC4H10.G4table","G4Table",10);
+
}
////////////////////////////////////////////////////////////////////////////////
@@ -85,7 +82,7 @@ void Analysis::TreatEvent(){
filename = filename.substr(0,minor_pos);
size_t major_pos = filename.rfind("_");
run_major = atoi(filename.substr(major_pos+1,4).c_str());
- // Get the Init information on beam position and energy
+ // Get the Init information on beam position and energy
// and apply by hand the experimental resolution
// This is because the beam diagnosis are not simulated
// PPAC position resolution on target is assumed to be 1mm
@@ -96,124 +93,127 @@ void Analysis::TreatEvent(){
// Beam energy is measured using F3 and F2 plastic TOF
//double BeamEnergy= myReaction-> GetBeamEnergy()* MeV;
//////////////////////////// LOOP on MUST2 Hit //////////////////
- for(unsigned int countMust2 = 0 ; countMust2 < M2->Si_E.size() ; countMust2++){
+ for(unsigned int countMust2 = 0 ; countMust2 < M2->Si_E.size() ; countMust2++){
/* //Part 0 : Get the usefull Data */
- // MUST2 */
- int X = M2->Si_X[countMust2];
- int Y = M2->Si_Y[countMust2];
- int TelescopeNumber = M2->TelescopeNumber[countMust2];
- Si_X_M2 = X ;
- Si_Y_M2 = Y ;
- // Beam Energy from Cav Time of Flight //
- double BeamSpeed = 138.898 - 14765.1/ ModularLeaf->GetCalibratedValue("T_CATS1_CAV") ; // mm/ns
- // Beam Energy before CATS1
- static double c2 = 299.792458*299.792458;// mm/ns
- double gamma = 1./sqrt(1-BeamSpeed*BeamSpeed/c2);
- BeamEnergy= 11200.962140*(gamma-1);
- double BeamAngle= BeamDirection.Angle(TVector3(0,0,1));
- double gammaCav = (BeamEnergy+11200.962140) / 11200.962140 ;
- double BeamSpeedCav = sqrt(c2*(1-1/(gammaCav*gammaCav)));
-//cout << ModularLeaf->GetCalibratedValue("T_CATS1_CAV") << " " << BeamSpeed << " " << BeamEnergy << " " << BeamEnergy/12. << endl;
- if(BeamEnergy>0 && TelescopeNumber<5){
- DetectorNumber = TelescopeNumber ;
-
- /* // Part 1 : Impact Angle */
- ThetaM2Surface = 0;
- ThetaNormalTarget = 0;
- if(XTarget>-1000 && YTarget>-1000){
- TVector3 BeamImpact(XTarget,YTarget,0);
- TVector3 HitDirection = M2 -> GetPositionOfInteraction(countMust2) - BeamImpact ;
- ThetaLab = HitDirection.Angle( BeamDirection );
-
- ThetaM2Surface = HitDirection.Angle(- M2 -> GetTelescopeNormal(countMust2) );
- ThetaNormalTarget = HitDirection.Angle( TVector3(0,0,1) ) ;
- X_M2 = M2 -> GetPositionOfInteraction(countMust2).X() ;
- Y_M2 = M2 -> GetPositionOfInteraction(countMust2).Y() ;
- Z_M2 = M2 -> GetPositionOfInteraction(countMust2).Z() ;
-
-
+ // MUST2 */
+ int X = M2->Si_X[countMust2];
+ int Y = M2->Si_Y[countMust2];
+ int TelescopeNumber = M2->TelescopeNumber[countMust2];
+ Si_X_M2 = X ;
+ Si_Y_M2 = Y ;
+
+ if(TelescopeNumber<9){
+ DetectorNumber = TelescopeNumber ;
+
+ /* // Part 1 : Impact Angle */
+ ThetaM2Surface = 0;
+ ThetaNormalTarget = 0;
+ if(XTarget>-1000 && YTarget>-1000){
+ TVector3 BeamImpact(XTarget,YTarget,0);
+ TVector3 HitDirection = M2 -> GetPositionOfInteraction(countMust2) - BeamImpact ;
+ ThetaLab = HitDirection.Angle( BeamDirection );
+
+ ThetaM2Surface = HitDirection.Angle(- M2 -> GetTelescopeNormal(countMust2) );
+ ThetaNormalTarget = HitDirection.Angle( TVector3(0,0,1) ) ;
+ X_M2 = M2 -> GetPositionOfInteraction(countMust2).X() ;
+ Y_M2 = M2 -> GetPositionOfInteraction(countMust2).Y() ;
+ Z_M2 = M2 -> GetPositionOfInteraction(countMust2).Z() ;
+
+ // Beam Energy from Cav Time of Flight //
+
+ // Beam speed from Beam Energy
+
+ // double BeamSpeed = 10.8727 + ModularLeaf->GetCalibratedValue("T_CATS1_CAV")*0.276825; // mm/ns
+ //double BeamSpeed = 5.17952 + ModularLeaf->GetCalibratedValue("T_CATS1_CAV")*0.305315; // mm/ns
+ double BeamSpeed = 11.0476 + ModularLeaf->GetCalibratedValue("T_CATS1_CAV")*0.278917; // mm/ns
+
+ // Beam Energy before CATS1
+ static double c2 = 299.792458*299.792458;// mm/ns
+ double gamma = 1./sqrt(1-BeamSpeed*BeamSpeed/c2);
+ BeamEnergy= 11200.962140*(gamma-1);
+ double BeamAngle= BeamDirection.Angle(TVector3(0,0,1));
+
// Beam Energy and speed after CATS1
- double BeamEnergyC1 = BeamMylar.Slow(BeamEnergy,1.2*micrometer,BeamAngle);
- BeamEnergyC1 = BeamIsobutane.Slow(BeamEnergyC1,cm/3.,BeamAngle);
- BeamEnergyC1 = BeamMylar.Slow(BeamEnergyC1,0.9*micrometer,BeamAngle);
- BeamEnergyC1 = BeamIsobutane.Slow(BeamEnergyC1,cm/3,BeamAngle);
- BeamEnergyC1 = BeamMylar.Slow(BeamEnergyC1,0.9*micrometer,BeamAngle);
- BeamEnergyC1 = BeamIsobutane.Slow(BeamEnergyC1,cm/3.,BeamAngle);
- BeamEnergyC1 = BeamMylar.Slow(BeamEnergyC1,1.2*micrometer,BeamAngle);
- double gammaC1 = (BeamEnergyC1+11200.962140) / 11200.962140 ;
- double BeamSpeedC1 = sqrt(c2*(1-1/(gammaC1*gammaC1)));
- TVector3 C1toC2 = TVector3(CATS->PositionX[1],CATS->PositionY[1],CATS->PositionZ[1])
- - TVector3(CATS->PositionX[0],CATS->PositionY[0],CATS->PositionZ[0]) ;
- TimeCorr = C1toC2.Mag()/BeamSpeedC1;
-
- // Beam Energy and speed after CATS2
- double BeamEnergyC2 = BeamMylar.Slow(BeamEnergyC1,1.2*micrometer,BeamAngle);
- BeamEnergyC2 = BeamIsobutane.Slow(BeamEnergyC2,cm/3.,BeamAngle);
- BeamEnergyC2 = BeamMylar.Slow(BeamEnergyC2,0.9*micrometer,BeamAngle);
- BeamEnergyC2 = BeamIsobutane.Slow(BeamEnergyC2,cm/3,BeamAngle);
- BeamEnergyC2 = BeamMylar.Slow(BeamEnergyC2,0.9*micrometer,BeamAngle);
- BeamEnergyC2 = BeamIsobutane.Slow(BeamEnergyC2,cm/3.,BeamAngle);
- BeamEnergyC2 = BeamMylar.Slow(BeamEnergyC2,1.2*micrometer,BeamAngle);
-
- double gammaC2 = (BeamEnergyC2+11200.962140) / 11200.962140;
- double BeamSpeedC2 = sqrt(c2*(1-1/(gammaC2*gammaC2)));
- TVector3 C2toTarget = BeamImpact-TVector3(CATS->PositionX[1],CATS->PositionY[1],CATS->PositionZ[1]);
- TimeCorr += C2toTarget.Mag()/BeamSpeedC2;
- // slow down beam inside the target
- BeamEnergy = BeamTarget.Slow(BeamEnergyC2,TargetThickness*0.5,BeamDirection.Angle(TVector3(0,0,1)));
- myReaction->SetBeamEnergy(BeamEnergy);
-
-
- ParticleLength=HitDirection.Mag();
-
- }
-
-
- else{
- BeamDirection = TVector3(-1000,-1000,-1000);
- ThetaM2Surface = -1000 ;
- ThetaNormalTarget = -1000 ;
- }
-
- /* // Part 2 : Impact Energy */
- Energy = ELab = E_M2 = 0;
- Si_E_M2 = M2->Si_E[countMust2];
- CsI_E_M2= M2->CsI_E[countMust2];
-
- /* // if CsI */
- /* /1* if(CsI_E_M2>0 ){ *1/ */
- /* /1* // The energy in CsI is calculate form dE/dx Table because *1/ */
- /* /1* // 20um resolution is poor *1/ */
- /* /1* Energy = *1/ */
- /* /1* LightSi.EvaluateEnergyFromDeltaE(Si_E_M2,300*micrometer, *1/ */
- /* /1* ThetaM2Surface, 0.01*MeV, *1/ */
- /* /1* 450.*MeV,0.001*MeV ,1000); *1/ */
- /* /1* E_M2=CsI_E_M2; *1/ */
- /* /1* } *1/ */
-
- /* /1* else *1/ */
-
-
-
-
- Energy = Si_E_M2;
-
- E_M2 += Si_E_M2;
-
- /* // Evaluate energy using the thickness */
- ELab = LightAl.EvaluateInitialEnergy( Energy,0.4*micrometer , ThetaM2Surface);
- /* // Target Correction */
- ELab = LightTarget.EvaluateInitialEnergy( ELab ,TargetThickness/2., ThetaNormalTarget);
-
- /* // Part 3 : Excitation Energy Calculation */
- Ex = myReaction -> ReconstructRelativistic( ELab , ThetaLab );
-
-
- /* // Part 4 : Theta CM Calculation */
- ThetaCM = myReaction -> EnergyLabToThetaCM( ELab , ThetaLab)/deg;
- ThetaLab=ThetaLab/deg;
-
- }
+ double BeamEnergyC1 = BeamMylar.Slow(BeamEnergy,1.2*micrometer,BeamAngle);
+ BeamEnergyC1 = BeamIsobutane.Slow(BeamEnergyC1,cm/3.,BeamAngle);
+ BeamEnergyC1 = BeamMylar.Slow(BeamEnergyC1,0.9*micrometer,BeamAngle);
+ BeamEnergyC1 = BeamIsobutane.Slow(BeamEnergyC1,cm/3,BeamAngle);
+ BeamEnergyC1 = BeamMylar.Slow(BeamEnergyC1,0.9*micrometer,BeamAngle);
+ BeamEnergyC1 = BeamIsobutane.Slow(BeamEnergyC1,cm/3.,BeamAngle);
+ BeamEnergyC1 = BeamMylar.Slow(BeamEnergyC1,1.2*micrometer,BeamAngle);
+
+ double gammaC1 = (BeamEnergyC1+11200.962140) / 11200.962140 ;
+ double BeamSpeedC1 = sqrt(c2*(1-1/(gammaC1*gammaC1)));
+ TVector3 C1toC2 = TVector3(CATS->PositionX[1],CATS->PositionY[1],CATS->PositionZ[1])
+ - TVector3(CATS->PositionX[0],CATS->PositionY[0],CATS->PositionZ[0]) ;
+ TimeCorr = C1toC2.Mag()/BeamSpeedC1;
+
+ // Beam Energy and speed after CATS2
+ double BeamEnergyC2 = BeamMylar.Slow(BeamEnergyC1,1.2*micrometer,BeamAngle);
+ BeamEnergyC2 = BeamIsobutane.Slow(BeamEnergyC2,cm/3.,BeamAngle);
+ BeamEnergyC2 = BeamMylar.Slow(BeamEnergyC2,0.9*micrometer,BeamAngle);
+ BeamEnergyC2 = BeamIsobutane.Slow(BeamEnergyC2,cm/3,BeamAngle);
+ BeamEnergyC2 = BeamMylar.Slow(BeamEnergyC2,0.9*micrometer,BeamAngle);
+ BeamEnergyC2 = BeamIsobutane.Slow(BeamEnergyC2,cm/3.,BeamAngle);
+ BeamEnergyC2 = BeamMylar.Slow(BeamEnergyC2,1.2*micrometer,BeamAngle);
+
+ double gammaC2 = (BeamEnergyC2+11200.962140) / 11200.962140;
+ double BeamSpeedC2 = sqrt(c2*(1-1/(gammaC2*gammaC2)));
+ TVector3 C2toTarget = BeamImpact-TVector3(CATS->PositionX[1],CATS->PositionY[1],CATS->PositionZ[1]);
+ TimeCorr += C2toTarget.Mag()/BeamSpeedC2;
+
+ // slow down beam inside the target
+ BeamEnergy = BeamCD2.Slow(BeamEnergyC2,TargetThickness*0.5,BeamDirection.Angle(TVector3(0,0,1)));
+ myReaction->SetBeamEnergy(BeamEnergy);
+
+ }
+
+
+ else{
+ BeamDirection = TVector3(-1000,-1000,-1000);
+ ThetaM2Surface = -1000 ;
+ ThetaNormalTarget = -1000 ;
+ }
+
+/* // Part 2 : Impact Energy */
+ Energy = ELab = E_M2 = 0;
+ Si_E_M2 = M2->Si_E[countMust2];
+ CsI_E_M2= M2->CsI_E[countMust2];
+
+/* // if CsI */
+/* /1* if(CsI_E_M2>0 ){ *1/ */
+/* /1* // The energy in CsI is calculate form dE/dx Table because *1/ */
+/* /1* // 20um resolution is poor *1/ */
+/* /1* Energy = *1/ */
+/* /1* He3Si.EvaluateEnergyFromDeltaE(Si_E_M2,300*micrometer, *1/ */
+/* /1* ThetaM2Surface, 0.01*MeV, *1/ */
+/* /1* 450.*MeV,0.001*MeV ,1000); *1/ */
+/* /1* E_M2=CsI_E_M2; *1/ */
+/* /1* } *1/ */
+
+/* /1* else *1/ */
+
+
+
+
+ Energy = Si_E_M2;
+
+ E_M2 += Si_E_M2;
+
+/* // Evaluate energy using the thickness */
+ ELab = He3Al.EvaluateInitialEnergy( Energy,0.4*micrometer , ThetaM2Surface);
+/* // Target Correction */
+ ELab = He3CD2.EvaluateInitialEnergy( ELab ,TargetThickness/2., ThetaNormalTarget);
+
+/* // Part 3 : Excitation Energy Calculation */
+ Ex = myReaction -> ReconstructRelativistic( ELab , ThetaLab );
+
+
+/* // Part 4 : Theta CM Calculation */
+ ThetaCM = myReaction -> EnergyLabToThetaCM( ELab , ThetaLab)/deg;
+ ThetaLab=ThetaLab/deg;
+
+ }
}//end loop MUST2
}
@@ -240,35 +240,35 @@ void Analysis::InitOutputBranch() {
RootOutput::getInstance()->GetTree()->Branch("RunMajor",&run_major,"RunMajor/I");
RootOutput::getInstance()->GetTree()->Branch("RunMinor",&run_minor,"RunMinor/I");
- /* RootOutput::getInstance()->GetTree()->Branch("ADC_CHIO_V15",&vADC_CHIO_V15,"ADC_CHIO_V15/s");
- RootOutput::getInstance()->GetTree()->Branch("ADC_VOIE_29",&vADC_VOIE_29,"ADC_VOIE_29/s");
- RootOutput::getInstance()->GetTree()->Branch("CHIO",&vCHIO,"CHIO/s");
- RootOutput::getInstance()->GetTree()->Branch("CONFDEC",&vCONFDEC,"CONFDEC/s");
- RootOutput::getInstance()->GetTree()->Branch("CONFDEC_AGAVA",&vCONFDEC_AGAVA,"CONFDEC_AGAVA/s");
- RootOutput::getInstance()->GetTree()->Branch("DATATRIG",&vDATATRIG,"DATATRIG/s");
- RootOutput::getInstance()->GetTree()->Branch("DATATRIG_CHIO",&vDATATRIG_CHIO,"DATATRIG_CHIO/s");
- RootOutput::getInstance()->GetTree()->Branch("E1D6",&vE1D6,"E1D6/s");
- RootOutput::getInstance()->GetTree()->Branch("E2D6",&vE2D6,"E2D6/s");
- RootOutput::getInstance()->GetTree()->Branch("ED4",&vED4,"ED4/s");
- RootOutput::getInstance()->GetTree()->Branch("EXL_HF",&vEXL_HF,"EXL_HF/s");
- RootOutput::getInstance()->GetTree()->Branch("GALD4X",&vGALD4X,"GALD4X/s");
- RootOutput::getInstance()->GetTree()->Branch("GALD4Y",&vGALD4Y,"GALD4Y/s");
- RootOutput::getInstance()->GetTree()->Branch("QCaviar",&vQCaviar,"QCaviar/s");
- RootOutput::getInstance()->GetTree()->Branch("QPlast",&vQPlast,"QPlast/s");
- RootOutput::getInstance()->GetTree()->Branch("TCAVHF",&vTCAVHF,"TCAVHF/s");
- RootOutput::getInstance()->GetTree()->Branch("TE1D6CAV",&vTE1D6CAV,"TE1D6CAV/s");
- RootOutput::getInstance()->GetTree()->Branch("TE1D6GAL",&vTE1D6GAL,"TE1D6GAL/s");
- RootOutput::getInstance()->GetTree()->Branch("TE1D6HF",&vTE1D6HF,"TE1D6HF/s");
- RootOutput::getInstance()->GetTree()->Branch("TED4HF",&vTED4HF,"TED4HF/s");
- RootOutput::getInstance()->GetTree()->Branch("TGALD4HF",&vTGALD4HF,"TGALD4HF/s");
- RootOutput::getInstance()->GetTree()->Branch("T_CATS1_2",&vT_CATS1_2,"T_CATS1_2/s");
- RootOutput::getInstance()->GetTree()->Branch("T_CATS1_CAV",&vT_CATS1_CAV,"T_CATS1_CAV/s");
- RootOutput::getInstance()->GetTree()->Branch("T_MUVI_CATS1",&vT_MUVI_CATS1,"T_MUVI_CATS1/s");
- RootOutput::getInstance()->GetTree()->Branch("T_PL_CATS1",&vT_PL_CATS1,"T_PL_CATS1/s");
- RootOutput::getInstance()->GetTree()->Branch("T_PL_CATS2",&vT_PL_CATS2,"T_PL_CATS2/s");
- RootOutput::getInstance()->GetTree()->Branch("T_PL_CHIO",&vT_PL_CHIO,"T_PL_CHIO/s");
- RootOutput::getInstance()->GetTree()->Branch("T_PLchCATS1",&vT_PLchCATS1,"T_PLchCATS1/s");
- */
+/* RootOutput::getInstance()->GetTree()->Branch("ADC_CHIO_V15",&vADC_CHIO_V15,"ADC_CHIO_V15/s");
+ RootOutput::getInstance()->GetTree()->Branch("ADC_VOIE_29",&vADC_VOIE_29,"ADC_VOIE_29/s");
+ RootOutput::getInstance()->GetTree()->Branch("CHIO",&vCHIO,"CHIO/s");
+ RootOutput::getInstance()->GetTree()->Branch("CONFDEC",&vCONFDEC,"CONFDEC/s");
+ RootOutput::getInstance()->GetTree()->Branch("CONFDEC_AGAVA",&vCONFDEC_AGAVA,"CONFDEC_AGAVA/s");
+ RootOutput::getInstance()->GetTree()->Branch("DATATRIG",&vDATATRIG,"DATATRIG/s");
+ RootOutput::getInstance()->GetTree()->Branch("DATATRIG_CHIO",&vDATATRIG_CHIO,"DATATRIG_CHIO/s");
+ RootOutput::getInstance()->GetTree()->Branch("E1D6",&vE1D6,"E1D6/s");
+ RootOutput::getInstance()->GetTree()->Branch("E2D6",&vE2D6,"E2D6/s");
+ RootOutput::getInstance()->GetTree()->Branch("ED4",&vED4,"ED4/s");
+ RootOutput::getInstance()->GetTree()->Branch("EXL_HF",&vEXL_HF,"EXL_HF/s");
+ RootOutput::getInstance()->GetTree()->Branch("GALD4X",&vGALD4X,"GALD4X/s");
+ RootOutput::getInstance()->GetTree()->Branch("GALD4Y",&vGALD4Y,"GALD4Y/s");
+ RootOutput::getInstance()->GetTree()->Branch("QCaviar",&vQCaviar,"QCaviar/s");
+ RootOutput::getInstance()->GetTree()->Branch("QPlast",&vQPlast,"QPlast/s");
+ RootOutput::getInstance()->GetTree()->Branch("TCAVHF",&vTCAVHF,"TCAVHF/s");
+ RootOutput::getInstance()->GetTree()->Branch("TE1D6CAV",&vTE1D6CAV,"TE1D6CAV/s");
+ RootOutput::getInstance()->GetTree()->Branch("TE1D6GAL",&vTE1D6GAL,"TE1D6GAL/s");
+ RootOutput::getInstance()->GetTree()->Branch("TE1D6HF",&vTE1D6HF,"TE1D6HF/s");
+ RootOutput::getInstance()->GetTree()->Branch("TED4HF",&vTED4HF,"TED4HF/s");
+ RootOutput::getInstance()->GetTree()->Branch("TGALD4HF",&vTGALD4HF,"TGALD4HF/s");
+ RootOutput::getInstance()->GetTree()->Branch("T_CATS1_2",&vT_CATS1_2,"T_CATS1_2/s");
+ RootOutput::getInstance()->GetTree()->Branch("T_CATS1_CAV",&vT_CATS1_CAV,"T_CATS1_CAV/s");
+ RootOutput::getInstance()->GetTree()->Branch("T_MUVI_CATS1",&vT_MUVI_CATS1,"T_MUVI_CATS1/s");
+ RootOutput::getInstance()->GetTree()->Branch("T_PL_CATS1",&vT_PL_CATS1,"T_PL_CATS1/s");
+ RootOutput::getInstance()->GetTree()->Branch("T_PL_CATS2",&vT_PL_CATS2,"T_PL_CATS2/s");
+ RootOutput::getInstance()->GetTree()->Branch("T_PL_CHIO",&vT_PL_CHIO,"T_PL_CHIO/s");
+ RootOutput::getInstance()->GetTree()->Branch("T_PLchCATS1",&vT_PLchCATS1,"T_PLchCATS1/s");
+*/
}
////////////////////////////////////////////////////////////////////////////////
@@ -276,37 +276,37 @@ void Analysis::InitInputBranch(){
/* RootInput:: getInstance()->GetChain()->SetBranchStatus("InitialConditions",true ); */
/* RootInput:: getInstance()->GetChain()->SetBranchStatus("fIC_*",true ); */
/* RootInput:: getInstance()->GetChain()->SetBranchStatus("InitialConditions",true ); */
-
+
RootInput:: getInstance()->GetChain()->SetBranchAddress("GATCONF",&vGATCONF);
- /* RootInput:: getInstance()->GetChain()->SetBranchAddress("ADC_CHIO_V15",&vADC_CHIO_V15);
- RootInput:: getInstance()->GetChain()->SetBranchAddress("ADC_VOIE_29",&vADC_VOIE_29);
- RootInput:: getInstance()->GetChain()->SetBranchAddress("CHIO",&vCHIO);
- RootInput:: getInstance()->GetChain()->SetBranchAddress("CONFDEC",&vCONFDEC);
- RootInput:: getInstance()->GetChain()->SetBranchAddress("CONFDEC_AGAVA",&vCONFDEC_AGAVA);
- RootInput:: getInstance()->GetChain()->SetBranchAddress("DATATRIG",&vDATATRIG);
- RootInput:: getInstance()->GetChain()->SetBranchAddress("DATATRIG_CHIO",&vDATATRIG_CHIO);
- RootInput:: getInstance()->GetChain()->SetBranchAddress("E1D6",&vE1D6);
- RootInput:: getInstance()->GetChain()->SetBranchAddress("E2D6",&vE2D6);
- RootInput:: getInstance()->GetChain()->SetBranchAddress("ED4",&vED4);
- RootInput:: getInstance()->GetChain()->SetBranchAddress("EXL_HF",&vEXL_HF);
- RootInput:: getInstance()->GetChain()->SetBranchAddress("GALD4X",&vGALD4X);
- RootInput:: getInstance()->GetChain()->SetBranchAddress("GALD4Y",&vGALD4Y);
- RootInput:: getInstance()->GetChain()->SetBranchAddress("QCaviar",&vQCaviar);
- RootInput:: getInstance()->GetChain()->SetBranchAddress("QPlast",&vQPlast);
- RootInput:: getInstance()->GetChain()->SetBranchAddress("TCAVHF",&vTCAVHF);
- RootInput:: getInstance()->GetChain()->SetBranchAddress("TE1D6CAV",&vTE1D6CAV);
- RootInput:: getInstance()->GetChain()->SetBranchAddress("TE1D6GAL",&vTE1D6GAL);
- RootInput:: getInstance()->GetChain()->SetBranchAddress("TE1D6HF",&vTE1D6HF);
- RootInput:: getInstance()->GetChain()->SetBranchAddress("TED4HF",&vTED4HF);
- RootInput:: getInstance()->GetChain()->SetBranchAddress("TGALD4HF",&vTGALD4HF);
- RootInput:: getInstance()->GetChain()->SetBranchAddress("T_CATS1_2",&vT_CATS1_2);
- RootInput:: getInstance()->GetChain()->SetBranchAddress("T_CATS1_CAV",&vT_CATS1_CAV);
- RootInput:: getInstance()->GetChain()->SetBranchAddress("T_MUVI_CATS1",&vT_MUVI_CATS1);
- RootInput:: getInstance()->GetChain()->SetBranchAddress("T_PL_CATS1",&vT_PL_CATS1);
- RootInput:: getInstance()->GetChain()->SetBranchAddress("T_PL_CATS2",&vT_PL_CATS2);
- RootInput:: getInstance()->GetChain()->SetBranchAddress("T_PL_CHIO",&vT_PL_CHIO);
- RootInput:: getInstance()->GetChain()->SetBranchAddress("T_PLchCATS1",&vT_PLchCATS1);
- */
+ /* RootInput:: getInstance()->GetChain()->SetBranchAddress("ADC_CHIO_V15",&vADC_CHIO_V15);
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("ADC_VOIE_29",&vADC_VOIE_29);
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("CHIO",&vCHIO);
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("CONFDEC",&vCONFDEC);
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("CONFDEC_AGAVA",&vCONFDEC_AGAVA);
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("DATATRIG",&vDATATRIG);
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("DATATRIG_CHIO",&vDATATRIG_CHIO);
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("E1D6",&vE1D6);
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("E2D6",&vE2D6);
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("ED4",&vED4);
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("EXL_HF",&vEXL_HF);
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("GALD4X",&vGALD4X);
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("GALD4Y",&vGALD4Y);
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("QCaviar",&vQCaviar);
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("QPlast",&vQPlast);
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("TCAVHF",&vTCAVHF);
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("TE1D6CAV",&vTE1D6CAV);
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("TE1D6GAL",&vTE1D6GAL);
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("TE1D6HF",&vTE1D6HF);
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("TED4HF",&vTED4HF);
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("TGALD4HF",&vTGALD4HF);
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("T_CATS1_2",&vT_CATS1_2);
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("T_CATS1_CAV",&vT_CATS1_CAV);
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("T_MUVI_CATS1",&vT_MUVI_CATS1);
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("T_PL_CATS1",&vT_PL_CATS1);
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("T_PL_CATS2",&vT_PL_CATS2);
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("T_PL_CHIO",&vT_PL_CHIO);
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("T_PLchCATS1",&vT_PLchCATS1);
+*/
}
////////////////////////////////////////////////////////////////////////////////
@@ -336,13 +336,13 @@ NPL::VAnalysis* Analysis::Construct(){
// Registering the construct method to the factory //
////////////////////////////////////////////////////////////////////////////////
extern "C"{
- class proxy{
- public:
- proxy(){
- NPL::AnalysisFactory::getInstance()->SetConstructor(Analysis::Construct);
- }
- };
+class proxy{
+ public:
+ proxy(){
+ NPL::AnalysisFactory::getInstance()->SetConstructor(Analysis::Construct);
+ }
+};
- proxy p;
+proxy p;
}
diff --git a/Projects/e748/Analysis.h b/Projects/e748/Analysis.h
index 4566a7d1b..6bac47cdc 100644
--- a/Projects/e748/Analysis.h
+++ b/Projects/e748/Analysis.h
@@ -74,10 +74,10 @@ class Analysis: public NPL::VAnalysis{
double BeamEnergy;
int run_major;
int run_minor;
- NPL::EnergyLoss LightTarget;
- NPL::EnergyLoss LightAl ;
- NPL::EnergyLoss LightSi ;
- NPL::EnergyLoss BeamTarget ;
+ NPL::EnergyLoss He3CD2 ;
+ NPL::EnergyLoss He3Al ;
+ NPL::EnergyLoss He3Si ;
+ NPL::EnergyLoss BeamCD2 ;
NPL::EnergyLoss BeamMylar ;
NPL::EnergyLoss BeamIsobutane ;
diff --git a/Projects/e748/configs/ConfigMust2.dat b/Projects/e748/configs/ConfigMust2.dat
old mode 100644
new mode 100755
diff --git a/Projects/e748/e748.detector b/Projects/e748/e748.detector
old mode 100644
new mode 100755
index a1821cf6b..90336a06a
--- a/Projects/e748/e748.detector
+++ b/Projects/e748/e748.detector
@@ -116,7 +116,7 @@ CATSDetector
ModularLeaf
DefaultValue= -1000
- Leafs= QPlast QCaviar CHIO T_CATS1_CAV T_CATS1_2 T_MUVI_CATS1 T_PL_CATS2 EXL_HF T_PL_CHIO T_PL_CATS E1D6
+ Leafs= QPlast QCaviar CHIO T_CATS1_CAV T_CATS1_2 T_MUVI_CATS1 T_PL_CATS2 EXL_HF T_PL_CHIO T_PL_CATS1
Chio_dig
Pos= -30. -30. 600. mm
--
GitLab