From bcf9f8ae1906ec203f992989756d25c28f6d4dc5 Mon Sep 17 00:00:00 2001
From: Adrien Matta <matta@lpccaen.in2p3.fr>
Date: Mon, 11 Sep 2023 09:48:02 +0200
Subject: [PATCH] * Fixing GetTheta calling in ShowResult of Gaspard benchmarck
---
Benchmarks/gaspard/ShowResults.C | 202 ++++++++++++++++---------------
1 file changed, 102 insertions(+), 100 deletions(-)
diff --git a/Benchmarks/gaspard/ShowResults.C b/Benchmarks/gaspard/ShowResults.C
index 053a03ba1..da6e53501 100644
--- a/Benchmarks/gaspard/ShowResults.C
+++ b/Benchmarks/gaspard/ShowResults.C
@@ -26,126 +26,131 @@
*****************************************************************************/
// C++ headers
-#include <iostream>
-#include <ctime>
#include <cstdlib>
+#include <ctime>
+#include <iostream>
using namespace std;
-// ROOT headers
-#include "TROOT.h"
-#include "TStyle.h"
-#include "TCanvas.h"
-#include "TSystem.h"
-#include "TFile.h"
-#include "TString.h"
-#include "TEllipse.h"
-#include "TLegend.h"
-#include "TTree.h"
-#include "TBranch.h"
-#include "TH1F.h"
-#include "TH2F.h"
-#include "TGraph.h"
-
-// nptool headers
-#include "TReactionConditions.h"
-#include "TInteractionCoordinates.h"
-#include "NPReaction.h"
+//// ROOT headers
+// #include "TROOT.h"
+// #include "TStyle.h"
+// #include "TCanvas.h"
+// #include "TSystem.h"
+// #include "TFile.h"
+// #include "TString.h"
+// #include "TEllipse.h"
+// #include "TLegend.h"
+// #include "TTree.h"
+// #include "TBranch.h"
+// #include "TH1F.h"
+// #include "TH2F.h"
+// #include "TGraph.h"
+//
+//// nptool headers
+// #include "TReactionConditions.h"
+// #include "TInteractionCoordinates.h"
+// #include "NPReaction.h"
using namespace NPL;
TCanvas* canvas1;
TCanvas* canvas2;
-void ShowResults(const char * fname = "benchmark_gaspard"){
- // for the style
+void ShowResults(const char* fname = "benchmark_gaspard") {
+ // for the style
gStyle->SetOptStat(0);
- gROOT->SetStyle("nptool");
- gROOT->ForceStyle(true);
+ gROOT->SetStyle("nptool");
+ gROOT->ForceStyle(true);
// Open output ROOT file from NPTool simulation run
TString path = gSystem->Getenv("NPTOOL");
path += "/Outputs/Simulation/";
TString inFileName = fname;
- if (!inFileName.Contains("root")) inFileName += ".root";
- TFile *inFile = new TFile(path + inFileName);
- if (!inFile->IsOpen()) exit(1);
- TTree *tree = (TTree*) inFile->Get("SimulatedTree");
+ if (!inFileName.Contains("root"))
+ inFileName += ".root";
+ TFile* inFile = new TFile(path + inFileName);
+ if (!inFile->IsOpen())
+ exit(1);
+ TTree* tree = (TTree*)inFile->Get("SimulatedTree");
// Connect the branches of the TTree and activate then if necessary
// TReactionConditions branch
- TReactionConditions *reacCond = 0;
+ TReactionConditions* reacCond = 0;
tree->SetBranchAddress("ReactionConditions", &reacCond);
tree->SetBranchStatus("ReactionConditions", 1);
-
+
// TInteractionCoordinates branch
- TInteractionCoordinates *interCoord = 0;
+ TInteractionCoordinates* interCoord = 0;
tree->SetBranchAddress("InteractionCoordinates", &interCoord);
tree->SetBranchStatus("InteractionCoordinates", 1);
// Declare histograms
// for incident beam
- TH2F *hEmittanceXY = new TH2F("hEmittanceXY", "Beam Position on Target (x, y)", 200, -10, 10, 200, -10, 10);
- TH2F *hEmittanceXTheta = new TH2F("hEmittanceXTheta", "Beam Emittance (#theta_{x}-90,x)", 200, -10, 10, 200, -2, 2);
- TH2F *hEmittanceYPhi = new TH2F("hEmittanceYPhi", "Beam Emittance (#phi_{y}-90,y)", 200, -10, 10, 200, -2, 2);
- TH1F *hIncidentTheta = new TH1F("hIncidentTheta", "Beam Incident Theta (Polar angle)", 100, 0, 2);
- TH1F *hIncidentPhi = new TH1F("hIncidentPhi", "Beam Incident Phi (Azimuthal angle)", 370, -185, 185);
- TH1F *hIncidentZ = new TH1F("hIncidentZ", "Beam z-position of interaction", 200, -0.007, 0.007);
+ TH2F* hEmittanceXY = new TH2F("hEmittanceXY", "Beam Position on Target (x, y)", 200, -10, 10, 200, -10, 10);
+ TH2F* hEmittanceXTheta = new TH2F("hEmittanceXTheta", "Beam Emittance (#theta_{x}-90,x)", 200, -10, 10, 200, -2, 2);
+ TH2F* hEmittanceYPhi = new TH2F("hEmittanceYPhi", "Beam Emittance (#phi_{y}-90,y)", 200, -10, 10, 200, -2, 2);
+ TH1F* hIncidentTheta = new TH1F("hIncidentTheta", "Beam Incident Theta (Polar angle)", 100, 0, 2);
+ TH1F* hIncidentPhi = new TH1F("hIncidentPhi", "Beam Incident Phi (Azimuthal angle)", 370, -185, 185);
+ TH1F* hIncidentZ = new TH1F("hIncidentZ", "Beam z-position of interaction", 200, -0.007, 0.007);
// for emitted particle
- TH1F *hEmittedThetaCM = new TH1F("hEmittedThetaCM", "Light Ejectile Theta CM", 180, 0, 180);
- TH1F *hEmittedThetaIF = new TH1F("hEmittedThetaIF", "Light Ejectile Theta (reaction frame)", 180, 0, 180);
- TH1F *hEmittedPhiIF = new TH1F("hEmittedPhiIF", "Light Ejectile Phi (reaction frame)", 360, 0, 360);
- TH1F *hEmittedThetaWF = new TH1F("hEmittedThetaWF", "Light Ejectile Theta (detector frame)", 180, 0, 180);
- TH1F *hEmittedPhiWF = new TH1F("hEmittedPhiWF", "Light Ejectile Phi (detector frame)", 360, 0, 360);
- TH1F *hControlPhi = new TH1F("hControlPhi", "Light Ejectile IncidentPhi + EmittedPhi", 360, 0, 360);
- TH2F *hEmittedETheta = new TH2F("hEmittedETheta", "Kinematic line", 1800, 0, 180, 1800, 0, 60);
-
+ TH1F* hEmittedThetaCM = new TH1F("hEmittedThetaCM", "Light Ejectile Theta CM", 180, 0, 180);
+ TH1F* hEmittedThetaIF = new TH1F("hEmittedThetaIF", "Light Ejectile Theta (reaction frame)", 180, 0, 180);
+ TH1F* hEmittedPhiIF = new TH1F("hEmittedPhiIF", "Light Ejectile Phi (reaction frame)", 360, 0, 360);
+ TH1F* hEmittedThetaWF = new TH1F("hEmittedThetaWF", "Light Ejectile Theta (detector frame)", 180, 0, 180);
+ TH1F* hEmittedPhiWF = new TH1F("hEmittedPhiWF", "Light Ejectile Phi (detector frame)", 360, 0, 360);
+ TH1F* hControlPhi = new TH1F("hControlPhi", "Light Ejectile IncidentPhi + EmittedPhi", 360, 0, 360);
+ TH2F* hEmittedETheta = new TH2F("hEmittedETheta", "Kinematic line", 1800, 0, 180, 1800, 0, 60);
// for emitted particle but gated on interaction in detector
- TH2F *hEmittedETheta_detected = new TH2F("hEmittedETheta_detected", "Kinematic line (interaction in det.)",1800, 0, 180, 1800, 0, 60);
- TH1F *hEmittedThetaIF_detected = new TH1F("hEmittedThetaIF_detected", "Light Ejectile Theta in reaction frame (interaction in det.)",180, 0, 180);
- TH1F *hEmittedThetaCM_detected = new TH1F("hEmittedThetaCM_detected", "Light Ejectile ThetaCM (interaction in det.)",180, 0, 180);
+ TH2F* hEmittedETheta_detected =
+ new TH2F("hEmittedETheta_detected", "Kinematic line (interaction in det.)", 1800, 0, 180, 1800, 0, 60);
+ TH1F* hEmittedThetaIF_detected =
+ new TH1F("hEmittedThetaIF_detected", "Light Ejectile Theta in reaction frame (interaction in det.)", 180, 0, 180);
+ TH1F* hEmittedThetaCM_detected =
+ new TH1F("hEmittedThetaCM_detected", "Light Ejectile ThetaCM (interaction in det.)", 180, 0, 180);
// Read the TTree
Int_t nentries = tree->GetEntries();
cout << endl << " TTree contains " << nentries << " events" << endl;
for (Int_t i = 0; i < nentries; i++) {
- if (i%10000 == 0 && i!=0) {
+ if (i % 10000 == 0 && i != 0) {
cout.precision(5);
- Double_t percent = (Double_t)i/nentries ;
- cout << "\r Progression: " << percent*100 << " %" << flush;
+ Double_t percent = (Double_t)i / nentries;
+ cout << "\r Progression: " << percent * 100 << " %" << flush;
}
- else if (i==nentries-1) cout << "\r Progression:" << " 100%" << endl;
+ else if (i == nentries - 1)
+ cout << "\r Progression:"
+ << " 100%" << endl;
// Get entry
tree->GetEntry(i);
// Fill histos
// incident beam
- hEmittanceXY -> Fill(reacCond->GetVertexPositionX(), reacCond->GetVertexPositionY());
- hIncidentZ -> Fill(reacCond->GetVertexPositionZ());
- hEmittanceXTheta -> Fill(reacCond->GetVertexPositionX(), reacCond->GetBeamEmittanceThetaX()-90);
- hEmittanceYPhi -> Fill(reacCond->GetVertexPositionY(), reacCond->GetBeamEmittancePhiY()-90);
- hIncidentTheta -> Fill(reacCond->GetBeamEmittanceTheta());
- hIncidentPhi -> Fill(reacCond->GetBeamEmittancePhi());
+ hEmittanceXY->Fill(reacCond->GetVertexPositionX(), reacCond->GetVertexPositionY());
+ hIncidentZ->Fill(reacCond->GetVertexPositionZ());
+ hEmittanceXTheta->Fill(reacCond->GetVertexPositionX(), reacCond->GetBeamEmittanceThetaX() - 90);
+ hEmittanceYPhi->Fill(reacCond->GetVertexPositionY(), reacCond->GetBeamEmittancePhiY() - 90);
+ hIncidentTheta->Fill(reacCond->GetBeamEmittanceTheta());
+ hIncidentPhi->Fill(reacCond->GetBeamEmittancePhi());
// ejected particle
- hEmittedThetaCM -> Fill(reacCond->GetThetaCM());
- hEmittedThetaIF -> Fill(reacCond->GetThetaLab_BeamFrame(0));
- hEmittedThetaWF -> Fill(reacCond->GetThetaLab_WorldFrame(0));
- hEmittedETheta -> Fill(reacCond->GetThetaLab_BeamFrame(0), reacCond->GetKineticEnergy(0));
+ hEmittedThetaCM->Fill(reacCond->GetThetaCM());
+ hEmittedThetaIF->Fill(reacCond->GetTheta_BeamFrame(0));
+ hEmittedThetaWF->Fill(reacCond->GetTheta_WorldFrame(0));
+ hEmittedETheta->Fill(reacCond->GetTheta_BeamFrame(0), reacCond->GetKineticEnergy(0));
if (interCoord->GetDetectedMultiplicity() > 0) {
- hEmittedETheta_detected -> Fill(reacCond->GetThetaLab_BeamFrame(0), reacCond->GetKineticEnergy(0));
- hEmittedThetaIF_detected -> Fill(reacCond->GetThetaLab_BeamFrame(0));
- hEmittedThetaCM_detected -> Fill(reacCond->GetThetaCM());
+ hEmittedETheta_detected->Fill(reacCond->GetTheta_BeamFrame(0), reacCond->GetKineticEnergy(0));
+ hEmittedThetaIF_detected->Fill(reacCond->GetTheta_BeamFrame(0));
+ hEmittedThetaCM_detected->Fill(reacCond->GetThetaCM());
}
}
-
// Display beam histograms
- canvas1 = new TCanvas("canvas1", "Incident beam properties from EventGenerator", 700,0, 500,750);
- canvas1->Divide(2,3);
+ canvas1 = new TCanvas("canvas1", "Incident beam properties from EventGenerator", 700, 0, 500, 750);
+ canvas1->Divide(2, 3);
canvas1->cd(1);
hEmittanceXTheta->Draw("colz");
@@ -155,7 +160,7 @@ void ShowResults(const char * fname = "benchmark_gaspard"){
canvas1->cd(2);
hEmittanceYPhi->Draw("colz");
hEmittanceYPhi->SetXTitle("Y (mm)");
- hEmittanceYPhi->SetYTitle("#Phi_{Y} (deg)");
+ hEmittanceYPhi->SetYTitle("#Phi_{Y} (deg)");
canvas1->cd(3);
hIncidentTheta->Draw();
@@ -164,16 +169,16 @@ void ShowResults(const char * fname = "benchmark_gaspard"){
canvas1->cd(4);
hIncidentPhi->Draw();
hIncidentPhi->SetXTitle("Incident #phi (deg)");
-
+
canvas1->cd(5);
hEmittanceXY->Draw("colz");
hEmittanceXY->SetXTitle("Position on Target X (mm)");
- hEmittanceXY->SetYTitle("Position on Target Y (mm)");
- TEllipse *target = new TEllipse(0,0,7.5,7.5);
+ hEmittanceXY->SetYTitle("Position on Target Y (mm)");
+ TEllipse* target = new TEllipse(0, 0, 7.5, 7.5);
target->SetFillStyle(0000);
target->SetLineStyle(2);
target->Draw("same");
- TEllipse *targetCenter = new TEllipse(1.0,2.0,0.2,0.2);
+ TEllipse* targetCenter = new TEllipse(1.0, 2.0, 0.2, 0.2);
targetCenter->Draw("same");
canvas1->cd(6);
@@ -181,24 +186,24 @@ void ShowResults(const char * fname = "benchmark_gaspard"){
hIncidentZ->SetXTitle("Z (mm)");
// Display detector histograms
- canvas2 = new TCanvas("canvas2", "Emitted particle properties",500,500);
- canvas2->Divide(2,2);
+ canvas2 = new TCanvas("canvas2", "Emitted particle properties", 500, 500);
+ canvas2->Divide(2, 2);
canvas2->cd(1);
hEmittedThetaCM->SetXTitle("#theta_{c.m.}");
hEmittedThetaCM->SetYTitle("counts / 1^{#circ}");
- hEmittedThetaCM->SetLineColor(kBlue-3);
- hEmittedThetaCM->SetFillColor(kBlue-3);
+ hEmittedThetaCM->SetLineColor(kBlue - 3);
+ hEmittedThetaCM->SetFillColor(kBlue - 3);
hEmittedThetaCM->GetYaxis()->SetTitleOffset(1.18);
hEmittedThetaCM->Draw();
-
- hEmittedThetaCM_detected->SetLineColor(kOrange-3);
- hEmittedThetaCM_detected->SetFillColor(kOrange-3);
+
+ hEmittedThetaCM_detected->SetLineColor(kOrange - 3);
+ hEmittedThetaCM_detected->SetFillColor(kOrange - 3);
hEmittedThetaCM_detected->Draw("same");
-
+
TLegend* leg = new TLegend(0.53, 0.68, 0.81, 0.9);
- leg->AddEntry(hEmittedThetaCM,"Emitted","f");
- leg->AddEntry(hEmittedThetaCM_detected,"Detected","f");
+ leg->AddEntry(hEmittedThetaCM, "Emitted", "f");
+ leg->AddEntry(hEmittedThetaCM_detected, "Detected", "f");
leg->Draw();
canvas2->cd(2);
@@ -206,48 +211,45 @@ void ShowResults(const char * fname = "benchmark_gaspard"){
hEmittedETheta_detected->SetYTitle("E_{Lab} (MeV)");
hEmittedETheta_detected->Draw();
NPL::Reaction r("132Sn(d,p)133Sn@1320");
- TGraph* Kine = r.GetKinematicLine3();
+ TGraph* Kine = r.GetKinematicLine3();
Kine->SetLineWidth(1);
Kine->SetLineStyle(2);
- Kine->SetLineColor(kOrange-3);
+ Kine->SetLineColor(kOrange - 3);
Kine->Draw("c same");
canvas2->cd(3);
hEmittedThetaIF->SetXTitle("#theta_{Lab}");
hEmittedThetaIF->SetYTitle("counts / 1^{#circ}");
- hEmittedThetaIF->SetLineColor(kBlue-3);
- hEmittedThetaIF->SetFillColor(kBlue-3);
+ hEmittedThetaIF->SetLineColor(kBlue - 3);
+ hEmittedThetaIF->SetFillColor(kBlue - 3);
hEmittedThetaIF->GetYaxis()->SetTitleOffset(1.18);
hEmittedThetaIF->Draw();
- hEmittedThetaIF_detected->SetLineColor(kOrange-3);
- hEmittedThetaIF_detected->SetFillColor(kOrange-3);
+ hEmittedThetaIF_detected->SetLineColor(kOrange - 3);
+ hEmittedThetaIF_detected->SetFillColor(kOrange - 3);
hEmittedThetaIF_detected->Draw("same");
canvas2->cd(4);
- TH1F *hEfficiency = new TH1F("hEfficiency", "Efficiency", 180, 0, 180);
+ TH1F* hEfficiency = new TH1F("hEfficiency", "Efficiency", 180, 0, 180);
hEfficiency->Divide(hEmittedThetaIF_detected, hEmittedThetaIF, 100, 1);
hEfficiency->GetXaxis()->SetTitle("#theta (deg)");
hEfficiency->GetYaxis()->SetTitle("#epsilon (%)");
hEfficiency->Draw();
-
TFile* referenceFile = new TFile("reference.root");
- TCanvas* canvas1_ref = (TCanvas*) referenceFile->FindObjectAny("canvas1_ref");
- canvas1_ref->SetTitle( Form("FROM reference.root: %s",canvas1_ref->GetTitle()) );
- TCanvas* canvas2_ref = (TCanvas*) referenceFile->FindObjectAny("canvas2_ref");
- canvas2_ref->SetTitle( Form("FROM reference.root: %s",canvas2_ref->GetTitle()) );
+ TCanvas* canvas1_ref = (TCanvas*)referenceFile->FindObjectAny("canvas1_ref");
+ canvas1_ref->SetTitle(Form("FROM reference.root: %s", canvas1_ref->GetTitle()));
+ TCanvas* canvas2_ref = (TCanvas*)referenceFile->FindObjectAny("canvas2_ref");
+ canvas2_ref->SetTitle(Form("FROM reference.root: %s", canvas2_ref->GetTitle()));
canvas1_ref->Draw();
canvas2_ref->Draw();
-
}
-
////////////////////////////////////////////////////////////////////////////////
// Use this method to overwrite the reference file only
// DO NOT USE UNLESS YOU WANT TO MAKE A CHANGE TO THE BENCHMARK
-void WriteGaspardReference(){
- TFile *outFile = new TFile("reference.root","RECREATE");
+void WriteGaspardReference() {
+ TFile* outFile = new TFile("reference.root", "RECREATE");
canvas1->SetName("canvas1_ref");
canvas2->SetName("canvas2_ref");
canvas1->Write();
--
GitLab