Skip to content
Snippets Groups Projects
Commit fd2bb3f7 authored by matta's avatar matta
Browse files

* Fix bug wiht NP reaction 

- Sin theta was not applied to the cross section
parent 302d0717
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Inputs/DetectorConfiguration/Sharc.detector
+ 1
− 1
View file @ fd2bb3f7
...@@ -39,7 +39,7 @@ Sharc ...@@ -39,7 +39,7 @@ Sharc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Upstream Box %Upstream Box
SharcBOX SharcBOX
Z= -30.775 Z= -34.3
ThicknessDector1= 500 ThicknessDector1= 500
ThicknessDector2= 500 ThicknessDector2= 500
ThicknessDector3= 500 ThicknessDector3= 500
......
NPAnalysis/macros/GeometricalEfficiency.C
+ 67
− 51
View file @ fd2bb3f7
...@@ -27,6 +27,7 @@ ...@@ -27,6 +27,7 @@
*****************************************************************************/ *****************************************************************************/
#include <iostream> #include <iostream>
#include <fstream>
#include <cmath> #include <cmath>
#include "TROOT.h" #include "TROOT.h"
...@@ -44,68 +45,83 @@ ...@@ -44,68 +45,83 @@
using namespace std ; using namespace std ;
void GeometricalEfficiency(const char * fname = "myResult") void GeometricalEfficiency(const char * fname = "myResult"){
{ // Open output ROOT file from NPTool simulation run
// Open output ROOT file from NPTool simulation run TString path = gSystem->Getenv("NPTOOL");
TString path = gSystem->Getenv("NPTOOL"); path += "/Outputs/Simulation/";
path += "/Outputs/Simulation/"; TString inFileName = fname;
TString inFileName = fname; if (!inFileName.Contains("root")) inFileName += ".root";
if (!inFileName.Contains("root")) inFileName += ".root"; TFile *inFile = new TFile(path + inFileName);
TFile *inFile = new TFile(path + inFileName); TTree *tree = (TTree*) inFile->Get("SimulatedTree");
TTree *tree = (TTree*) inFile->Get("SimulatedTree");
// Connect the branches of the TTree and activate then if necessary
// Connect the branches of the TTree and activate then if necessary // TInitialConditions branch
// TInitialConditions branch TInitialConditions *initCond = 0;
TInitialConditions *initCond = 0; tree->SetBranchAddress("InitialConditions", &initCond);
tree->SetBranchAddress("InitialConditions", &initCond); tree->SetBranchStatus("InitialConditions", true);
tree->SetBranchStatus("InitialConditions", true); // TInteractionCoordinates branch
// TInteractionCoordinates branch TInteractionCoordinates *interCoord = 0;
TInteractionCoordinates *interCoord = 0; tree->SetBranchAddress("InteractionCoordinates", &interCoord);
tree->SetBranchAddress("InteractionCoordinates", &interCoord); tree->SetBranchStatus("InteractionCoordinates", true);
tree->SetBranchStatus("InteractionCoordinates", true);
// Prepare histograms
// Prepare histograms TH1F *hDetecTheta = new TH1F("hDetecTheta", "DetecTheta", 180, 0, 180);
TH1F *hDetecTheta = new TH1F("hDetecTheta", "DetecTheta", 180, 0, 180); TH1F *hDetecThetaCM = new TH1F("hDetecThetaCM", "hDetecThetaCM", 90, 0, 180);
TH1F *hEmittTheta = new TH1F("hEmittTheta", "EmittTheta", 180, 0, 180); TH1F *hEmittTheta = new TH1F("hEmittTheta", "EmittTheta", 180, 0, 180);
TH1F *hEmittThetaCM = new TH1F("hEmittThetaCM", "hEmittThetaCM", 180, 0, 180);
// Read the TTree
Int_t nentries = tree->GetEntries(); // Read the TTree
// cout << "TTree contains " << nentries << " events" << endl; Int_t nentries = tree->GetEntries();
for (Int_t i = 0; i < nentries; i++) { // cout << "TTree contains " << nentries << " events" << endl;
//if (i%1000 == 0) cout << "Entry " << i << endl; for (Int_t i = 0; i < nentries; i++) {
tree->GetEntry(i); //if (i%1000 == 0) cout << "Entry " << i << endl;
// Fill histos tree->GetEntry(i);
hEmittTheta->Fill(initCond->GetThetaLab_WorldFrame(0)); // Fill histos
hEmittTheta->Fill(initCond->GetThetaLab_WorldFrame(0));
if (interCoord->GetDetectedMultiplicity() > 0) hEmittTheta->Fill(initCond->GetThetaCM(0));
hDetecTheta->Fill(interCoord->GetDetectedAngleTheta(0));
}
TCanvas* c1 = new TCanvas("c1", "c1");
// Compute relative efficiency in %
TH1F *efficiency = new TH1F("hEfficiency", "Efficiency", 180, 0, 180);
// efficiency->SetTitle("Efficiency GASPARD (Spheric version);#Theta [deg];#epsilon [%]");
efficiency->SetTitle("Efficiency GASPARD;#Theta [deg];#epsilon [%]");
efficiency->Divide(hDetecTheta,hEmittTheta,100);
efficiency->SetMaximum(110);
efficiency->Draw();
TCanvas* c2 = new TCanvas("c2", "c2"); if (interCoord->GetDetectedMultiplicity() > 0){
hEmittTheta->Draw(); hDetecTheta->Fill(interCoord->GetDetectedAngleTheta(0));
hDetecTheta->SetLineColor(kRed); hDetecThetaCM->Fill(initCond->GetThetaCM(0));
hDetecTheta->Draw("same"); }
}
TCanvas* c1 = new TCanvas("c1", "c1");
// Compute relative efficiency in %
TH1F *efficiency = new TH1F("hEfficiency", "Efficiency", 180, 0, 180);
// efficiency->SetTitle("Efficiency GASPARD (Spheric version);#Theta [deg];#epsilon [%]");
efficiency->SetTitle("Efficiency;#Theta [deg];#epsilon [%]");
efficiency->Divide(hDetecTheta,hEmittTheta,100);
efficiency->SetMaximum(110);
efficiency->Draw();
TCanvas* c2 = new TCanvas("c2", "c2");
hEmittTheta->Draw();
hDetecTheta->SetLineColor(kRed);
hDetecTheta->Draw("same");
TCanvas* c3 = new TCanvas("c3", "c3"); TCanvas* c3 = new TCanvas("c3", "c3");
TH1F* SolidA = new TH1F(*hDetecTheta); TH1F* SolidA = new TH1F(*hDetecTheta);
SolidA->Sumw2(); SolidA->Sumw2();
TF1* C = new TF1("C",Form("%i /(4*%f)",nentries,M_PI),0,180); TF1* C = new TF1("C",Form("%i /(4*%f)",nentries,M_PI),0,180);
SolidA->Divide(C,1); SolidA->Divide(C,1);
SolidA->Rebin(2); // SolidA->Rebin(2);
SolidA->Draw(); SolidA->Draw();
TF1* f = new TF1("f",Form("2 * %f * sin(x*%f/180.) *2*%f/180.",M_PI,M_PI,M_PI),0,180); TF1* f = new TF1("f",Form("2 * %f * sin(x*%f/180.) *2*%f/180.",M_PI,M_PI,M_PI),0,180);
f->Draw("SAME"); f->Draw("SAME");
SolidA->GetXaxis()->SetTitle("#theta_{Lab} (deg)"); SolidA->GetXaxis()->SetTitle("#theta_{Lab} (deg)");
SolidA->GetYaxis()->SetTitle("d#Omega (sr) "); SolidA->GetYaxis()->SetTitle("d#Omega (sr) ");
c3->Update(); c3->Update();
ofstream output;
output.open("data.txt");
for(int i = 1 ; i < 90 ; i++)
output << hDetecThetaCM->GetBinCenter(i) << " "
<< hDetecThetaCM->GetBinContent(i)<< endl ;
output.close();
} }
NPLib/Physics/NPReaction.cxx
+ 11
− 1
View file @ fd2bb3f7
...@@ -52,6 +52,10 @@ ...@@ -52,6 +52,10 @@
using namespace NPUNITS; using namespace NPUNITS;
#endif #endif
// ROOT
#include"TF1.h"
ClassImp(Reaction) ClassImp(Reaction)
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
...@@ -372,7 +376,13 @@ void Reaction::ReadConfigurationFile(string Path){ ...@@ -372,7 +376,13 @@ void Reaction::ReadConfigurationFile(string Path){
string FileName,HistName; string FileName,HistName;
ReactionFile >> FileName >> HistName; ReactionFile >> FileName >> HistName;
if(fVerboseLevel==1) cout << "Reading Cross Section file: " << FileName << endl; if(fVerboseLevel==1) cout << "Reading Cross Section file: " << FileName << endl;
SetCrossSectionHist( Read1DProfile(FileName, HistName )); TH1F* CStemp = Read1DProfile(FileName, HistName );
// multiply CStemp by sin(theta)
TF1* fsin = new TF1("sin",Form("1/(sin(x*%f/180.))",M_PI),0,180);
CStemp->Divide(fsin,1);
SetCrossSectionHist(CStemp);
delete fsin;
} }
else if (DataBuffer.compare(0, 17, "HalfOpenAngleMin=") == 0) { else if (DataBuffer.compare(0, 17, "HalfOpenAngleMin=") == 0) {
......
NPSimulation/Simulation.cc
+ 92
− 84
View file @ fd2bb3f7
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#include "G4RunManager.hh" #include "G4RunManager.hh"
#include "G4PhysListFactory.hh"
// UI // UI
#include "G4UImanager.hh" #include "G4UImanager.hh"
#include "G4UIterminal.hh" #include "G4UIterminal.hh"
...@@ -32,84 +32,92 @@ ...@@ -32,84 +32,92 @@
#include "NPOptionManager.h" #include "NPOptionManager.h"
#include "RootOutput.h" #include "RootOutput.h"
int main(int argc, char** argv) int main(int argc, char** argv){
{ // Initialize NPOptionManager object
// Initialize NPOptionManager object NPOptionManager* OptionManager = NPOptionManager::getInstance(argc, argv);
NPOptionManager* OptionManager = NPOptionManager::getInstance(argc, argv); // Test if input files are found. If not, exit
// Test if input files are found. If not, exit if (OptionManager->IsDefault("EventGenerator")) OptionManager->SendErrorAndExit("EventGenerator");
if (OptionManager->IsDefault("EventGenerator")) OptionManager->SendErrorAndExit("EventGenerator"); if (OptionManager->IsDefault("DetectorConfiguration")) OptionManager->SendErrorAndExit("DetectorConfiguration");
if (OptionManager->IsDefault("DetectorConfiguration")) OptionManager->SendErrorAndExit("DetectorConfiguration"); // case when input files are here
// case when input files are here G4String EventGeneratorFileName = OptionManager->GetReactionFile();
G4String EventGeneratorFileName = OptionManager->GetReactionFile(); G4String DetectorFileName = OptionManager->GetDetectorFile();
G4String DetectorFileName = OptionManager->GetDetectorFile();
// my Verbose output class
// my Verbose output class G4VSteppingVerbose::SetInstance(new SteppingVerbose);
G4VSteppingVerbose::SetInstance(new SteppingVerbose);
// Construct the default run manager
// Construct the default run manager G4RunManager* runManager = new G4RunManager;
G4RunManager* runManager = new G4RunManager;
// set mandatory initialization classes
// set mandatory initialization classes DetectorConstruction* detector = new DetectorConstruction();
DetectorConstruction* detector = new DetectorConstruction(); runManager->SetUserInitialization(detector);
runManager->SetUserInitialization(detector);
PhysicsList* physicsList = new PhysicsList();
PhysicsList* physics = new PhysicsList(); runManager->SetUserInitialization(physicsList);
runManager->SetUserInitialization(physics);
PrimaryGeneratorAction* primary = new PrimaryGeneratorAction(detector); // Test for Built in physics list
// G4PhysListFactory *physListFactory = new G4PhysListFactory();
// Initialize Geant4 kernel //G4VUserPhysicsList *physicsList =
runManager->Initialize(); //physListFactory->GetReferencePhysList("QGSP_BERT");
physics->MyOwnConstruction();
runManager->SetUserInitialization(physicsList);
///////////////////////////////////////////////////////////////
///////////////// Initializing the Root Output ////////////////
/////////////////////////////////////////////////////////////// PrimaryGeneratorAction* primary = new PrimaryGeneratorAction(detector);
RootOutput::getInstance("Simulation/" + OptionManager->GetOutputFile());
// Initialize Geant4 kernel
/////////////////////////////////////////////////////////////// runManager->Initialize();
////////////// Reading Detector Configuration ///////////////// physicsList->MyOwnConstruction();
///////////////////////////////////////////////////////////////
detector->ReadConfigurationFile(DetectorFileName); ///////////////////////////////////////////////////////////////
///////////////// Initializing the Root Output ////////////////
/////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////
////////////////////// Reading Reaction /////////////////////// RootOutput::getInstance("Simulation/" + OptionManager->GetOutputFile());
///////////////////////////////////////////////////////////////
primary->ReadEventGeneratorFile(EventGeneratorFileName); ///////////////////////////////////////////////////////////////
runManager->SetUserAction(primary); ////////////// Reading Detector Configuration /////////////////
///////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////// detector->ReadConfigurationFile(DetectorFileName);
////////////////// Starting the Event Action //////////////////
/////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////
EventAction* event_action = new EventAction() ; ////////////////////// Reading Reaction ///////////////////////
event_action->SetDetector(detector) ; ///////////////////////////////////////////////////////////////
runManager->SetUserAction(event_action) ; primary->ReadEventGeneratorFile(EventGeneratorFileName);
runManager->SetUserAction(primary);
///////////////////////////////////////////////////////////////
/////// Get the pointer to the User Interface manager //////// ///////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////// ////////////////// Starting the Event Action //////////////////
// G4UImanager* UI = G4UImanager::GetUIpointer(); ///////////////////////////////////////////////////////////////
EventAction* event_action = new EventAction() ;
/////////////////////////////////////////////////////////////// event_action->SetDetector(detector) ;
/////////// Define UI terminal for interactive mode /////////// runManager->SetUserAction(event_action) ;
///////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////
/////// Get the pointer to the User Interface manager ////////
///////////////////////////////////////////////////////////////
// G4UImanager* UI = G4UImanager::GetUIpointer();
///////////////////////////////////////////////////////////////
/////////// Define UI terminal for interactive mode ///////////
///////////////////////////////////////////////////////////////
#ifdef G4VIS_USE #ifdef G4VIS_USE
G4VisManager* visManager = new G4VisExecutive("Quiet"); G4VisManager* visManager = new G4VisExecutive("Quiet");
visManager->Initialize(); visManager->Initialize();
#endif #endif
/* G4UIsession* session = 0; /* G4UIsession* session = 0;
#ifdef G4UI_USE_TCSH #ifdef G4UI_USE_TCSH
session = new G4UIterminal(new G4UItcsh); session = new G4UIterminal(new G4UItcsh);
#else #else
session = new G4UIterminal(); session = new G4UIterminal();
#endif #endif
UI->ApplyCommand("/control/execute vis.mac"); UI->ApplyCommand("/control/execute vis.mac");
session->SessionStart(); session->SessionStart();
delete session; delete session;
*/ */
// interactive mode : define UI session // interactive mode : define UI session
// Get the pointer to the User Interface manager // Get the pointer to the User Interface manager
G4UImanager* UImanager = G4UImanager::GetUIpointer(); G4UImanager* UImanager = G4UImanager::GetUIpointer();
...@@ -123,18 +131,18 @@ int main(int argc, char** argv) ...@@ -123,18 +131,18 @@ int main(int argc, char** argv)
ui->SessionStart(); ui->SessionStart();
delete ui; delete ui;
#endif #endif
#ifdef G4VIS_USE #ifdef G4VIS_USE
delete visManager; delete visManager;
#endif #endif
/////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////
////////////////////// Job termination //////////////////////// ////////////////////// Job termination ////////////////////////
/////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////
// delete primary; delete detector; // delete primary; delete detector;
delete runManager; delete runManager;
NPOptionManager::getInstance()->Destroy(); NPOptionManager::getInstance()->Destroy();
RootOutput::getInstance()->Destroy(); RootOutput::getInstance()->Destroy();
return 0; return 0;
} }
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment