From 8e20575c309c122a7f233ad3394fa186bb5d2e29 Mon Sep 17 00:00:00 2001
From: Adrien Matta <matta@lpccaen.in2p3.fr>
Date: Fri, 26 Jun 2020 10:15:33 +0200
Subject: [PATCH] * fixing small bugs in vertex determination of simulation
---
NPLib/Core/NPOptionManager.cxx | 3 +
NPLib/Core/NPOptionManager.h | 7 ++-
NPSimulation/Process/BeamReaction.cc | 42 ++++---------
NPSimulation/Process/BeamReaction.hh | 1 +
Projects/MUGAST/Analysis.cxx | 78 ++++++++++++++++++------
Projects/MUGAST/Analysis.h | 9 ++-
Projects/MUGAST/ShowResults.C | 90 ++++++++++++----------------
7 files changed, 129 insertions(+), 101 deletions(-)
diff --git a/NPLib/Core/NPOptionManager.cxx b/NPLib/Core/NPOptionManager.cxx
index c4e1573dc..8a42e82e6 100644
--- a/NPLib/Core/NPOptionManager.cxx
+++ b/NPLib/Core/NPOptionManager.cxx
@@ -168,6 +168,8 @@ void NPOptionManager::ReadTheInputArgument(int argc, char** argv){
else if (argument == "--circular") {fCircularTree = true;}
+ else if (argument == "--definition" && argc >= i + 1) {std::string def= argv[++i];fDefinition.insert(def);}
+
else{
SendErrorAndExit(argument.c_str());
}
@@ -401,6 +403,7 @@ void NPOptionManager::DisplayHelp(){
std::cout << "\t--event-generator　-E <arg>\tSet arg as the event generator file" << std::endl ;
std::cout << "\t--output　-O <arg>\t\tSet arg as the Output File Name (output tree)" << std::endl ;
std::cout << "\t--tree-name <arg>\t\tSet arg as the Output Tree Name " << std::endl ;
+ std::cout << "\t--definition <definition> \tAdd <definition> to the list of definition" << std::endl ;
std::cout << "\t--verbose -V <arg>\t\tSet the verbose level, 0 for nothing, 1 for normal printout."<<std::endl;
std::cout << "\t\t\t\t\tError and warning are not affected" << std::endl ;
std::cout << std::endl << "NPAnalysis only:"<<std::endl;
diff --git a/NPLib/Core/NPOptionManager.h b/NPLib/Core/NPOptionManager.h
index fda9c185b..a628d8fe0 100644
--- a/NPLib/Core/NPOptionManager.h
+++ b/NPLib/Core/NPOptionManager.h
@@ -26,6 +26,7 @@
// C++ headers
#include <iostream>
#include <string>
+#include <set>
class NPOptionManager{
public:
@@ -116,7 +117,10 @@ class NPOptionManager{
void SetDetectorFile(const std::string& name) {fDetectorFileName = name;CheckDetectorConfiguration();}
void SetRunToReadFile(const std::string& name) {fRunToReadFileName = name;}
void SetVerboseLevel(int VerboseLevel) {fVerboseLevel = VerboseLevel;}
-
+
+ public: // user definition
+ bool HasDefinition(std::string def) {return(fDefinition.find(def)!=fDefinition.end());}
+
private:
// default values
std::string fDefaultReactionFileName;
@@ -153,6 +157,7 @@ class NPOptionManager{
std::string fSharedLibExtension; // lib extension is platform dependent
std::string fG4MacroPath; // Path to a geant4 macro to execute at start of nps
bool fG4BatchMode; // Execute geant4 in batch mode, running the given macro
+ std::set<std::string> fDefinition; // a set of user defined definition
};
#endif
diff --git a/NPSimulation/Process/BeamReaction.cc b/NPSimulation/Process/BeamReaction.cc
index caa7369ab..bd5ca564c 100644
--- a/NPSimulation/Process/BeamReaction.cc
+++ b/NPSimulation/Process/BeamReaction.cc
@@ -39,6 +39,7 @@ NPS::BeamReaction::BeamReaction(G4String modelName, G4Region* envelope)
m_PreviousEnergy = 0;
m_PreviousLength = 0;
m_PreviousDirection.set(0,0,0);
+ m_shoot=false;
}
////////////////////////////////////////////////////////////////////////////////
@@ -98,8 +99,7 @@ void NPS::BeamReaction::ReadConfiguration() {
}
////////////////////////////////////////////////////////////////////////////////
-G4bool
-NPS::BeamReaction::IsApplicable(const G4ParticleDefinition& particleType) {
+G4bool NPS::BeamReaction::IsApplicable(const G4ParticleDefinition& particleType) {
if (!m_active)
return false;
@@ -115,7 +115,6 @@ NPS::BeamReaction::IsApplicable(const G4ParticleDefinition& particleType) {
G4bool NPS::BeamReaction::ModelTrigger(const G4FastTrack& fastTrack) {
//cout<< "--------- MODEL TRIG ---------"<<endl;
- static bool shoot = false;
static double rand = 0;
const G4Track* PrimaryTrack = fastTrack.GetPrimaryTrack();
G4ThreeVector V = PrimaryTrack->GetMomentum().unit();
@@ -128,56 +127,40 @@ G4bool NPS::BeamReaction::ModelTrigger(const G4FastTrack& fastTrack) {
double out = solid->DistanceToOut(P, -V);
double ratio = in / (out + in);
- //cout.precision(10);
- //cout<<"in:"<<in<<endl;
- //cout<<"out:"<<out<<endl;
- //cout<<"ratio:"<<ratio<<endl;
-
m_Parent_ID = PrimaryTrack->GetParentID();
// process reserved to the beam
if(m_Parent_ID!=0)
return false;
- if (out == 0) { // first step of current event
- //cout<< "FIRST/////////////////"<<endl;
+ if(!m_shoot){
rand = G4RandFlat::shoot();
m_PreviousLength = m_StepSize;
m_PreviousEnergy = PrimaryTrack->GetKineticEnergy();
m_PreviousDirection = PrimaryTrack->GetMomentumDirection();
+
// Clear Previous Event
m_ReactionConditions->Clear();
- shoot = true;
+ m_shoot=true;
}
- else if (((in-m_StepSize) <= 1E-9) && shoot) { // last step
- //cout<< "LAST"<<endl;
+ else if (((in-m_StepSize) <= 1E-20) && m_shoot) { // last step
return true;
}
- //cout<< "rand:"<<rand<<std::scientific<<endl;
// If the condition is met, the event is generated
if (ratio < rand) {
// Reset the static for next event
- // shoot = false;
if(m_ReactionType=="QFSReaction"){
- if ( shoot && m_QFS.IsAllowed() ) {
- shoot = false;
+ if ( m_shoot && m_QFS.IsAllowed() ) {
return true;
- } else {
- return false;
- }
+ }
}
else if(m_ReactionType=="TwoBodyReaction"){
- if ( shoot && m_Reaction.IsAllowed(PrimaryTrack->GetKineticEnergy()) ) {
- shoot = false;
+ if ( m_shoot && m_Reaction.IsAllowed(PrimaryTrack->GetKineticEnergy()) ) {
return true;
- } else {
- return false;
- }
- }else{
- return false;
+ }
}
}
@@ -185,7 +168,6 @@ G4bool NPS::BeamReaction::ModelTrigger(const G4FastTrack& fastTrack) {
// so it can be used in the next one
if (!PrimaryTrack->GetStep()->IsLastStepInVolume()) {
m_PreviousLength = PrimaryTrack->GetStep()->GetStepLength();
- //cout<< "PreviousLength="<<m_PreviousLength<<endl;
m_PreviousEnergy = PrimaryTrack->GetKineticEnergy();
m_PreviousDirection = PrimaryTrack->GetMomentumDirection();
}
@@ -196,8 +178,8 @@ G4bool NPS::BeamReaction::ModelTrigger(const G4FastTrack& fastTrack) {
////////////////////////////////////////////////////////////////////////////////
void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
G4FastStep& fastStep) {
-
- //cout<< "--------- DO IT ---------"<<endl;
+
+ m_shoot=false;
// Get the track info
const G4Track* PrimaryTrack = fastTrack.GetPrimaryTrack();
diff --git a/NPSimulation/Process/BeamReaction.hh b/NPSimulation/Process/BeamReaction.hh
index 78ebc0fd8..7f1f519ee 100644
--- a/NPSimulation/Process/BeamReaction.hh
+++ b/NPSimulation/Process/BeamReaction.hh
@@ -52,6 +52,7 @@ namespace NPS{
double m_PreviousLength;
G4ThreeVector m_PreviousDirection;
bool m_active;// is the process active
+ bool m_shoot;
double m_StepSize;
int m_Parent_ID;
diff --git a/Projects/MUGAST/Analysis.cxx b/Projects/MUGAST/Analysis.cxx
index 30a945f89..4c86b827c 100644
--- a/Projects/MUGAST/Analysis.cxx
+++ b/Projects/MUGAST/Analysis.cxx
@@ -36,16 +36,29 @@ Analysis::~Analysis(){
////////////////////////////////////////////////////////////////////////////////
void Analysis::Init() {
-
+ if(NPOptionManager::getInstance()->HasDefinition("simulation")){
+ cout << "Considering input data as simulation" << endl;
+ simulation=true;
+ }
+ else{
+ cout << "Considering input data as real" << endl;
+ simulation=false;
+ }
agata_zShift=51*mm;
// initialize input and output branches
+ if(simulation){
+ Initial = new TInitialConditions();
+ ReactionConditions = new TReactionConditions();
+ }
+
InitOutputBranch();
InitInputBranch();
// get MUST2 and Gaspard objects
M2 = (TMust2Physics*) m_DetectorManager -> GetDetector("M2Telescope");
MG = (TMugastPhysics*) m_DetectorManager -> GetDetector("Mugast");
- CATS = (TCATSPhysics*) m_DetectorManager->GetDetector("CATSDetector");
+ if(!simulation)
+ CATS = (TCATSPhysics*) m_DetectorManager->GetDetector("CATSDetector");
// get reaction information
reaction.ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
@@ -65,7 +78,12 @@ void Analysis::Init() {
LightSi = NPL::EnergyLoss(light+"_Si.G4table","G4Table",100);
BeamCD2 = NPL::EnergyLoss(beam+"_"+TargetMaterial+".G4table","G4Table",100);
+ FinalBeamEnergy = BeamCD2.Slow(OriginalBeamEnergy,TargetThickness*0.5,0);
+ cout << "Original beam energy: " << OriginalBeamEnergy << " MeV Mid-target beam energy: " << FinalBeamEnergy << "MeV " << endl;
+ reaction.SetBeamEnergy(FinalBeamEnergy);
+
if(WindowsThickness){
+ cout << "Cryogenic target with windows" << endl;
BeamWindow= new NPL::EnergyLoss(beam+"_"+WindowsMaterial+".G4table","G4Table",100);
LightWindow= new NPL::EnergyLoss(light+"_"+WindowsMaterial+".G4table","G4Table",100);
}
@@ -97,10 +115,20 @@ void Analysis::Init() {
void Analysis::TreatEvent() {
// Reinitiate calculated variable
ReInitValue();
- double XTarget = CATS->GetPositionOnTarget().X();
- double YTarget = CATS->GetPositionOnTarget().Y();
- TVector3 BeamDirection = CATS->GetBeamDirection();
-
+ double XTarget, YTarget;
+ TVector3 BeamDirection;
+ if(!simulation){
+ XTarget = CATS->GetPositionOnTarget().X();
+ YTarget = CATS->GetPositionOnTarget().Y();
+ BeamDirection = CATS->GetBeamDirection();
+ }
+ else{
+ XTarget = 0;
+ YTarget = 0;
+ BeamDirection = TVector3(0,0,1);
+ OriginalELab = ReactionConditions->GetKineticEnergy(0);
+ OriginalThetaLab = ReactionConditions->GetTheta(0);
+ }
BeamImpact = TVector3(XTarget,YTarget,0);
// determine beam energy for a randomized interaction point in target
// 1% FWHM randominastion (E/100)/2.35
@@ -151,7 +179,7 @@ void Analysis::TreatEvent() {
// Evaluate energy using the thickness
ELab = LightAl.EvaluateInitialEnergy( Energy ,0.4*micrometer , ThetaM2Surface);
// Target Correction
- ELab = LightTarget.EvaluateInitialEnergy( ELab ,TargetThickness/2., ThetaNormalTarget);
+ ELab = LightTarget.EvaluateInitialEnergy( ELab ,TargetThickness*0.5, ThetaNormalTarget);
if(LightWindow)
ELab = LightWindow->EvaluateInitialEnergy( ELab ,WindowsThickness, ThetaNormalTarget);
@@ -192,8 +220,10 @@ void Analysis::TreatEvent() {
// Part 2 : Impact Energy
Energy = ELab = 0;
Energy = MG->GetEnergyDeposit(countMugast);
+
+ // ELab = LightAl.EvaluateInitialEnergy( Energy ,0.4*micrometer , ThetaMGSurface);
// Target Correction
- ELab = LightTarget.EvaluateInitialEnergy( Energy ,TargetThickness*0.5, ThetaNormalTarget);
+ ELab = LightTarget.EvaluateInitialEnergy( Energy,TargetThickness*0.5, ThetaNormalTarget);
if(LightWindow)
ELab = LightWindow->EvaluateInitialEnergy( ELab ,WindowsThickness, ThetaNormalTarget);
@@ -207,7 +237,7 @@ void Analysis::TreatEvent() {
ThetaLab=ThetaLab/deg;
}//end loop Mugast
-
+
////////////////////////////////////////////////////////////////////////////
///////////////////////////////// LOOP on AGATA ////////////////////////////
////////////////////////////////////////////////////////////////////////////
@@ -230,7 +260,7 @@ void Analysis::TreatEvent() {
GammaLV.Boost(-beta);
// Get EDC
EDC=GammaLV.Energy();
- }
+ }
}
@@ -271,6 +301,10 @@ void Analysis::InitOutputBranch() {
RootOutput::getInstance()->GetTree()->Branch("coreTS",coreTS,"coreTS[nbCores]/l");
RootOutput::getInstance()->GetTree()->Branch("coreE0",coreE0,"coreE0[nbCores]/F");
//
+ if(simulation){
+ RootOutput::getInstance()->GetTree()->Branch("OriginalELab",&OriginalELab,"OriginalELab/D");
+ RootOutput::getInstance()->GetTree()->Branch("OriginalThetaLab",&OriginalThetaLab,"OriginalThetaLab/D");
+ }
}
@@ -292,6 +326,14 @@ void Analysis::InitInputBranch(){
RootInput::getInstance()->GetChain()->SetBranchAddress("coreId",coreId);
RootInput::getInstance()->GetChain()->SetBranchAddress("coreTS",coreTS);
RootInput::getInstance()->GetChain()->SetBranchAddress("coreE0",coreE0);
+ if(simulation){
+ RootInput:: getInstance()->GetChain()->SetBranchStatus("InitialConditions",true );
+ RootInput:: getInstance()->GetChain()->SetBranchStatus("fIC_*",true );
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("InitialConditions",&Initial);
+ RootInput:: getInstance()->GetChain()->SetBranchStatus("ReactionConditions",true );
+ RootInput:: getInstance()->GetChain()->SetBranchStatus("fRC_*",true );
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("ReactionConditions",&ReactionConditions);
+ }
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::ReInitValue(){
@@ -318,13 +360,13 @@ NPL::VAnalysis* Analysis::Construct(){
// Registering the construct method to the factory //
////////////////////////////////////////////////////////////////////////////////
extern "C"{
-class proxy_analysis{
- public:
- proxy_analysis(){
- NPL::AnalysisFactory::getInstance()->SetConstructor(Analysis::Construct);
- }
-};
-
-proxy_analysis p_analysis;
+ class proxy_analysis{
+ public:
+ proxy_analysis(){
+ NPL::AnalysisFactory::getInstance()->SetConstructor(Analysis::Construct);
+ }
+ };
+
+ proxy_analysis p_analysis;
}
diff --git a/Projects/MUGAST/Analysis.h b/Projects/MUGAST/Analysis.h
index 65afe5a68..d744d46d1 100644
--- a/Projects/MUGAST/Analysis.h
+++ b/Projects/MUGAST/Analysis.h
@@ -26,6 +26,8 @@
#include"NPReaction.h"
#include"RootOutput.h"
#include"RootInput.h"
+#include "TInitialConditions.h"
+#include "TReactionConditions.h"
#include "TMust2Physics.h"
#include "TMugastPhysics.h"
#include "TCATSPhysics.h"
@@ -54,6 +56,9 @@ class Analysis: public NPL::VAnalysis{
double ELab;
double ThetaLab;
double ThetaCM;
+ double OriginalELab;
+ double OriginalThetaLab;
+
NPL::Reaction reaction;
// Energy loss table: the G4Table are generated by the simulation
NPL::EnergyLoss LightTarget;
@@ -111,6 +116,8 @@ class Analysis: public NPL::VAnalysis{
TMust2Physics* M2;
TMugastPhysics* MG;
TCATSPhysics* CATS;
-
+ bool simulation;
+ TInitialConditions* Initial;
+ TReactionConditions* ReactionConditions;
};
#endif
diff --git a/Projects/MUGAST/ShowResults.C b/Projects/MUGAST/ShowResults.C
index 15bbb392d..d5fe10231 100644
--- a/Projects/MUGAST/ShowResults.C
+++ b/Projects/MUGAST/ShowResults.C
@@ -53,66 +53,54 @@ using namespace std;
using namespace NPL;
void ShowResults(){
- // get tree
- TFile *f = new TFile("../../Outputs/Analysis/PhysicsTree.root");
- TTree *t = (TTree*) f->Get("PhysicsTree");
+ // get tree
+ TFile* f = new TFile("../../Outputs/Analysis/PhysicsTree.root");
+ TTree* t = (TTree*) f->Get("PhysicsTree");
- // draw kinematic information
- // canvas
- TCanvas *c1 = new TCanvas("c1", "kinematic information", 600, 600);
- c1->Draw();
- // kinematic line
- TH2F *hk = new TH2F("hk", "hk", 180, 0, 180, 200, 0, 60);
- hk->GetXaxis()->SetTitle("#Theta_{lab} (deg)");
- hk->GetYaxis()->SetTitle("E_{p} (MeV)");
- cout << "counts " << t->Draw("ELab:ThetaLab>>hk","Ex>0&&Ex<6","colz") << endl;
- NPL::Reaction* reaction = new NPL::Reaction();
- reaction->ReadConfigurationFile("28Mg.reaction");
- reaction->GetKinematicLine3()->Draw("c");
+ // draw kinematic information
+ // canvas
+ TCanvas *c1 = new TCanvas("c1", "Results", 1000, 1000);
+ c1->Divide(2,2);
+ c1->cd(1);
+ // kinematic line
+ TH2F* hk = new TH2F("hk", "hk", 180*3, 0, 180, 1000, 0, 60);
+ t->Draw("ELab:ThetaLab>>hk","","col");
+ hk->GetXaxis()->SetTitle("#Theta_{lab} (deg)");
+ hk->GetYaxis()->SetTitle("E_{p} (MeV)");
+ NPL::Reaction* reaction = new NPL::Reaction();
+ reaction->ReadConfigurationFile("16Odp17O_870keV_12.reaction");
+ reaction->GetKinematicLine3()->Draw("c");
- new TCanvas();
+ c1->cd(2);
+ TH1F* hEx = new TH1F("hEx", "hEx",1000, -1, 5);
+ t->Draw("Ex>>hEx","ThetaLab>100 && ThetaLab<156","col");
+ hEx->GetXaxis()->SetTitle("Ex (MeV)");
+ hEx->GetYaxis()->SetTitle("counts/100 keV");
+
+ c1->cd(3);
TH1F *hCM = new TH1F("hCM", "hCM", 36, 0, 180);
- t->Draw("ThetaCM>>hCM","Ex>0&&Ex<6","",2900);
+ t->Draw("ThetaCM>>hCM","Ex>0&&Ex<6","");
for(int i = 0 ; i < hCM->GetNbinsX();i++){
if(hCM->GetBinCenter(i)==37.5 || hCM->GetBinCenter(i)==97.5|| hCM->GetBinCenter(i)==167.5|| hCM->GetBinCenter(i)==42.5){
hCM->SetBinContent(i,0);
-
- }
-
+
}
- gPad->SetLogy();
+ }
- TFile* file= new TFile("effMUGAST.root");
- TH1* hOmega = (TH1*)file->FindObjectAny("hDetecThetaCM");
- hOmega->Rebin(5);
- hCM->Sumw2();
- hCM->Divide(hOmega);
- TGraph* g= new TGraph("22.jjj");
- g->Draw("c");
- hCM->Scale(g->Eval(32.5)/hCM->GetBinContent(hCM->FindBin(32.5)));
- TGraph* g2= new TGraph("22.l2");
- g2->Draw("c");
- TGraph* g3= new TGraph("22.l3");
- g3->Draw("c");
+
+ TCanvas *c2 = new TCanvas("c2", "Control", 1000, 1000);
+ c2->Divide(2,2);
+ c2->cd(1);
+ TH1F* hcT = new TH1F("hcT", "hcT", 180*3, -1,1);
+ t->Draw("ThetaLab-OriginalThetaLab>>hcT","","col");
+ TLine* lT = new TLine(0,0,180,180);
+ //lT->Draw();
+ c2->cd(2);
+ TH1F* hcE = new TH1F("hcE", "hcE", 1000, -1, 1);
+ t->Draw("ELab-OriginalELab>>hcE","","col");
+ TLine* lE = new TLine(0,0,60,60);
+ //lE->Draw();
}
-
-void CountingRates(Double_t ibeam = 1e5, Double_t ubt = 30){
- // load event generator file
- NPL::Reaction* reaction = new NPL::Reaction();
- reaction->ReadConfigurationFile("28Mg.reaction");
-// reaction->ReadConfigurationFile("11Be_d3He.reaction");
- // get angular distribution
- TH1F *dsig = reaction->GetCrossSectionHist();
- dsig->Draw();
- // calculate total cross section
- Double_t stot = reaction->GetTotalCrossSection();
- cout << "total cross section = " << reaction->GetTotalCrossSection() << " mb" << endl;
-
- // get target thickness
-// NPL::DetectorManager* myDetector = new NPL::DetectorManager();
-// myDetector->ReadConfigurationFile("MUGAST_Manu.detector");
-
-}
--
GitLab