From 87fc406052808a704c5b33d93b72ee74749fb2c9 Mon Sep 17 00:00:00 2001
From: fysquare <freddy.flavigny@laposte.net>
Date: Mon, 27 Apr 2020 18:54:32 +0200
Subject: [PATCH] Adding an Example (5) for QFS (Quasi-Free Scattering)
---
Examples/Example5/CheckSimu.C | 170 ++++++++++++++++++++++++
Examples/Example5/Example5.detector | 35 +++++
Examples/Example5/Example5.mac | 1 +
Examples/Example5/Example5.reaction | 29 ++++
Examples/Example5/PhysicsListOption.txt | 1 +
Examples/Example5/sim.sh | 2 +
6 files changed, 238 insertions(+)
create mode 100644 Examples/Example5/CheckSimu.C
create mode 100644 Examples/Example5/Example5.detector
create mode 100644 Examples/Example5/Example5.mac
create mode 100755 Examples/Example5/Example5.reaction
create mode 100644 Examples/Example5/PhysicsListOption.txt
create mode 100755 Examples/Example5/sim.sh
diff --git a/Examples/Example5/CheckSimu.C b/Examples/Example5/CheckSimu.C
new file mode 100644
index 000000000..237fe38a1
--- /dev/null
+++ b/Examples/Example5/CheckSimu.C
@@ -0,0 +1,170 @@
+#include <iostream>
+#include <ctime>
+#include <cstdlib>
+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"
+#include "NPQFS.h"
+
+using namespace NPL;
+
+TCanvas* canvas2;
+
+void CheckSimu(const char * fname = "Example5"){
+ // for the style
+ gStyle->SetOptStat(0);
+ gROOT->SetStyle("nptool");
+ gROOT->ForceStyle(true);
+ gStyle->SetPalette(1);
+
+ // 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");
+
+ // Connect the branches of the TTree and activate then if necessary
+ // TReactionConditions branch
+ TReactionConditions *reacCond = 0;
+ tree->SetBranchAddress("ReactionConditions", &reacCond);
+ tree->SetBranchStatus("ReactionConditions", 1);
+
+ // TInteractionCoordinates branch
+ TInteractionCoordinates *interCoord = 0;
+ tree->SetBranchAddress("InteractionCoordinates", &interCoord);
+ tree->SetBranchStatus("InteractionCoordinates", 1);
+
+ // Declare histograms
+ // for emitted particle
+ TH1F *hEmThetaCM = new TH1F("hEmThetaCM", "Light Ejectile Theta CM",180, 0, 180);
+ TH1F *hEmTheta1IF = new TH1F("hEmTheta1IF", "Ejectile1 Theta (reaction frame)", 180, 0, 180);
+ TH1F *hEmPhi1IF = new TH1F("hEmPhi1IF", "Ejectile1 Phi (reaction frame)", 360, 0, 360);
+ TH1F *hEmTheta1WF = new TH1F("hEmTheta1WF", "Ejectile1 Theta (world frame)", 180, 0, 180);
+ TH1F *hEmPhi1WF = new TH1F("hEmPhi1WF", "Ejectile1 Phi (world frame)", 360, 0, 360);
+ TH1F *hEmTheta2IF = new TH1F("hEmTheta2IF", "Ejectile2 Theta (reaction frame)", 180, 0, 180);
+ TH1F *hEmPhi2IF = new TH1F("hEmPhi2IF", "Ejectile2 Phi (reaction frame)", 360, 0, 360);
+ TH1F *hEmTheta2WF = new TH1F("hEmTheta2WF", "Ejectile2 Theta (world frame)", 180, 0, 180);
+ TH1F *hEmPhi2WF = new TH1F("hEmPhi2WF", "Ejectile2 Phi (world frame)", 360, 0, 360);
+
+ TH1F *hEmOpAngle = new TH1F("hEmOpAngle", "Opening angle betw. (1) and (2)", 100, 0, 100);
+ TH1F *hdPhi = new TH1F("hdPhi", "Phi angle difference between (1) and (2)", 200, 0, 200);
+ TH2F *hEmE1Theta1 = new TH2F("hEmE1Theta1", "Kinematics (1)", 900, 0, 90, 1000, 0, 500);
+ TH2F *hEmE2Theta2 = new TH2F("hEmE1Theta1", "Kinematics (2)", 900, 0, 90, 1000, 0, 500);
+
+ TH2F *hEmE1VsE2 = new TH2F("hEmE1VsE2", " E1 VS E2 (reaction frame)", 300, 0, 300,300,0,300);
+ TH2F *hEmTheta1VsTheta2 = new TH2F("hEmTheta1VsTheta2", " Theta1 VS Theta2 (reaction frame)", 360, 0, 90,360,0,90);
+ TH2F *hEmPhi1VsPhi2 = new TH2F("hEmPhi1VsPhi2", " Phi1 VS Phi2 (reaction frame)", 360, 0, 360,360,0,360);
+
+ // 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) {
+ cout.precision(5);
+ Double_t percent = (Double_t)i/nentries ;
+ cout << "\r Progression: " << percent*100 << " %" << flush;
+ }
+ else if (i==nentries-1) cout << "\r Progression:" << " 100%" << endl;
+
+ // Get entry
+ tree->GetEntry(i);
+
+ // Fill histos
+ // ejected particles
+ hEmThetaCM -> Fill(reacCond->GetThetaCM());
+
+ hEmTheta1IF -> Fill(reacCond->GetThetaLab_BeamFrame(0));
+ hEmTheta1WF -> Fill(reacCond->GetThetaLab_WorldFrame(0));
+ hEmE1Theta1 -> Fill(reacCond->GetThetaLab_BeamFrame(0), reacCond->GetKineticEnergy(0));
+ hEmTheta2IF -> Fill(reacCond->GetThetaLab_BeamFrame(1));
+ hEmTheta2WF -> Fill(reacCond->GetThetaLab_WorldFrame(1));
+ hEmE2Theta2 -> Fill(reacCond->GetThetaLab_BeamFrame(1), reacCond->GetKineticEnergy(1));
+
+ hEmTheta1VsTheta2 -> Fill(reacCond->GetTheta(1), reacCond->GetTheta(0));
+ hEmPhi1VsPhi2 -> Fill(reacCond->GetPhi(1), reacCond->GetPhi(0));
+ hEmE1VsE2 -> Fill(reacCond->GetKineticEnergy(1), reacCond->GetKineticEnergy(0));
+
+ double theta1 = reacCond->GetThetaLab_BeamFrame(0)*TMath::Pi()/180.;
+ double theta2 = reacCond->GetThetaLab_BeamFrame(1)*TMath::Pi()/180.;
+ double phi1 = reacCond->GetPhiLab_BeamFrame(0)*TMath::Pi()/180.;
+ double phi2 = reacCond->GetPhiLab_BeamFrame(1)*TMath::Pi()/180.;
+ double Opang = acos( sin(theta1)*sin(theta2)*cos(phi1-phi2) +
+ cos(theta1)* cos(theta2) );
+ hEmOpAngle->Fill(Opang*180./TMath::Pi());
+
+ double df = fabs(phi1-phi2)*180./TMath::Pi();
+
+ if(df>0 && df <= 180) hdPhi->Fill(df);
+ else hdPhi->Fill(360 - df);
+ }
+
+
+ // Display emmitted paricles histograms
+ canvas2 = new TCanvas("canvas2", "Emmited particles properties in reaction frame",1000,600);
+ canvas2->Divide(3,2);
+
+ canvas2->cd(1);
+ hEmThetaCM->SetXTitle("#theta_{c.m.}");
+ hEmThetaCM->SetYTitle("counts / 1^{#circ}");
+ hEmThetaCM->GetYaxis()->SetTitleOffset(1.18);
+ hEmThetaCM->GetYaxis()->SetRangeUser(0,400);
+ hEmThetaCM->Draw();
+
+ canvas2->cd(2);
+ hEmE1Theta1->SetXTitle("#theta_{1}");
+ hEmE1Theta1->SetYTitle("E_{1} (MeV)");
+ hEmE1Theta1->Draw("colz");
+
+ canvas2->cd(3);
+ hdPhi->Draw();
+ hdPhi->SetXTitle("#phi_{1} - #phi_{2} (deg)");
+ hdPhi->SetYTitle("Counts");
+
+ canvas2->cd(4);
+ hEmTheta1VsTheta2->Draw("colz");
+ hEmTheta1VsTheta2->SetXTitle("#theta_{1} (deg)");
+ hEmTheta1VsTheta2->SetYTitle("#theta_{2} (deg)");
+ NPL::QFS qfs;
+ qfs.ReadConfigurationFile("Example5.reaction");
+ qfs.SetMomentumSigma(0.);
+ TGraph* Kine = qfs.GetTheta2VsTheta1(1);
+ Kine->SetLineWidth(2);
+ Kine->SetLineColor(kRed);
+ Kine->SetLineStyle(2);
+ Kine->Draw("csame");
+
+ canvas2->cd(5);
+ hEmPhi1VsPhi2->Draw("colz");
+ hEmPhi1VsPhi2->SetXTitle("#phi_{1} (deg)");
+ hEmPhi1VsPhi2->SetYTitle("#phi_{2} (deg)");
+
+ canvas2->cd(6);
+ hEmOpAngle->Draw();
+ hEmOpAngle->SetXTitle("Opening angle (1-2) (deg)");
+ hEmOpAngle->SetYTitle("Counts");
+
+}
+
+
diff --git a/Examples/Example5/Example5.detector b/Examples/Example5/Example5.detector
new file mode 100644
index 000000000..d0c21ac27
--- /dev/null
+++ b/Examples/Example5/Example5.detector
@@ -0,0 +1,35 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Target
+ Thickness= 1 mm
+ Radius= 30 mm
+ Material= H2
+ Angle= 0 deg
+ X= 0 mm
+ Y= 0 mm
+ Z= 0 mm
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ScintillatorPlastic
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Plastic
+ X= 0 mm
+ Y= 0 mm
+ Z= 2000 mm
+ Shape= Square
+ Height= 500 mm
+ Width= 500 mm
+ Thickness= 0.001 mm
+ Scintillator= BC400
+ LeadThickness= 0 mm
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Plastic
+ X= 0 mm
+ Y= 0 mm
+ Z= 2255 mm
+ Shape= Square
+ Height= 500 mm
+ Width= 500 mm
+ Thickness= 500 mm
+ Scintillator= Pb
+ LeadThickness= 0 mm
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
diff --git a/Examples/Example5/Example5.mac b/Examples/Example5/Example5.mac
new file mode 100644
index 000000000..3a7cd691d
--- /dev/null
+++ b/Examples/Example5/Example5.mac
@@ -0,0 +1 @@
+/run/beamOn 50000
diff --git a/Examples/Example5/Example5.reaction b/Examples/Example5/Example5.reaction
new file mode 100755
index 000000000..ff84b32b4
--- /dev/null
+++ b/Examples/Example5/Example5.reaction
@@ -0,0 +1,29 @@
+%%%%%%%%%%%%%%%%%% QFS Example in progress %%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Beam
+ Particle= 12C
+ Energy= 4800 MeV
+ SigmaEnergy= 0 MeV
+ SigmaThetaX= 0 deg
+ SigmaPhiY= 0 deg
+ SigmaX= 0 mm
+ SigmaY= 0 mm
+ MeanThetaX= 0 deg
+ MeanPhiY= 0 deg
+ MeanX= 0 mm
+ MeanY= 0 mm
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+QFSReaction
+ Beam= 12C
+ Target= 1H
+ Scattered= 1H
+ KnockedOut= 1H
+ Heavy= 11B
+ ExcitationEnergyBeam= 0.0 MeV
+ ExcitationEnergyHeavy= 0.0 MeV
+ MomentumSigma= 50.0
+ ShootHeavy= 1
+ ShootLight= 1
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/Examples/Example5/PhysicsListOption.txt b/Examples/Example5/PhysicsListOption.txt
new file mode 100644
index 000000000..602b8c776
--- /dev/null
+++ b/Examples/Example5/PhysicsListOption.txt
@@ -0,0 +1 @@
+DefaultCutOff 100000000000
diff --git a/Examples/Example5/sim.sh b/Examples/Example5/sim.sh
new file mode 100755
index 000000000..7d32d779e
--- /dev/null
+++ b/Examples/Example5/sim.sh
@@ -0,0 +1,2 @@
+npsimulation -D Example5.detector -E Example5.reaction -O Example5 -B Example5.mac
+
--
GitLab