Skip to content
Snippets Groups Projects
Commit 6fe017a4 authored by Morfouace's avatar Morfouace
Browse files

Merge branch 'NPTool.2.dev' of https://github.com/adrien-matta/nptool into NPTool.2.dev

parents 55cdcefe f1eb934f
No related branches found
No related tags found
No related merge requests found
......@@ -43,13 +43,13 @@ void Analysis::Init(){
LightCD2 = EnergyLoss("deuteron_CD2.G4table","G4Table",10);
LightAl = EnergyLoss("deuteron_Al.G4table","G4Table",10);
// LightAl = EnergyLoss("alpha_Al.G4table","G4Table",10);
// LightAl = EnergyLoss("alpha_Al.G4table","G4Table",10);
BeamCD2 = EnergyLoss("Si28_CD2.G4table","G4Table",10);
myReaction = new NPL::Reaction();
myReaction->ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
TargetThickness = m_DetectorManager->GetTargetThickness()*micrometer;
// TargetThickness = 0;
// TargetThickness = 0;
OriginalBeamEnergy = myReaction->GetBeamEnergy();
Rand = TRandom3();
DetectorNumber = 0 ;
......@@ -72,31 +72,43 @@ void Analysis::Init(){
myReaction->SetBeamEnergy(finalEnergy);
cout << "Set Beam energy to: " << finalEnergy << " MeV" << endl;
// Load cut for PAD cal
TFile* file ;
file = new TFile("cuts/deuton_d10.root","READ");
cut_deuton_d10 = (TCutG*) file->FindObjectAny("deuton_d10");
file = new TFile("cuts/proton_d10.root","READ");
cut_proton_d10 = (TCutG*) file->FindObjectAny("proton_d10");
file = new TFile("cuts/deuton_d12.root","READ");
cut_deuton_d12 = (TCutG*) file->FindObjectAny("deuton_d12");
file = new TFile("cuts/proton_d12.root","READ");
cut_proton_d12 = (TCutG*) file->FindObjectAny("proton_d12");
// Load cut for second Pad cal
/* TFile* file ;
file = new TFile("cuts/cut_elgs_det10.root","READ");
cut_deuton_d10 = (TCutG*) file->FindObjectAny("cut_elgs_det10");
file = new TFile("cuts/cut_elgs_det12.root","READ");
cut_deuton_d12 = (TCutG*) file->FindObjectAny("cut_elgs_det12");
file = new TFile("cuts/Cut_dd.root","READ");
cut_dd = (TCutG*) file->FindObjectAny("Cut_dd");
PADFile.open("PADEvent.txt");
*/
/* // Load the cut for alignement
for(unsigned int i = 0 ; i < 8 ; i++){
TFile* file = new TFile(Form("cuts/Ex5_D%i.root",i+1),"READ");
cut_ex5.push_back((TCutG*) file->FindObjectAny(Form("Ex5_D%i",i+1)));
if(! file->FindObjectAny(Form("Ex5_D%i",i+1))){
// Load cut for PAD cal
/* TFile* file ;
file = new TFile("cuts/deuton_d10.root","READ");
cut_deuton_d10 = (TCutG*) file->FindObjectAny("deuton_d10");
file = new TFile("cuts/proton_d10.root","READ");
cut_proton_d10 = (TCutG*) file->FindObjectAny("proton_d10");
file = new TFile("cuts/deuton_d12.root","READ");
cut_deuton_d12 = (TCutG*) file->FindObjectAny("deuton_d12");
file = new TFile("cuts/proton_d12.root","READ");
cut_proton_d12 = (TCutG*) file->FindObjectAny("proton_d12");
PADFile.open("PADEvent.txt");
*/
/* // Load the cut for alignement
for(unsigned int i = 0 ; i < 8 ; i++){
TFile* file = new TFile(Form("cuts/Ex5_D%i.root",i+1),"READ");
cut_ex5.push_back((TCutG*) file->FindObjectAny(Form("Ex5_D%i",i+1)));
if(! file->FindObjectAny(Form("Ex5_D%i",i+1))){
cout << "Fail to Load " << Form("Ex5_D%i",i+1) << endl;
exit(1);
}
}
}
}
box_pos.open("BoxPos.txt");
qqq_pos.open("QQQPos.txt");
*/
box_pos.open("BoxPos.txt");
qqq_pos.open("QQQPos.txt");
*/
}
////////////////////////////////////////////////////////////////////////////////
......@@ -113,7 +125,8 @@ void Analysis::TreatEvent(){
TVector3 HitDirection = Sharc -> GetPositionOfInteraction(0)-TargetPosition;
ThetaLab = HitDirection.Angle( BeamDirection );
ThetaNormalTarget = HitDirection.Angle( TVector3(0,0,1) ) ;
ThetaDetector = HitDirection.Angle(Sharc->GetDetectorNormal(0));
ThetaDetector = HitDirection.Angle(-Sharc->GetDetectorNormal(0));
Distance = HitDirection.Mag();
/************************************************/
/************************************************/
......@@ -129,46 +142,72 @@ void Analysis::TreatEvent(){
// Target Correction
ELab = LightCD2.EvaluateInitialEnergy( Energy ,TargetThickness*0.5, ThetaNormalTarget);
/************************************************/
int DetectorNumber = Sharc->DetectorNumber[0];
if(DetectorNumber==10 && Sharc->PAD_E[0]>0){
if(cut_proton_d10->IsInside(Sharc->PAD_E[0],-Sharc->Strip_E[0]*cos(ThetaDetector)))
PADFile << "10 proton " << Sharc->Strip_E[0] << " " << Sharc->PAD_E[0] << " " << ThetaDetector << endl;
else if(cut_deuton_d10->IsInside(Sharc->PAD_E[0],-Sharc->Strip_E[0]*cos(ThetaDetector)))
PADFile << "10 deuton " << Sharc->Strip_E[0] << " " << Sharc->PAD_E[0] << " " << ThetaDetector << endl;
}
/* int DetectorNumber = Sharc->DetectorNumber[0];
if(DetectorNumber==10 && Sharc->PAD_E[0]>0){
if(cut_proton_d10->IsInside(Sharc->PAD_E[0],Sharc->Strip_E[0]*cos(ThetaDetector))){
PADFile << "10 proton " << Sharc->Strip_E[0] << " " << Sharc->PAD_E[0] << " " << ThetaDetector << endl;
}
else if(cut_deuton_d10->IsInside(Sharc->PAD_E[0],Sharc->Strip_E[0]*cos(ThetaDetector)))
PADFile << "10 deuton " << Sharc->Strip_E[0] << " " << Sharc->PAD_E[0] << " " << ThetaDetector << endl;
}
if(DetectorNumber==12 && Sharc->PAD_E[0]>0){
if(cut_proton_d12->IsInside(Sharc->PAD_E[0],Sharc->Strip_E[0]*cos(ThetaDetector)))
PADFile << "12 proton " << Sharc->Strip_E[0] << " " << Sharc->PAD_E[0] << " " << ThetaDetector << endl;
else if(cut_deuton_d12->IsInside(Sharc->PAD_E[0], Sharc->Strip_E[0]*cos(ThetaDetector)))
PADFile << "12 deuton " << Sharc->Strip_E[0] << " " << Sharc->PAD_E[0] << " " << ThetaDetector << endl;
}
*/
if(DetectorNumber==12 && Sharc->PAD_E[0]>0){
if(cut_proton_d12->IsInside(Sharc->PAD_E[0],-Sharc->Strip_E[0]*cos(ThetaDetector)))
PADFile << "12 proton " << Sharc->Strip_E[0] << " " << Sharc->PAD_E[0] << " " << ThetaDetector << endl;
else if(cut_deuton_d12->IsInside(Sharc->PAD_E[0], -Sharc->Strip_E[0]*cos(ThetaDetector)))
PADFile << "12 deuton " << Sharc->Strip_E[0] << " " << Sharc->PAD_E[0] << " " << ThetaDetector << endl;
}
/************************************************/
/* int DetectorNumber = Sharc->DetectorNumber[0];
if(cut_dd->IsInside(Sharc->PAD_E[0],Sharc->Strip_E[0]*cos(ThetaDetector))){
if(DetectorNumber == 10){
if(cut_deuton_d10->IsInside(ThetaLab/deg,ELab)){
double a = 9.23351e-13 ;
double b = 8.35809e-06 ;
double c = 5.77842e-01-Sharc->PAD_E[0];
double delta = b*b-4*a*c;
double RawE = (-b + sqrt(delta))/(2*a);
PADFile << "10 deuton " << Sharc->Strip_E[0] << " " << RawE << " " << ThetaLab << " " << ThetaNormalTarget << endl;
}
}
else if(DetectorNumber == 12){
if(cut_deuton_d12->IsInside(ThetaLab/deg,ELab)){
double a = 1.63514e-12;
double b = 7.52050e-06;
double c = 5.93687e-01-Sharc->PAD_E[0];
double delta = b*b-4*a*c;
double RawE = (-b + sqrt(delta))/(2*a);
PADFile << "12 deuton " << Sharc->Strip_E[0] << " " << RawE<< " " << ThetaLab << " " << ThetaNormalTarget << endl;
}
}
/************************************************/
/************************************************/
/* int DetectorNumber = Sharc->DetectorNumber[0];
bool checkT = false;
if(Tigress->AddBack_DC.size()>0){
if(Tigress->AddBack_DC[0]/1000.> 4.6 && Tigress->AddBack_DC[0]/1000.< 5.4)
checkT=true;
}
}*/
if(checkT){
if(DetectorNumber<5){
if(cut_ex5[DetectorNumber-1]->IsInside(ThetaLab/deg,ELab))
/************************************************/
/* int DetectorNumber = Sharc->DetectorNumber[0];
bool checkT = false;
if(Tigress->AddBack_DC.size()>0){
if(Tigress->AddBack_DC[0]/1000.> 4.6 && Tigress->AddBack_DC[0]/1000.< 5.4)
checkT=true;
}
if(checkT){
if(DetectorNumber<5){
if(cut_ex5[DetectorNumber-1]->IsInside(ThetaLab/deg,ELab))
qqq_pos << DetectorNumber << " " << Energy << " " << HitDirection.X() << " " << HitDirection.Y() << " " << HitDirection.Z() << endl;
}
}
else if(DetectorNumber<9){
if(cut_ex5[DetectorNumber-1]->IsInside(ThetaLab/deg,ELab))
else if(DetectorNumber<9){
if(cut_ex5[DetectorNumber-1]->IsInside(ThetaLab/deg,ELab))
box_pos << DetectorNumber << " " << Energy << " " << HitDirection.X() << " " << HitDirection.Y() << " " << HitDirection.Z() << endl;
}
}
*/
}
}
*/
/************************************************/
/************************************************/
......@@ -196,6 +235,7 @@ void Analysis::InitOutputBranch(){
RootOutput::getInstance()->GetTree()->Branch("ThetaLab",&ThetaLab,"ThetaLab/D");
RootOutput::getInstance()->GetTree()->Branch("ThetaCM",&ThetaCM,"ThetaCM/D");
RootOutput::getInstance()->GetTree()->Branch("ThetaDetector",&ThetaDetector,"ThetaDetector/D");
RootOutput::getInstance()->GetTree()->Branch("Distance",&Distance,"Distance/D");
RootOutput::getInstance()->GetTree()->Branch("RunNumber",&RunNumber,"RunNumber/I");
RootOutput::getInstance()->GetTree()->Branch("RunNumberMinor",&RunNumberMinor,"RunNumberMinor/I");
......
......@@ -85,13 +85,14 @@ TInitialConditions* myInit ;
double ThetaDetector;
double Si_E_Sharc ;
double E_Sharc ;
double Distance;
TSharcPhysics* Sharc;
TTigressPhysics* Tigress;
std::vector<TCutG*> cut_ex5;
std::ofstream box_pos;
std::ofstream qqq_pos;
TCutG* cut_dd;
TCutG* cut_deuton_d10;
TCutG* cut_proton_d10;
TCutG* cut_deuton_d12;
......
......@@ -48,6 +48,8 @@ void TTrifoilData::Clear(){
fTrifoil_Waveform.clear();
fTrifoil_TimeCFD.clear();
fTrifoil_TimeLED.clear();
fTrifoil_TimeStamp.clear();
}
/////////////////////////
void TTrifoilData::Dump() const{
......
......@@ -39,7 +39,8 @@ private:
vector<TH1F> fTrifoil_Waveform;
vector<Double_t> fTrifoil_TimeCFD;
vector<Double_t> fTrifoil_TimeLED;
vector<Double_t> fTrifoil_TimeStamp;
public:
TTrifoilData();
virtual ~TTrifoilData();
......@@ -52,11 +53,13 @@ public:
inline void SetWaveform(const TH1F &Waveform) {fTrifoil_Waveform.push_back(Waveform);}
inline void SetTimeCFD(const Double_t &TimeCFD) {fTrifoil_TimeCFD.push_back(TimeCFD);}
inline void SetTimeLED(const Double_t &TimeLED) {fTrifoil_TimeLED.push_back(TimeLED);}
inline void SetTimeStamp(const Double_t &TimeStamp) {fTrifoil_TimeStamp.push_back(TimeStamp);}
///////////////////// GETTERS ////////////////////////
inline TH1F* GetWaveform(const unsigned int &i) {return &(fTrifoil_Waveform[i]);}
inline Double_t GetTimeCFD(const unsigned int &i) {return fTrifoil_TimeCFD[i];}
inline Double_t GetTimeLED(const unsigned int &i) {return fTrifoil_TimeLED[i];}
inline Double_t GetTimeStamp(const unsigned int &i) {return fTrifoil_TimeStamp[i];}
inline unsigned int GetMultiplicity() {return fTrifoil_Waveform.size();}
ClassDef(TTrifoilData,1) // TrifoilData structure
......
......@@ -52,13 +52,15 @@ void TTrifoilPhysics::BuildSimplePhysicalEvent(){
///////////////////////////////////////////////////////////////////////////
void TTrifoilPhysics::BuildPhysicalEvent(){
unsigned int mysize = m_EventData->GetMultiplicity();
double base,maxi,mytime;
for (unsigned int i = 0 ; i < mysize ; i++){
TH1F* h = m_EventData->GetWaveform(i);
double base = h->GetBinContent(h->GetMinimumBin());
double maxi = h->GetBinContent(h->GetMaximumBin());
base = h->GetBinContent(h->GetMinimumBin());
maxi = h->GetBinContent(h->GetMaximumBin());
if(maxi>2000 && base>-300){
Time.push_back(h->GetMaximumBin());
mytime = (h->GetMaximumBin());
TimeStamp.push_back(m_EventData->GetTimeStamp(i)*10);
Time.push_back(mytime);
Energy.push_back(maxi);
}
}
......@@ -68,6 +70,7 @@ void TTrifoilPhysics::BuildPhysicalEvent(){
void TTrifoilPhysics::Clear(){
Energy.clear() ;
Time.clear() ;
TimeStamp.clear();
}
///////////////////////////////////////////////////////////////////////////
......
......@@ -48,9 +48,10 @@ class TTrifoilPhysics : public TObject, public NPL::VDetector
public:
// EventType is True if the wave form analysis return a valid Trifoil event
vector<double> Energy;
vector<double> Energy;
vector<double> Time ;
vector<double> TimeStamp ;
public: // Innherited from VDetector Class
// Read stream at ConfigFile to pick-up parameters of detector (Position,...) using Token
......
......@@ -75,7 +75,13 @@ void EventGeneratorBeam::GenerateEvent(G4Event* anEvent){
if( anEvent->GetEventID()==0){
// Define the particle to be shoot
m_particle = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(m_Beam->GetZ(), m_Beam->GetA() ,m_Beam->GetExcitationEnergy());
if(m_Beam->GetZ()==0 && m_Beam->GetA()==1){
m_particle = G4ParticleTable::GetParticleTable()->FindParticle("neutron");
}
else
m_particle = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(m_Beam->GetZ(), m_Beam->GetA() ,m_Beam->GetExcitationEnergy());
}
///////////////////////////////////////////////////////////////////////
......
......@@ -243,6 +243,13 @@ G4Material* MaterialManager::GetMaterialFromLibrary(string Name){
m_Material[Name]=material;
return material;
}
else if(Name == "NaturalUranium"){
G4Material* material = new G4Material(Name, 19.1*g/cm3,1);
material->AddElement(GetElementFromLibrary("U"),1);
m_Material[Name]=material;
return material;
}
else if(Name == "CsI_Scintillator"){
G4Material* material = new G4Material(Name, 4.51*g/cm3,2);
......
......@@ -44,12 +44,15 @@
#include "G4eIonisation.hh"
#include "G4eBremsstrahlung.hh"
#include "G4eplusAnnihilation.hh"
#include "G4NeutronHPFission.hh"
#include "G4MuIonisation.hh"
#include "G4MuMultipleScattering.hh"
#include "G4MuBremsstrahlung.hh"
#include "G4MuPairProduction.hh"
//neutron
#include "G4HadronFissionProcess.hh"
#include "G4hIonisation.hh"
#include "G4ionIonisation.hh"
#include "G4hMultipleScattering.hh"
......@@ -79,7 +82,9 @@
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
PhysicsList::PhysicsList(){
// ie: no secondaries
defaultCutValue = 1000 * pc;
//defaultCutValue = 1000 * pc;
defaultCutValue = 1*mm;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
PhysicsList::~PhysicsList(){
......@@ -91,22 +96,22 @@ void PhysicsList::ConstructParticle(){
// for all particles which you want to use.
// This ensures that objects of these particle types will be
// created in the program.
//Usefull to test geometry
G4Geantino::GeantinoDefinition();
//Usefull for Gamma
ConstructBosons();
//Usefull for betaDecay
ConstructLeptons();
//Needed by G4 (4.9.2) to run on mac os X ;-)
ConstructMesons();
//usefull for p and n
ConstructBaryons();
//Usefull of course :p
ConstructIons();
}
......@@ -160,28 +165,28 @@ void PhysicsList::ConstructIons(){
G4Deuteron::DeuteronDefinition() ;
G4Triton::TritonDefinition() ;
G4Alpha::AlphaDefinition() ;
G4IonConstructor iConstructor ;
iConstructor.ConstructParticle() ;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void PhysicsList::ConstructProcess(){
AddTransportation() ;
ConstructEM() ;
SetCuts() ;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void PhysicsList::ConstructEM(){
theParticleIterator->reset();
while ((*theParticleIterator)()) {
G4ParticleDefinition* particle = theParticleIterator->value() ;
G4ProcessManager* pmanager = particle->GetProcessManager() ;
G4String particleName = particle->GetParticleName() ;
if (particleName == "gamma") {
// gamma
//standard Geant4
......@@ -196,7 +201,7 @@ void PhysicsList::ConstructEM(){
pmanager->AddProcess(new G4eIonisation , -1, 2, 2) ;
pmanager->AddProcess(new G4eBremsstrahlung , -1, -1, 3) ;
}
else if (particleName == "e+") {
//positron
pmanager->AddProcess(new G4eMultipleScattering , -1, 1, 1 );
......@@ -205,54 +210,61 @@ void PhysicsList::ConstructEM(){
pmanager->AddProcess(new G4eBremsstrahlung , -1, -1, 3 );
pmanager->AddProcess(new G4eplusAnnihilation , 0, -1, 4 );
}
else if (particleName == "mu+" ||
particleName == "mu-") {
particleName == "mu-") {
}
else if ( particleName == "GenericIon"
|| particleName == "He3"
|| particleName == "alpha"
|| particleName == "deuteron"
|| particleName == "triton") {
|| particleName == "He3"
|| particleName == "alpha"
|| particleName == "deuteron"
|| particleName == "triton") {
pmanager->AddProcess(new G4hMultipleScattering(), -1, 1, 1) ;
G4ionIonisation* iI = new G4ionIonisation ;
iI->ActivateStoppingData(true) ;
pmanager->AddProcess(iI , -1, 2, 2) ;
}
else if (particleName == "neutron"){
G4NeutronHPFission* fissionModel = new G4NeutronHPFission();
G4HadronFissionProcess* fissionProcess = new G4HadronFissionProcess();
fissionProcess->RegisterMe(fissionModel);
pmanager->AddDiscreteProcess(fissionProcess);
}
else if ((!particle->IsShortLived()) &&
(particle->GetPDGCharge() != 0.0) &&
(particleName != "chargedgeantino")) {
(particle->GetPDGCharge() != 0.0) &&
(particleName != "chargedgeantino")) {
G4hIonisation* hI = new G4hIonisation ;
pmanager->AddProcess(new G4hMultipleScattering , -1, 1, 1) ;
pmanager->AddProcess(hI , -1, 2, 2) ;
}
}
G4EmProcessOptions opt ;
opt.SetSubCutoff(true) ;
opt.SetMinEnergy(0.001*eV) ;
opt.SetMaxEnergy(50.*GeV) ;
opt.SetDEDXBinning(5000) ;
opt.SetLambdaBinning(5000) ;
//energy loss
opt.SetLinearLossLimit(1.e-9);
opt.SetStepFunction(0.001, 10.*um);
//Multiple scattering
opt.SetMscLateralDisplacement(true);
opt.SetLossFluctuations(true);
//opt.SetMscStepLimitation(fMinimal);
//opt.SetMscStepLimitation(fUseDistanceToBoundary);
opt.SetMscStepLimitation(fUseSafety);
//ionization
opt.SetSubCutoff(false);
//energy loss
opt.SetLinearLossLimit(1.e-9);
opt.SetStepFunction(0.001, 10.*um);
//Multiple scattering
opt.SetMscLateralDisplacement(true);
opt.SetLossFluctuations(true);
//opt.SetMscStepLimitation(fMinimal);
//opt.SetMscStepLimitation(fUseDistanceToBoundary);
opt.SetMscStepLimitation(fUseSafety);
//ionization
opt.SetSubCutoff(false);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
......@@ -263,20 +275,20 @@ void PhysicsList::SetCuts(){
// " G4VUserPhysicsList::SetCutsWithDefault" method sets
// the default cut value for all particle types
SetCutsWithDefault();
// for gamma-rays
/* G4double em_cuts = 0.01*mm; // Use 10cm to avoid generating delta-rays
//G4double em_cuts = 10.*cm; // Use 10cm to avoid generating delta-rays
SetCutValue(em_cuts, "gamma");
SetCutValue(em_cuts, "e-");
SetCutValue(em_cuts, "e+");*/
//G4double em_cuts = 10.*cm; // Use 10cm to avoid generating delta-rays
SetCutValue(em_cuts, "gamma");
SetCutValue(em_cuts, "e-");
SetCutValue(em_cuts, "e+");*/
// Retrieve verbose level
SetVerboseLevel(temp);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void PhysicsList::ConstructDecay(){
// Add Decay Process
G4Decay* theDecayProcess = new G4Decay() ;
theParticleIterator->reset() ;
......@@ -291,7 +303,7 @@ void PhysicsList::ConstructDecay(){
}
}
//end Add Decay Process
}
void PhysicsList::MyOwnConstruction(){
......
......@@ -501,6 +501,10 @@ G4double Target::SlowDownBeam(G4ParticleDefinition* Beam,
G4double ZInteraction,
G4double IncidentTheta){
if(Beam->GetParticleName()=="neutron"){
return IncidentEnergy;
}
G4double ThicknessBeforeInteraction =
abs(ZInteraction - 0.5*m_EffectiveThickness) / cos(m_TargetAngle);
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment