diff --git a/Benchmarks/gaspard/ShowResults.C b/Benchmarks/gaspard/ShowResults.C
index 053a03ba11fbb7666dd9dc1ab2015c7e8d38428f..da6e5350144a16896ccec0507c32ed34d4862ad9 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();
diff --git a/NPLib/Detectors/Exogam/TExogamData.cxx b/NPLib/Detectors/Exogam/TExogamData.cxx
index 448d049ce88ad70235f60a432e82bb4dd7db1e22..514a5059f828a26b88b709b0123f45d379b8671b 100644
--- a/NPLib/Detectors/Exogam/TExogamData.cxx
+++ b/NPLib/Detectors/Exogam/TExogamData.cxx
@@ -43,21 +43,21 @@ TExogamData::~TExogamData()
void TExogamData::Clear()
{
- fEXO_E.clear();
- fEXO_E_CrystalNbr.clear();
- fEXO_E_TS.clear();
- fEXO_HG.clear();
- fEXO_HG_CrystalNbr.clear();
- fEXO_HG_TS.clear();
- fEXO_TDC.clear();
- fEXO_TDC_CrystalNbr.clear();
- fEXO_TDC_TS.clear();
- fEXO_Outer.clear();
- fEXO_Outer_SubCrystalNbr.clear();
- fEXO_BGO.clear();
- fEXO_BGO_CrystalNbr.clear();
- fEXO_CSI.clear();
- fEXO_CSI_CrystalNbr.clear();
+ fExoE.clear();
+ fExoE_CrystalNbr.clear();
+ fExoE_TS.clear();
+ fExoHG.clear();
+ fExoHG_CrystalNbr.clear();
+ fExoHG_TS.clear();
+ fExoTDC.clear();
+ fExoTDC_CrystalNbr.clear();
+ fExoTDC_TS.clear();
+ fExoOuter.clear();
+ fExoOuter_SubCrystalNbr.clear();
+ fExoBGO.clear();
+ fExoBGO_CrystalNbr.clear();
+ fExoCsI.clear();
+ fExoCsI_CrystalNbr.clear();
}
@@ -66,28 +66,28 @@ void TExogamData::Dump() const
{
cout << "XXXXXXXXXXXXXXXXXXXXXXXX New Event XXXXXXXXXXXXXXXXX" << endl;
- cout << "Inner6MV Mult = " << fEXO_E.size() << endl;
- for (UShort_t i = 0; i < fEXO_E.size(); i++) {
- cout << "Energy: " << fEXO_E[i] << " Cristal Numb: " << fEXO_E_CrystalNbr[i] << " TimeStamp: " << fEXO_E_TS[i] << endl;
+ cout << "Inner6MV Mult = " << fExoE.size() << endl;
+ for (UShort_t i = 0; i < fExoE.size(); i++) {
+ cout << "Energy: " << fExoE[i] << " Cristal Numb: " << fExoE_CrystalNbr[i] << " TimeStamp: " << fExoE_TS[i] << endl;
}
- cout << "Inner20MV Mult = " << fEXO_HG.size() << endl;
- for (UShort_t i = 0; i < fEXO_HG.size(); i++) {
- cout << "Energy: " << fEXO_HG[i] << " Cristal Numb: " << fEXO_HG_CrystalNbr[i] << " TimeStamp: " << fEXO_HG_TS[i] << endl;
+ cout << "Inner20MV Mult = " << fExoHG.size() << endl;
+ for (UShort_t i = 0; i < fExoHG.size(); i++) {
+ cout << "Energy: " << fExoHG[i] << " Cristal Numb: " << fExoHG_CrystalNbr[i] << " TimeStamp: " << fExoHG_TS[i] << endl;
}
- cout << "OutersV Mult = " << fEXO_Outer.size() << endl;
- for (UShort_t i = 0; i < fEXO_Outer.size(); i++) {
- cout << "Energy: " << fEXO_Outer[i] << " Cristal Numb: " << fEXO_Outer_SubCrystalNbr[i] << endl;
+ cout << "OutersV Mult = " << fExoOuter.size() << endl;
+ for (UShort_t i = 0; i < fExoOuter.size(); i++) {
+ cout << "Energy: " << fExoOuter[i] << " Cristal Numb: " << fExoOuter_SubCrystalNbr[i] << endl;
}
- cout << "DeltaTV Mult = " << fEXO_TDC.size() << endl;
- for (UShort_t i = 0; i < fEXO_TDC.size(); i++) {
- cout << "Energy: " << fEXO_TDC[i] << " Cristal Numb: " << fEXO_TDC_CrystalNbr[i] << " TimeStamp: " << fEXO_TDC_TS[i] << endl;
+ cout << "DeltaTV Mult = " << fExoTDC.size() << endl;
+ for (UShort_t i = 0; i < fExoTDC.size(); i++) {
+ cout << "Energy: " << fExoTDC[i] << " Cristal Numb: " << fExoTDC_CrystalNbr[i] << " TimeStamp: " << fExoTDC_TS[i] << endl;
}
- cout << "BGOV Mult = " << fEXO_BGO.size() << endl;
- for (UShort_t i = 0; i < fEXO_BGO.size(); i++) {
- cout << "Energy: " << fEXO_BGO[i] << " Cristal Numb: " << fEXO_BGO_CrystalNbr[i] << endl;
+ cout << "BGOV Mult = " << fExoBGO.size() << endl;
+ for (UShort_t i = 0; i < fExoBGO.size(); i++) {
+ cout << "Energy: " << fExoBGO[i] << " Cristal Numb: " << fExoBGO_CrystalNbr[i] << endl;
}
- cout << "CSIV Mult = " << fEXO_CSI.size() << endl;
- for (UShort_t i = 0; i < fEXO_CSI.size(); i++) {
- cout << "Energy: " << fEXO_CSI[i] << " Cristal Numb: " << fEXO_CSI_CrystalNbr[i] << endl;
+ cout << "CsIV Mult = " << fExoCsI.size() << endl;
+ for (UShort_t i = 0; i < fExoCsI.size(); i++) {
+ cout << "Energy: " << fExoCsI[i] << " Cristal Numb: " << fExoCsI_CrystalNbr[i] << endl;
}
}
diff --git a/NPLib/Detectors/Exogam/TExogamData.h b/NPLib/Detectors/Exogam/TExogamData.h
index 38d242b64523eea117d4dbff42fa73d5d64e52e1..10e3ba2c513baa0002f69c3d41c44d682024f480 100644
--- a/NPLib/Detectors/Exogam/TExogamData.h
+++ b/NPLib/Detectors/Exogam/TExogamData.h
@@ -28,26 +28,26 @@
class TExogamData : public TObject {
private:
- std::vector<UShort_t> fEXO_E;
- std::vector<UShort_t> fEXO_E_CrystalNbr;
- std::vector<ULong64_t> fEXO_E_TS;
+ std::vector<UShort_t> fExoE;
+ std::vector<UShort_t> fExoE_CrystalNbr;
+ std::vector<ULong64_t> fExoE_TS;
- std::vector<UShort_t> fEXO_HG; // Same as Energy but with High Gain (for Low Energy events)
- std::vector<UShort_t> fEXO_HG_CrystalNbr;
- std::vector<ULong64_t> fEXO_HG_TS;
+ std::vector<UShort_t> fExoHG; // Same as Energy but with High Gain (for Low Energy events)
+ std::vector<UShort_t> fExoHG_CrystalNbr;
+ std::vector<ULong64_t> fExoHG_TS;
- std::vector<UShort_t> fEXO_TDC; // Internal TDC of EXOGAM
- std::vector<UShort_t> fEXO_TDC_CrystalNbr;
- std::vector<ULong64_t> fEXO_TDC_TS;
+ std::vector<UShort_t> fExoTDC; // Internal TDC of EXOGAM
+ std::vector<UShort_t> fExoTDC_CrystalNbr;
+ std::vector<ULong64_t> fExoTDC_TS;
- std::vector<UShort_t> fEXO_Outer;
- std::vector<UShort_t> fEXO_Outer_SubCrystalNbr;
+ std::vector<UShort_t> fExoOuter;
+ std::vector<UShort_t> fExoOuter_SubCrystalNbr;
- std::vector<UShort_t> fEXO_BGO;
- std::vector<UShort_t> fEXO_BGO_CrystalNbr;
+ std::vector<UShort_t> fExoBGO;
+ std::vector<UShort_t> fExoBGO_CrystalNbr;
- std::vector<UShort_t> fEXO_CSI;
- std::vector<UShort_t> fEXO_CSI_CrystalNbr;
+ std::vector<UShort_t> fExoCsI;
+ std::vector<UShort_t> fExoCsI_CrystalNbr;
public:
TExogamData();
@@ -58,49 +58,49 @@ class TExogamData : public TObject {
void Dump() const;
///////////////////// SETTERS ////////////////////////
- inline void SetInner6MV(const UShort_t& Energy, const UShort_t& DetNumb, const ULong64_t& TimeStamp) {
- fEXO_E.push_back(Energy);
- fEXO_E_CrystalNbr.push_back(DetNumb);
- fEXO_E_TS.push_back(TimeStamp);
+ inline void SetExo(const UShort_t& Energy, const UShort_t& DetNumb, const ULong64_t& TimeStamp) {
+ fExoE.push_back(Energy);
+ fExoE_CrystalNbr.push_back(DetNumb);
+ fExoE_TS.push_back(TimeStamp);
}
- inline void SetInner20MV(const UShort_t& Energy, const UShort_t& DetNumb, const ULong64_t& TimeStamp) {
- fEXO_HG.push_back(Energy);
- fEXO_HG_CrystalNbr.push_back(DetNumb);
- fEXO_HG_TS.push_back(TimeStamp);
+ inline void SetExoHG(const UShort_t& Energy, const UShort_t& DetNumb, const ULong64_t& TimeStamp) {
+ fExoHG.push_back(Energy);
+ fExoHG_CrystalNbr.push_back(DetNumb);
+ fExoHG_TS.push_back(TimeStamp);
}
- inline void SetDeltaTV(const UShort_t& Time, const UShort_t& DetNumb, const ULong64_t& TimeStamp) {
- fEXO_TDC.push_back(Time);
- fEXO_TDC_CrystalNbr.push_back(DetNumb);
- fEXO_TDC_TS.push_back(TimeStamp);
+ inline void SetExoDelta(const UShort_t& Time, const UShort_t& DetNumb, const ULong64_t& TimeStamp) {
+ fExoTDC.push_back(Time);
+ fExoTDC_CrystalNbr.push_back(DetNumb);
+ fExoTDC_TS.push_back(TimeStamp);
}
- inline void SetOutersV(const UShort_t& Energy, const UShort_t& OutersNumb) {
- fEXO_Outer.push_back(Energy);
- fEXO_Outer_SubCrystalNbr.push_back(OutersNumb);
+ inline void SetExoOuter(const UShort_t& Energy, const UShort_t& OutersNumb) {
+ fExoOuter.push_back(Energy);
+ fExoOuter_SubCrystalNbr.push_back(OutersNumb);
}
- inline void SetBGOV(const UShort_t& Energy, const UShort_t& BGONumb) {
- fEXO_BGO.push_back(Energy);
- fEXO_BGO_CrystalNbr.push_back(BGONumb);
+ inline void SetExoBGO(const UShort_t& Energy, const UShort_t& BGONumb) {
+ fExoBGO.push_back(Energy);
+ fExoBGO_CrystalNbr.push_back(BGONumb);
}
- inline void SetCSIV(const UShort_t& Energy, const UShort_t& CSINumb) {
- fEXO_CSI.push_back(Energy);
- fEXO_CSI_CrystalNbr.push_back(CSINumb);
+ inline void SetExoCsI(const UShort_t& Energy, const UShort_t& CsINumb) {
+ fExoCsI.push_back(Energy);
+ fExoCsI_CrystalNbr.push_back(CsINumb);
}
///////////////////// GETTERS ////////////////////////
- inline UShort_t GetEXO_E(const UShort_t& i)const { return fEXO_E[i]; }
- inline UShort_t GetEXO_E_CrystalNbr(const UShort_t& i)const { return fEXO_E_CrystalNbr[i]; }
- inline ULong64_t GetEXO_E_TS(const UShort_t& i)const { return fEXO_E_TS[i]; }
- inline UShort_t GetEXO_HG(const UShort_t& i)const { return fEXO_HG[i]; }
- inline UShort_t GetEXO_HG_CrystalNbr(const UShort_t& i)const { return fEXO_HG_CrystalNbr[i]; }
- inline ULong64_t GetEXO_HG_TS(const UShort_t& i)const { return fEXO_HG_TS[i]; }
- inline UShort_t GetEXO_TDC(const UShort_t& i)const { return fEXO_TDC[i]; }
- inline UShort_t GetEXO_TDC_CrystalNbr(const UShort_t& i)const { return fEXO_TDC_CrystalNbr[i]; }
- inline ULong64_t GetEXO_TDC_TS(const UShort_t& i)const { return fEXO_TDC_TS[i]; }
- inline UShort_t GetEXO_Outer(const UShort_t& i)const { return fEXO_Outer[i]; }
- inline UShort_t GetEXO_Outer_SubCrystalNbr(const UShort_t& i)const { return fEXO_Outer_SubCrystalNbr[i]; }
- inline UShort_t GetEXO_BGO(const UShort_t& i)const { return fEXO_BGO[i]; }
- inline UShort_t GetEXO_BGO_CrystalNbr(const UShort_t& i)const { return fEXO_BGO_CrystalNbr[i]; }
- inline UShort_t GetEXO_CSI(const UShort_t& i)const { return fEXO_CSI[i]; }
- inline UShort_t GetEXO_CSI_CrystalNbr(const UShort_t& i)const { return fEXO_CSI_CrystalNbr[i]; }
+ inline UShort_t GetExoE(UShort_t& i) { return fExoE[i]; }
+ inline UShort_t GetExoE_CrystalNbr(UShort_t& i) { return fExoE_CrystalNbr[i]; }
+ inline ULong64_t GetExoE_TS(UShort_t& i) { return fExoE_TS[i]; }
+ inline UShort_t GetExoHG(UShort_t& i) { return fExoHG[i]; }
+ inline UShort_t GetExoHG_CrystalNbr(UShort_t& i) { return fExoHG_CrystalNbr[i]; }
+ inline ULong64_t GetExoHG_TS(UShort_t& i) { return fExoHG_TS[i]; }
+ inline UShort_t GetExoTDC(UShort_t& i) { return fExoTDC[i]; }
+ inline UShort_t GetExoTDC_CrystalNbr(UShort_t& i) { return fExoTDC_CrystalNbr[i]; }
+ inline ULong64_t GetExoTDC_TS(UShort_t& i) { return fExoTDC_TS[i]; }
+ inline UShort_t GetExoOuter(UShort_t& i) { return fExoOuter[i]; }
+ inline UShort_t GetExoOuter_SubCrystalNbr(UShort_t& i) { return fExoOuter_SubCrystalNbr[i]; }
+ inline UShort_t GetExoBGO(UShort_t& i) { return fExoBGO[i]; }
+ inline UShort_t GetExoBGO_CrystalNbr(UShort_t& i) { return fExoBGO_CrystalNbr[i]; }
+ inline UShort_t GetExoCsI(UShort_t& i) { return fExoCsI[i]; }
+ inline UShort_t GetExoCsI_CrystalNbr(UShort_t& i) { return fExoCsI_CrystalNbr[i]; }
ClassDef(TExogamData, 1) // ExogamData structure
};
diff --git a/NPLib/Detectors/Exogam/TExogamPhysics.cxx b/NPLib/Detectors/Exogam/TExogamPhysics.cxx
index 58ea53a90b9a04c7c64b08164514dbe759c27777..c9d907539d557cc63d4b3cb14cc47ebfffe18ced 100644
--- a/NPLib/Detectors/Exogam/TExogamPhysics.cxx
+++ b/NPLib/Detectors/Exogam/TExogamPhysics.cxx
@@ -20,56 +20,48 @@
*****************************************************************************/
#include "TExogamPhysics.h"
using namespace EXOGAM_LOCAL;
-
+
// STL
-#include <sstream>
-#include <iostream>
#include <cmath>
+#include <iostream>
+#include <sstream>
#include <stdlib.h>
// NPL
-#include "RootInput.h"
#include "NPDetectorFactory.h"
-#include "RootOutput.h"
-#include "NPVDetector.h"
#include "NPOptionManager.h"
+#include "NPVDetector.h"
+#include "RootInput.h"
+#include "RootOutput.h"
// ROOT
#include "TChain.h"
-
///////////////////////////////////////////////////////////////////////////
ClassImp(TExogamPhysics)
-///////////////////////////////////////////////////////////////////////////
-TExogamPhysics::TExogamPhysics()
-{
- EventMultiplicity = 0 ;
- ECC_Multiplicity = 0 ;
- GOCCE_Multiplicity = 0 ;
- NumberOfHitClover = 0 ;
- NumberOfHitCristal = 0 ;
- m_Spectra = NULL;
- NumberOfClover=0;
-
- PreTreatedData = new TExogamData ;
- EventData = new TExogamData ;
- EventPhysics = this ;
- NumberOfClover = 0 ;
- CloverMult = 0 ;
-
+ ///////////////////////////////////////////////////////////////////////////
+ TExogamPhysics::TExogamPhysics() {
+ EventMultiplicity = 0;
+ ECC_Multiplicity = 0;
+ GOCCE_Multiplicity = 0;
+ NumberOfHitClover = 0;
+ NumberOfHitCristal = 0;
+ m_Spectra = NULL;
+ NumberOfClover = 0;
+
+ PreTreatedData = new TExogamData;
+ EventData = new TExogamData;
+ EventPhysics = this;
+ NumberOfClover = 0;
+ CloverMult = 0;
}
-
+
///////////////////////////////////////////////////////////////////////////
-void TExogamPhysics::BuildSimplePhysicalEvent()
-{
-
- BuildPhysicalEvent();
-}
+void TExogamPhysics::BuildSimplePhysicalEvent() { BuildPhysicalEvent(); }
///////////////////////////////////////////////////////////////////////////
-void TExogamPhysics::PreTreat()
-{
+void TExogamPhysics::PreTreat() {
/*ClearPreTreatedData();
- //E
+ //E
for(unsigned int i = 0 ; i < EventData -> GetECCEMult(); i++) {
UShort_t cristal_E = 10000 ; UShort_t cristal_T = 2000;
@@ -77,83 +69,88 @@ void TExogamPhysics::PreTreat()
{
int clover = EventData -> GetECCEClover(i);
int cristal = EventData -> GetECCECristal(i);
-
- if(EventData -> GetECCEEnergy(i) < 3000) cristal_E = CalibrationManager::getInstance()-> ApplyCalibration("EXOGAM/Cl"+ NPL::itoa(clover)+"_Cr"+ NPL::itoa(cristal)+"_Elow", EventData -> GetECCEEnergy(i));
- else cristal_E = CalibrationManager::getInstance()-> ApplyCalibration("EXOGAM/Cl"+ NPL::itoa(clover)+"_Cr"+ NPL::itoa(cristal)+"_Ehigh", EventData -> GetECCEEnergy(i));
-
+ if(EventData -> GetECCEEnergy(i) < 3000) cristal_E = CalibrationManager::getInstance()->
+ ApplyCalibration("EXOGAM/Cl"+ NPL::itoa(clover)+"_Cr"+ NPL::itoa(cristal)+"_Elow", EventData -> GetECCEEnergy(i));
+ else cristal_E = CalibrationManager::getInstance()->
+ ApplyCalibration("EXOGAM/Cl"+ NPL::itoa(clover)+"_Cr"+ NPL::itoa(cristal)+"_Ehigh", EventData -> GetECCEEnergy(i));
+
+
if(cristal_E > Threshold_ECC)
- {
-
- PreTreatedData->SetECCEClover ( clover ) ;
- PreTreatedData->SetECCECristal( cristal ) ;
- PreTreatedData->SetECCEEnergy ( cristal_E ) ;
-
- bool checkT = false;
- for(unsigned int k = 0; k < EventData -> GetECCTMult(); k++){
- if(clover == EventData -> GetECCTClover(k) && cristal == EventData -> GetECCTCristal(k)){
- // cout << EventData -> GetECCTTime(k) << endl;
-
- if(EventData -> GetECCTTime(k) < 16383) cristal_T = CalibrationManager::getInstance()-> ApplyCalibration("EXOGAM/Cl"+ NPL::itoa(clover)+"_Cr"+ NPL::itoa(cristal)+"_T", EventData -> GetECCTTime(k));
- else cristal_T = 2500;
-
- //if(cristal_T >5000 && cristal_T !=25000 ) cout << "PreTreat " << cristal_T << " " << EventData -> GetECCTTime(k) << " " << clover << " " << cristal << " " << EventData->GetECCTMult() << endl;
-
- checkT=true;
- PreTreatedData->SetECCTClover (clover ) ;
- PreTreatedData->SetECCTCristal( cristal ) ;
- PreTreatedData->SetECCTTime ( cristal_T ) ;
-
- ECC_Multiplicity ++;
- GOCCE_Multiplicity++;
- }
-
- }
-
- if(!checkT) {
- PreTreatedData->SetECCTClover (clover ) ;
- PreTreatedData->SetECCTCristal( cristal ) ;
- PreTreatedData->SetECCTTime ( -1000 ) ;
- }
-
-
- }
+ {
+
+ PreTreatedData->SetECCEClover ( clover ) ;
+ PreTreatedData->SetECCECristal( cristal ) ;
+ PreTreatedData->SetECCEEnergy ( cristal_E ) ;
+
+ bool checkT = false;
+ for(unsigned int k = 0; k < EventData -> GetECCTMult(); k++){
+ if(clover == EventData -> GetECCTClover(k) && cristal == EventData -> GetECCTCristal(k)){
+ // cout << EventData -> GetECCTTime(k) << endl;
+
+ if(EventData -> GetECCTTime(k) < 16383) cristal_T = CalibrationManager::getInstance()->
+ ApplyCalibration("EXOGAM/Cl"+ NPL::itoa(clover)+"_Cr"+ NPL::itoa(cristal)+"_T", EventData -> GetECCTTime(k)); else
+ cristal_T = 2500;
+
+ //if(cristal_T >5000 && cristal_T !=25000 ) cout << "PreTreat " << cristal_T << " " << EventData ->
+ GetECCTTime(k) << " " << clover << " " << cristal << " " << EventData->GetECCTMult() << endl;
+
+ checkT=true;
+ PreTreatedData->SetECCTClover (clover ) ;
+ PreTreatedData->SetECCTCristal( cristal ) ;
+ PreTreatedData->SetECCTTime ( cristal_T ) ;
+
+ ECC_Multiplicity ++;
+ GOCCE_Multiplicity++;
+ }
+
+ }
+
+ if(!checkT) {
+ PreTreatedData->SetECCTClover (clover ) ;
+ PreTreatedData->SetECCTCristal( cristal ) ;
+ PreTreatedData->SetECCTTime ( -1000 ) ;
+ }
+
+
+ }
}
}
//cout << PreTreatedData-> GetECCTMult() << " " << PreTreatedData-> GetECCEMult() << endl;
-
+
//GOCCE
- //E
-
+ //E
+
for(unsigned int i = 0 ; i < EventData -> GetGOCCEEMult(); i++) {
- UShort_t segment_E = 25000;
+ UShort_t segment_E = 25000;
//if(IsValidChannel)
{
int clover = EventData -> GetGOCCEEClover(i);
int cristal = EventData -> GetGOCCEECristal(i);
int segment = EventData -> GetGOCCEESegment(i);
-
+
if(EventData -> GetGOCCEEEnergy(i) > RawThreshold_GOCCE)
- {
- segment_E = CalibrationManager::getInstance()->ApplyCalibration("EXOGAM/Cl"+ NPL::itoa(clover)+"_Cr"+ NPL::itoa(cristal)+"_Seg"+ NPL::itoa(segment)+"_E", EventData -> GetGOCCEEEnergy(i));
-
- if(segment_E > Threshold_GOCCE)
- {
- PreTreatedData->SetGOCCEEClover ( clover ) ;
- PreTreatedData->SetGOCCEECristal( cristal ) ;
- PreTreatedData->SetGOCCEESegment( segment ) ;
- PreTreatedData->SetGOCCEEEnergy ( segment_E ) ;
-
- }
- }
+ {
+ segment_E = CalibrationManager::getInstance()->ApplyCalibration("EXOGAM/Cl"+ NPL::itoa(clover)+"_Cr"+
+ NPL::itoa(cristal)+"_Seg"+ NPL::itoa(segment)+"_E", EventData -> GetGOCCEEEnergy(i));
+
+ if(segment_E > Threshold_GOCCE)
+ {
+ PreTreatedData->SetGOCCEEClover ( clover ) ;
+ PreTreatedData->SetGOCCEECristal( cristal ) ;
+ PreTreatedData->SetGOCCEESegment( segment ) ;
+ PreTreatedData->SetGOCCEEEnergy ( segment_E ) ;
+
+ }
+ }
else
- {
-
- }
+ {
+
+ }
}
}
@@ -162,110 +159,111 @@ void TExogamPhysics::PreTreat()
*/
}
///////////////////////////////////////////////////////////////////////////
-
-void TExogamPhysics::BuildPhysicalEvent()
-{
+
+void TExogamPhysics::BuildPhysicalEvent() {
/*PreTreat();
-
- if(PreTreatedData -> GetECCEMult() != PreTreatedData -> GetECCTMult()) cout << PreTreatedData -> GetECCEMult() << " " << PreTreatedData -> GetECCTMult() << endl;
-
+ if(PreTreatedData -> GetECCEMult() != PreTreatedData -> GetECCTMult()) cout << PreTreatedData -> GetECCEMult() << " "
+ << PreTreatedData -> GetECCTMult() << endl;
+
+
for(unsigned int i = 0 ; i < PreTreatedData -> GetECCEMult(); i++) {
-
+
// cout << i << " " << cristal_E << endl;
- // if(PreTreatedData->GetECCTTime(i) > 0)
+ // if(PreTreatedData->GetECCTTime(i) > 0)
{
ECC_E.push_back(PreTreatedData->GetECCEEnergy(i));
ECC_T.push_back(PreTreatedData->GetECCTTime(i));
ECC_CloverNumber.push_back(PreTreatedData->GetECCEClover(i));
ECC_CristalNumber.push_back(PreTreatedData->GetECCECristal(i));
- // cout << "BuildPhys " << PreTreatedData->GetECCEClover(i) << " " << PreTreatedData->GetECCECristal(i)<< " " << PreTreatedData->GetECCTTime(i) << " " << endl;
+ // cout << "BuildPhys " << PreTreatedData->GetECCEClover(i) << " " << PreTreatedData->GetECCECristal(i)<< " " <<
+ PreTreatedData->GetECCTTime(i) << " " << endl;
}
}
-
+
for(unsigned int j = 0 ; j < PreTreatedData -> GetGOCCEEMult(); j++) {
GOCCE_E.push_back(PreTreatedData->GetGOCCEEEnergy(j));
GOCCE_CloverNumber.push_back(PreTreatedData->GetGOCCEEClover(j));
GOCCE_CristalNumber.push_back(PreTreatedData->GetGOCCEECristal(j));
GOCCE_SegmentNumber.push_back(PreTreatedData->GetGOCCEESegment(j));
}
-
+
//int NumberOfHitClover = 0;
-
+
int DetectorID = -1;
-
+
for( unsigned short i = 0 ; i < PreTreatedData->GetECCEMult() ; i++ )
- {
+ {
// cout << PreTreatedData->GetECCEClover(i) << endl;
if( PreTreatedData->GetECCEClover(i) != DetectorID)
- {
- if(i==0)
- {
- NumberOfHitClover++;
- }
- else if(PreTreatedData->GetECCEClover(i)!= PreTreatedData->GetECCEClover(i-1) )
- {
- NumberOfHitClover++;
- }
- }
- if(NumberOfHitClover == 4) break;
+ {
+ if(i==0)
+ {
+ NumberOfHitClover++;
+ }
+ else if(PreTreatedData->GetECCEClover(i)!= PreTreatedData->GetECCEClover(i-1) )
+ {
+ NumberOfHitClover++;
+ }
+ }
+ if(NumberOfHitClover == 4) break;
//clover_mult -> Fill(NumberOfHitClover);
-
+
}
-
+
//cout << "NumberOfHitClover " << NumberOfHitClover << endl;
-
+
map<int, vector<int> > MapCristal;
- map<int, vector<int> > MapSegment;
-
+ map<int, vector<int> > MapSegment;
+
map<int, vector<int> > :: iterator it; // iterator used with MapCristal
map<int, vector<int> > :: iterator at; // iterator used with MapSegment
vector<int> PositionOfCristal_Buffer_ECC;
vector<int> PositionOfSegment_Buffer_GOCCE;
-
+
//Fill map Cristal
for(int clo = 0; clo < NumberOfClover; clo++)
{
for(unsigned int k = 0; k < ECC_CloverNumber.size(); k++)
- {
- if(ECC_CloverNumber.at(k) == clo) // && ECC_CristalNumber.at(k)== cri )
- PositionOfCristal_Buffer_ECC.push_back(k);
- }
+ {
+ if(ECC_CloverNumber.at(k) == clo) // && ECC_CristalNumber.at(k)== cri )
+ PositionOfCristal_Buffer_ECC.push_back(k);
+ }
if(PositionOfCristal_Buffer_ECC.size() != 0) MapCristal[clo] = PositionOfCristal_Buffer_ECC;
-
+
PositionOfCristal_Buffer_ECC.clear();
-
+
}
-
+
//Fill map Segment
for(int clo = 0; clo < NumberOfClover; clo++)
{
for(int cri = 0; cri < 4 ; cri++)
- {
- // for(int seg = 0; seg < 4 ; seg++)
- {
- for(unsigned int m = 0; m < GOCCE_CloverNumber.size(); m++)
- {
- if(GOCCE_CloverNumber.at(m) == clo && GOCCE_CristalNumber.at(m) == cri)// && GOCCE_SegmentNumber.at(m) == seg)
- {
- // PositionOfSegment_Buffer_GOCCE.push_back(4*clo+cri);
- PositionOfSegment_Buffer_GOCCE.push_back(m);
- }
- }
- }
- if(PositionOfSegment_Buffer_GOCCE.size() != 0) MapSegment[4*clo+cri] = PositionOfSegment_Buffer_GOCCE;
-
- PositionOfSegment_Buffer_GOCCE.clear();
- }
+ {
+ // for(int seg = 0; seg < 4 ; seg++)
+ {
+ for(unsigned int m = 0; m < GOCCE_CloverNumber.size(); m++)
+ {
+ if(GOCCE_CloverNumber.at(m) == clo && GOCCE_CristalNumber.at(m) == cri)// && GOCCE_SegmentNumber.at(m) == seg)
+ {
+ // PositionOfSegment_Buffer_GOCCE.push_back(4*clo+cri);
+ PositionOfSegment_Buffer_GOCCE.push_back(m);
+ }
}
+ }
+ if(PositionOfSegment_Buffer_GOCCE.size() != 0) MapSegment[4*clo+cri] = PositionOfSegment_Buffer_GOCCE;
- // Treatment
+ PositionOfSegment_Buffer_GOCCE.clear();
+ }
+ }
+
+ // Treatment
for(int clo = 0; clo < NumberOfClover ; clo++)
{
double E = 0; double T = 0;
@@ -276,202 +274,195 @@ void TExogamPhysics::BuildPhysicalEvent()
int Emax = 0, Emin = 1000000;
int Tmin = 0, Tmax = 0;
- //ADD-BACK
- it = MapCristal.find(clo);
-
+ //ADD-BACK
+ it = MapCristal.find(clo);
+
int cristal_cond = 0;
-
- if(it != MapCristal.end())
- {
- vector<int> PositionOfCristal = it -> second;
-
- mult_cristal = PositionOfCristal.size();
- //if(mult_cristal!=0) cristal_mult -> Fill(mult_cristal);
-
- // ADD-BACK
- //cout << "boucle" << endl;
-
- for(unsigned int k = 0; k < PositionOfCristal.size(); k++)
- {
- int indice = PositionOfCristal.at(k);
-
- cristal_cond += ECC_CristalNumber.at(indice);
- // cout << ECC_CristalNumber.at(k) << " " ECC_E.at(k) << endl;
-
- if(mult_cristal < 3)
- {
- E+= ECC_E.at(indice);
-
- if(ECC_E.at(indice) < Emin) {
- cristal_Emin = ECC_CristalNumber.at(indice);
- Emin = ECC_E.at(indice);
- Tmin = ECC_T.at(indice);
- }
-
- if(ECC_E.at(indice) > Emax) {
- cristal_Emax = ECC_CristalNumber.at(indice);
- Emax = ECC_E.at(indice);
- Tmax = ECC_T.at(indice);
- }
- }
-
- else // case of multiplicity = 3 or 4
- {
- E = -1; cristal_Emax = -1; cristal_Emin = -1; Tmax = -1; Tmin = -1;
- }
-
- // cout << ECC_E.at(indice) << " " << Emax << " " << cristal_Emax << " " << Emin << " " << cristal_Emin << endl;
-
- }
-
- if( (mult_cristal==1) || (mult_cristal ==2 && cristal_cond %2 == 1) )
- {
- // cout << cristal_cond << endl;
-
- //cristal = cristal_Emax; T = Tmax;
- //cout << Emax << " " << cristal_Emax << " " << Emin << " " << cristal_Emin << endl;
-
- if(E > 500) { cristal = cristal_Emax; T = Tmax; }
- else { cristal = cristal_Emin; T = Tmin; }
-
-
- // DOPPLER CORRECTION
-
- at = MapSegment.find(4*clo+cristal);
- segment = -1;
-
- if(at != MapSegment.end())
- {
- vector<int> PositionOfSegment = at -> second; // position of segment k in the vector
-
- int segment_max = -1, E_temp = -1;
-
- for(unsigned int m = 0; m < PositionOfSegment.size(); m++) // loop on hit segments of cristal cri of clover clo
- {
- int indice = PositionOfSegment.at(m);
-
- if(GOCCE_E.at(indice) > 0 && GOCCE_CristalNumber.at(indice) == cristal)
- {
- if( GOCCE_E.at(indice) > E_temp )
- {
- segment_max = GOCCE_SegmentNumber.at(indice) ;
- E_temp = GOCCE_E.at(indice);
- }
- }
- }
- segment = segment_max;
- }
-
- }
-
-
- if(E > 0 && cristal != -1 && segment != -1)
- {
- TotalEnergy_lab.push_back(E);
- Time.push_back(T);
- CloverNumber.push_back(clo);
- CristalNumber.push_back(cristal);
- SegmentNumber.push_back(segment);
-
- double theta = GetSegmentAngleTheta(clo, cristal, segment);
-
- Theta.push_back(theta);
-
- double doppler_E = DopplerCorrection(E, theta);
- DopplerCorrectedEnergy.push_back(doppler_E);
-
- // cout << E << " " << clo << " " << cristal << " " << segment << " " << theta << " " << doppler_E << endl;
-
- }
-
- } // end of condition over CristalMap
+
+ if(it != MapCristal.end())
+ {
+ vector<int> PositionOfCristal = it -> second;
+
+ mult_cristal = PositionOfCristal.size();
+ //if(mult_cristal!=0) cristal_mult -> Fill(mult_cristal);
+
+ // ADD-BACK
+ //cout << "boucle" << endl;
+
+ for(unsigned int k = 0; k < PositionOfCristal.size(); k++)
+ {
+ int indice = PositionOfCristal.at(k);
+
+ cristal_cond += ECC_CristalNumber.at(indice);
+ // cout << ECC_CristalNumber.at(k) << " " ECC_E.at(k) << endl;
+
+ if(mult_cristal < 3)
+ {
+ E+= ECC_E.at(indice);
+
+ if(ECC_E.at(indice) < Emin) {
+ cristal_Emin = ECC_CristalNumber.at(indice);
+ Emin = ECC_E.at(indice);
+ Tmin = ECC_T.at(indice);
+ }
+
+ if(ECC_E.at(indice) > Emax) {
+ cristal_Emax = ECC_CristalNumber.at(indice);
+ Emax = ECC_E.at(indice);
+ Tmax = ECC_T.at(indice);
+ }
+ }
+
+ else // case of multiplicity = 3 or 4
+ {
+ E = -1; cristal_Emax = -1; cristal_Emin = -1; Tmax = -1; Tmin = -1;
+ }
+
+ // cout << ECC_E.at(indice) << " " << Emax << " " << cristal_Emax << " " << Emin << " " << cristal_Emin << endl;
+
+ }
+
+ if( (mult_cristal==1) || (mult_cristal ==2 && cristal_cond %2 == 1) )
+ {
+ // cout << cristal_cond << endl;
+
+ //cristal = cristal_Emax; T = Tmax;
+ //cout << Emax << " " << cristal_Emax << " " << Emin << " " << cristal_Emin << endl;
+
+ if(E > 500) { cristal = cristal_Emax; T = Tmax; }
+ else { cristal = cristal_Emin; T = Tmin; }
+
+
+ // DOPPLER CORRECTION
+
+ at = MapSegment.find(4*clo+cristal);
+ segment = -1;
+
+ if(at != MapSegment.end())
+ {
+ vector<int> PositionOfSegment = at -> second; // position of segment k in the vector
+
+ int segment_max = -1, E_temp = -1;
+
+ for(unsigned int m = 0; m < PositionOfSegment.size(); m++) // loop on hit segments of cristal cri of
+ clover clo
+ {
+ int indice = PositionOfSegment.at(m);
+
+ if(GOCCE_E.at(indice) > 0 && GOCCE_CristalNumber.at(indice) == cristal)
+ {
+ if( GOCCE_E.at(indice) > E_temp )
+ {
+ segment_max = GOCCE_SegmentNumber.at(indice) ;
+ E_temp = GOCCE_E.at(indice);
+ }
+ }
+ }
+ segment = segment_max;
+ }
+
+ }
+
+
+ if(E > 0 && cristal != -1 && segment != -1)
+ {
+ TotalEnergy_lab.push_back(E);
+ Time.push_back(T);
+ CloverNumber.push_back(clo);
+ CristalNumber.push_back(cristal);
+ SegmentNumber.push_back(segment);
+
+ double theta = GetSegmentAngleTheta(clo, cristal, segment);
+
+ Theta.push_back(theta);
+
+ double doppler_E = DopplerCorrection(E, theta);
+ DopplerCorrectedEnergy.push_back(doppler_E);
+
+ // cout << E << " " << clo << " " << cristal << " " << segment << " " << theta << " " << doppler_E << endl;
+
+ }
+
+ } // end of condition over CristalMap
} // loop over NumberOfClover
CloverMult = GetClover_Mult();
- //cout << "Exogam fine" << endl;
+ //cout << "Exogam fine" << endl;
*/
-}
-
+}
-double TExogamPhysics::DopplerCorrection(double E, double Theta)
-{
- double Pi = 3.141592654 ;
+double TExogamPhysics::DopplerCorrection(double E, double Theta) {
+ double Pi = 3.141592654;
TString filename = "configs/beta.txt";
ifstream file;
- //cout << filename << endl;
+ // cout << filename << endl;
file.open(filename);
- if(!file) cout << filename << " was not opened" << endl;
+ if (!file)
+ cout << filename << " was not opened" << endl;
double E_corr = 0;
- double beta = 0.;
- file>>beta;
- double gamma = 1./ sqrt(1-beta*beta);
+ double beta = 0.;
+ file >> beta;
+ double gamma = 1. / sqrt(1 - beta * beta);
- E_corr = gamma * E * ( 1. - beta * cos(Theta*Pi/180.));
-
- return(E_corr);
+ E_corr = gamma * E * (1. - beta * cos(Theta * Pi / 180.));
+ return (E_corr);
}
-
///////////////////////////////////////////////////////////////////////////
-
-void TExogamPhysics::Clear()
-{
- EventMultiplicity = 0;
- ECC_Multiplicity = 0;
+void TExogamPhysics::Clear() {
+ EventMultiplicity = 0;
+ ECC_Multiplicity = 0;
GOCCE_Multiplicity = 0;
- NumberOfHitClover = 0;
- NumberOfHitCristal = 0;
+ NumberOfHitClover = 0;
+ NumberOfHitCristal = 0;
- ECC_CloverNumber .clear() ;
- ECC_CristalNumber .clear() ;
- GOCCE_CloverNumber .clear() ;
- GOCCE_CristalNumber .clear() ;
- GOCCE_SegmentNumber .clear() ;
+ ECC_CloverNumber.clear();
+ ECC_CristalNumber.clear();
+ GOCCE_CloverNumber.clear();
+ GOCCE_CristalNumber.clear();
+ GOCCE_SegmentNumber.clear();
// ECC
- ECC_E.clear() ;
+ ECC_E.clear();
ECC_T.clear();
// GOCCE
- GOCCE_E.clear() ;
-
- CristalNumber.clear() ;
- SegmentNumber.clear() ;
- CloverNumber .clear() ;
-
- TotalEnergy_lab .clear() ;
- Time .clear() ;
- DopplerCorrectedEnergy.clear() ;
- Position .clear() ;
- Theta .clear() ;
-
+ GOCCE_E.clear();
+ CristalNumber.clear();
+ SegmentNumber.clear();
+ CloverNumber.clear();
+ TotalEnergy_lab.clear();
+ Time.clear();
+ DopplerCorrectedEnergy.clear();
+ Position.clear();
+ Theta.clear();
}
///////////////////////////////////////////////////////////////////////////
-//// Innherited from VDetector Class ////
-
+//// Innherited from VDetector Class ////
+
// Read stream at ConfigFile to pick-up parameters of detector (Position,...) using Token
-void TExogamPhysics::ReadConfiguration(NPL::InputParser parser){
+void TExogamPhysics::ReadConfiguration(NPL::InputParser parser) {
vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("Exogam");
- if(NPOptionManager::getInstance()->GetVerboseLevel())
- cout << "//// " << blocks.size() << " detectors found " << endl;
+ if (NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << "//// " << blocks.size() << " detectors found " << endl;
vector<string> token = {"ANGLE_FILE"};
- for(unsigned int i = 0 ; i < blocks.size() ; i++){
- if(blocks[i]->HasTokenList(token)){
+ for (unsigned int i = 0; i < blocks.size(); i++) {
+ if (blocks[i]->HasTokenList(token)) {
string AngleFile = blocks[i]->GetString("ANGLE_FILE");
AddClover(AngleFile);
}
- else{
+ else {
cout << "ERROR: check your input file formatting " << endl;
exit(1);
}
@@ -479,122 +470,111 @@ void TExogamPhysics::ReadConfiguration(NPL::InputParser parser){
}
///////////////////////////////////////////////////////////////////////////
-void TExogamPhysics::InitSpectra(){
- m_Spectra = new TExogamSpectra(NumberOfClover);
-}
+void TExogamPhysics::InitSpectra() { m_Spectra = new TExogamSpectra(NumberOfClover); }
///////////////////////////////////////////////////////////////////////////
-void TExogamPhysics::FillSpectra(){
- m_Spectra -> FillRawSpectra(EventData);
- m_Spectra -> FillPreTreatedSpectra(PreTreatedData);
- m_Spectra -> FillPhysicsSpectra(EventPhysics);
+void TExogamPhysics::FillSpectra() {
+ m_Spectra->FillRawSpectra(EventData);
+ m_Spectra->FillPreTreatedSpectra(PreTreatedData);
+ m_Spectra->FillPhysicsSpectra(EventPhysics);
}
///////////////////////////////////////////////////////////////////////////
-void TExogamPhysics::CheckSpectra(){
- m_Spectra->CheckSpectra();
-}
+void TExogamPhysics::CheckSpectra() { m_Spectra->CheckSpectra(); }
///////////////////////////////////////////////////////////////////////////
-void TExogamPhysics::ClearSpectra(){
+void TExogamPhysics::ClearSpectra() {
// To be done
}
///////////////////////////////////////////////////////////////////////////
-map< string , TH1*> TExogamPhysics::GetSpectra() {
- if(m_Spectra)
+map<string, TH1*> TExogamPhysics::GetSpectra() {
+ if (m_Spectra)
return m_Spectra->GetMapHisto();
- else{
- map< string , TH1*> empty;
+ else {
+ map<string, TH1*> empty;
return empty;
}
-}
+}
//////////////////////////////////////////////////////////////////////////
-void TExogamPhysics::AddClover(string AngleFile)
-{
+void TExogamPhysics::AddClover(string AngleFile) {
ifstream file;
// TString filename = Form("posBaptiste/angles_exogam_clover%d.txt",NumberOfClover);
// TString filename = Form("posz42_simu50mm/angles_exogam_clover%d.txt",NumberOfClover);
// TString filename = Form("posz42_exp_stat_demiring/angles_exogam_clover%d.txt",NumberOfClover);
-
+
string path = "configs/";
TString filename = path + AngleFile;
-
+
cout << filename << endl;
file.open(filename);
- if(!file) cout << filename << " was not opened" << endl;
+ if (!file)
+ cout << filename << " was not opened" << endl;
+
+ vector<double> Angles;
+ vector<vector<double>> Segment_angles;
+ vector<vector<vector<double>>> Cristal_angles;
- vector <double> Angles;
- vector < vector <double> > Segment_angles;
- vector < vector < vector <double> > > Cristal_angles;
-
Cristal_angles.clear();
- double angle; string buffer;
+ double angle;
+ string buffer;
- for(int i = 0; i < 4; i++)
- {
- Segment_angles.clear();
-
- for(int j = 0; j < 4; j++)
- {
- Angles.clear();
-
- for(int k = 0; k < 2; k++)
- {
- file >> buffer >> angle;
-
- Angles.push_back(angle); // Theta (k = 0) Phi (k = 1)
-
- //cout << angle << endl;
- if(Angles.size()==2)
- cout << "Clover " << NumberOfClover << ": Theta=" << Angles[0] << " Phi=" << Angles[1]<< endl;
-
- }
-
- Segment_angles.push_back(Angles);
- }
-
- Cristal_angles.push_back(Segment_angles);
+ for (int i = 0; i < 4; i++) {
+ Segment_angles.clear();
+
+ for (int j = 0; j < 4; j++) {
+ Angles.clear();
+
+ for (int k = 0; k < 2; k++) {
+ file >> buffer >> angle;
+
+ Angles.push_back(angle); // Theta (k = 0) Phi (k = 1)
+
+ // cout << angle << endl;
+ if (Angles.size() == 2)
+ cout << "Clover " << NumberOfClover << ": Theta=" << Angles[0] << " Phi=" << Angles[1] << endl;
+ }
+
+ Segment_angles.push_back(Angles);
}
+ Cristal_angles.push_back(Segment_angles);
+ }
+
Clover_Angles_Theta_Phi.push_back(Cristal_angles);
file.close();
NumberOfClover++;
-
-}
-
+}
// Add Parameter to the CalibrationManger
-void TExogamPhysics::AddParameterToCalibrationManager()
-{
-
+void TExogamPhysics::AddParameterToCalibrationManager() {
+
CalibrationManager* Cal = CalibrationManager::getInstance();
-
- for(int i = 0 ; i < NumberOfClover ; i++)
- {
- for( int j = 0 ; j < 4 ; j++)
- {
- Cal->AddParameter("EXOGAM", "Cl"+ NPL::itoa(i)+"_Cr"+ NPL::itoa(j)+"_Elow" ,"EXOGAM_Cl"+ NPL::itoa(i)+"_Cr"+ NPL::itoa(j)+"_Elow");
- Cal->AddParameter("EXOGAM", "Cl"+ NPL::itoa(i)+"_Cr"+ NPL::itoa(j)+"_Ehigh","EXOGAM_Cl"+ NPL::itoa(i)+"_Cr"+ NPL::itoa(j)+"_Ehigh");
- Cal->AddParameter("EXOGAM", "Cl"+ NPL::itoa(i)+"_Cr"+ NPL::itoa(j)+"_T","EXOGAM_Cl"+ NPL::itoa(i)+"_Cr"+ NPL::itoa(j)+"_T") ;
-
- for( int k = 0 ; k < 4 ; k++)
- {
- Cal->AddParameter("EXOGAM", "Cl"+ NPL::itoa(i)+"_Cr"+ NPL::itoa(j)+"_Seg"+ NPL::itoa(k)+"_E","EXOGAM_Cl"+ NPL::itoa(i)+"_Cr"+ NPL::itoa(j)+"_Seg"+ NPL::itoa(k)+"_E") ;
- }
- }
+
+ for (int i = 0; i < NumberOfClover; i++) {
+ for (int j = 0; j < 4; j++) {
+ Cal->AddParameter("EXOGAM", "Cl" + NPL::itoa(i) + "_Cr" + NPL::itoa(j) + "_Elow",
+ "EXOGAM_Cl" + NPL::itoa(i) + "_Cr" + NPL::itoa(j) + "_Elow");
+ Cal->AddParameter("EXOGAM", "Cl" + NPL::itoa(i) + "_Cr" + NPL::itoa(j) + "_Ehigh",
+ "EXOGAM_Cl" + NPL::itoa(i) + "_Cr" + NPL::itoa(j) + "_Ehigh");
+ Cal->AddParameter("EXOGAM", "Cl" + NPL::itoa(i) + "_Cr" + NPL::itoa(j) + "_T",
+ "EXOGAM_Cl" + NPL::itoa(i) + "_Cr" + NPL::itoa(j) + "_T");
+
+ for (int k = 0; k < 4; k++) {
+ Cal->AddParameter("EXOGAM", "Cl" + NPL::itoa(i) + "_Cr" + NPL::itoa(j) + "_Seg" + NPL::itoa(k) + "_E",
+ "EXOGAM_Cl" + NPL::itoa(i) + "_Cr" + NPL::itoa(j) + "_Seg" + NPL::itoa(k) + "_E");
+ }
}
+ }
}
-
// Activated associated Branches and link it to the private member DetectorData address
// In this method mother Branches (Detector) AND daughter leaf (fDetector_parameter) have to be activated
-void TExogamPhysics::InitializeRootInputRaw()
-{
- TChain* inputChain = RootInput::getInstance()->GetChain() ;
- inputChain->SetBranchStatus( "EXOGAM" , true ) ;
- inputChain->SetBranchStatus( "fEXO_*" , true ) ;
- inputChain->SetBranchAddress( "EXOGAM" , &EventData ) ;
+void TExogamPhysics::InitializeRootInputRaw() {
+ TChain* inputChain = RootInput::getInstance()->GetChain();
+ inputChain->SetBranchStatus("EXOGAM", true);
+ inputChain->SetBranchStatus("fEXO_*", true);
+ inputChain->SetBranchAddress("EXOGAM", &EventData);
/*
TList* outputList = RootOutput::getInstance()->GetList();
@@ -610,82 +590,74 @@ void TExogamPhysics::InitializeRootInputRaw()
// In this method mother Branches (Detector) AND daughter leaf (parameter) have to be activated
void TExogamPhysics::InitializeRootInputPhysics() {
TChain* inputChain = RootInput::getInstance()->GetChain();
- inputChain->SetBranchStatus( "EventMultiplicty" , true );
- inputChain->SetBranchStatus( "ECC_Multiplicity" , true );
- inputChain->SetBranchStatus( "GOCCE_Multiplicity" , true );
- inputChain->SetBranchStatus( "ECC_CloverNumber" , true );
- inputChain->SetBranchStatus( "ECC_CristalNumber" , true );
- inputChain->SetBranchStatus( "GOCCE_CloverNumber" , true );
- inputChain->SetBranchStatus( "GOCCE_CristalNumber" , true );
- inputChain->SetBranchStatus( "GOCCE_SegmentNumber" , true );
- inputChain->SetBranchStatus( "ECC_E" , true );
- inputChain->SetBranchStatus( "ECC_T" , true );
- inputChain->SetBranchStatus( "GOCCE_E" , true );
- inputChain->SetBranchStatus( "CristalNumber" , true );
- inputChain->SetBranchStatus( "SegmentNumber" , true );
- inputChain->SetBranchStatus( "CloverNumber" , true );
- inputChain->SetBranchStatus( "CloverMult" , true );
- inputChain->SetBranchStatus( "TotalEnergy_lab" , true );
- inputChain->SetBranchStatus( "Time" , true );
- inputChain->SetBranchStatus( "DopplerCorrectedEnergy" , true );
- inputChain->SetBranchStatus( "Position" , true );
- inputChain->SetBranchStatus( "Theta" , true );
- inputChain->SetBranchAddress( "EXOGAM" , &EventPhysics );
-
+ inputChain->SetBranchStatus("EventMultiplicty", true);
+ inputChain->SetBranchStatus("ECC_Multiplicity", true);
+ inputChain->SetBranchStatus("GOCCE_Multiplicity", true);
+ inputChain->SetBranchStatus("ECC_CloverNumber", true);
+ inputChain->SetBranchStatus("ECC_CristalNumber", true);
+ inputChain->SetBranchStatus("GOCCE_CloverNumber", true);
+ inputChain->SetBranchStatus("GOCCE_CristalNumber", true);
+ inputChain->SetBranchStatus("GOCCE_SegmentNumber", true);
+ inputChain->SetBranchStatus("ECC_E", true);
+ inputChain->SetBranchStatus("ECC_T", true);
+ inputChain->SetBranchStatus("GOCCE_E", true);
+ inputChain->SetBranchStatus("CristalNumber", true);
+ inputChain->SetBranchStatus("SegmentNumber", true);
+ inputChain->SetBranchStatus("CloverNumber", true);
+ inputChain->SetBranchStatus("CloverMult", true);
+ inputChain->SetBranchStatus("TotalEnergy_lab", true);
+ inputChain->SetBranchStatus("Time", true);
+ inputChain->SetBranchStatus("DopplerCorrectedEnergy", true);
+ inputChain->SetBranchStatus("Position", true);
+ inputChain->SetBranchStatus("Theta", true);
+ inputChain->SetBranchAddress("EXOGAM", &EventPhysics);
}
/////////////////////////////////////////////////////////////////////
// Create associated branches and associated private member DetectorPhysics address
-void TExogamPhysics::InitializeRootOutput()
-{
- TTree* outputTree = RootOutput::getInstance()->GetTree() ;
- outputTree->Branch( "EXOGAM" , "TExogamPhysics" , &EventPhysics ) ;
+void TExogamPhysics::InitializeRootOutput() {
+ TTree* outputTree = RootOutput::getInstance()->GetTree();
+ outputTree->Branch("EXOGAM", "TExogamPhysics", &EventPhysics);
// control histograms if needed
- /*
+ /*
TList* outputList = RootOutput::getInstance()->GetList();
controle = new TH1F("controle","histo de controle",20,0,20);
outputList->Add(controle);
*/
-
}
-
///////////////////////////////////////////////////////////////////////////
-namespace EXOGAM_LOCAL
-{
+namespace EXOGAM_LOCAL {
// tranform an integer to a string
- string itoa(int value)
- {
+ string itoa(int value) {
std::ostringstream o;
-
+
if (!(o << value))
- return "" ;
-
+ return "";
+
return o.str();
}
-}
-
- /////////////////////////////
+} // namespace EXOGAM_LOCAL
+
+/////////////////////////////
////////////////////////////////////////////////////////////////////////////////
// Construct Method to be pass to the DetectorFactory //
////////////////////////////////////////////////////////////////////////////////
-NPL::VDetector* TExogamPhysics::Construct(){
- return (NPL::VDetector*) new TExogamPhysics();
-}
+NPL::VDetector* TExogamPhysics::Construct() { return (NPL::VDetector*)new TExogamPhysics(); }
////////////////////////////////////////////////////////////////////////////////
// Registering the construct method to the factory //
////////////////////////////////////////////////////////////////////////////////
-extern "C"{
-class proxy_exogam{
- public:
- proxy_exogam(){
- NPL::DetectorFactory::getInstance()->AddToken("Exogam","Exogam");
- NPL::DetectorFactory::getInstance()->AddDetector("Exogam",TExogamPhysics::Construct);
- }
+extern "C" {
+class proxy_exogam {
+ public:
+ proxy_exogam() {
+ NPL::DetectorFactory::getInstance()->AddToken("Exogam", "Exogam");
+ NPL::DetectorFactory::getInstance()->AddDetector("Exogam", TExogamPhysics::Construct);
+ }
};
proxy_exogam p;
-}
\ No newline at end of file
+}
diff --git a/NPLib/Detectors/PISTA/TPISTAPhysics.cxx b/NPLib/Detectors/PISTA/TPISTAPhysics.cxx
index 95084a5c2c5857868c975c092a5b8ef12c287429..aed56416356ee6c5e6e50ff3c2f6966437401a84 100644
--- a/NPLib/Detectors/PISTA/TPISTAPhysics.cxx
+++ b/NPLib/Detectors/PISTA/TPISTAPhysics.cxx
@@ -110,7 +110,7 @@ void TPISTAPhysics::AddDetector(TVector3 A, TVector3 B, TVector3 C, TVector3 D){
//Strip_1_1 = A + u*(StripPitchX / 2.) + v*(StripPitchY / 2.);
Strip_1_1 = A + u*deltaX + u*(ContractedStripPitchX / 2.) + v*(NumberOfStripsY*StripPitchY - StripPitchY / 2.);
- Strip_1_1 = Strip_1_1 + 2.9*u + 2.0*v;
+ //Strip_1_1 = Strip_1_1 + 2.9*u + 2.0*v;
TVector3 StripPos;
for(int i=0; i<NumberOfStripsX; i++){
diff --git a/NPLib/Detectors/Vendeta/TVendetaPhysics.cxx b/NPLib/Detectors/Vendeta/TVendetaPhysics.cxx
index 4353ebb26b8b9a41c0f1bea0ae8d7041cc034e96..94ce084b7d9a76edb64c4af33fb5bb49fa442316 100644
--- a/NPLib/Detectors/Vendeta/TVendetaPhysics.cxx
+++ b/NPLib/Detectors/Vendeta/TVendetaPhysics.cxx
@@ -1,24 +1,24 @@
/*****************************************************************************
- * Copyright (C) 2009-2022 this file is part of the NPTool Project *
- * *
- * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
- * For the list of contributors see $NPTOOL/Licence/Contributors *
- *****************************************************************************/
+ * Copyright (C) 2009-2022 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
/*****************************************************************************
- * Original Author: Pierre Morfouace contact address: pierre.morfouace@cea.fr *
- * *
- * Creation Date : February 2022 *
- * Last update : *
- *---------------------------------------------------------------------------*
- * Decription: *
- * This class hold Vendeta Treated data *
- * *
- *---------------------------------------------------------------------------*
- * Comment: *
- * *
- * *
- *****************************************************************************/
+ * Original Author: Pierre Morfouace contact address: pierre.morfouace@cea.fr *
+ * *
+ * Creation Date : February 2022 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * This class hold Vendeta Treated data *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
#include "TVendetaPhysics.h"
@@ -42,39 +42,39 @@ using namespace std;
ClassImp(TVendetaPhysics)
- ///////////////////////////////////////////////////////////////////////////
+ ///////////////////////////////////////////////////////////////////////////
TVendetaPhysics::TVendetaPhysics()
- : m_EventData(new TVendetaData),
- m_PreTreatedData(new TVendetaData),
- m_EventPhysics(this),
- m_Spectra(0),
- m_E_RAW_Threshold(0), // adc channels
- m_E_Threshold(0), // MeV
- m_AnodeNumber(1),
- m_NumberOfDetectors(0) {
- }
+ : m_EventData(new TVendetaData),
+ m_PreTreatedData(new TVendetaData),
+ m_EventPhysics(this),
+ m_Spectra(0),
+ m_E_RAW_Threshold(0), // adc channels
+ m_E_Threshold(0), // MeV
+ m_AnodeNumber(1),
+ m_NumberOfDetectors(0) {
+ }
///////////////////////////////////////////////////////////////////////////
/// A usefull method to bundle all operation to add a detector
void TVendetaPhysics::AddDetector(TVector3 ){
- // In That simple case nothing is done
- // Typically for more complex detector one would calculate the relevant
- // positions (stripped silicon) or angles (gamma array)
- m_NumberOfDetectors++;
+ // In That simple case nothing is done
+ // Typically for more complex detector one would calculate the relevant
+ // positions (stripped silicon) or angles (gamma array)
+ m_NumberOfDetectors++;
}
///////////////////////////////////////////////////////////////////////////
void TVendetaPhysics::AddDetector(double R, double Theta, double Phi){
- // Compute the TVector3 corresponding
- TVector3 Pos(R*sin(Theta)*cos(Phi),R*sin(Theta)*sin(Phi),R*cos(Theta));
- // Call the cartesian method
- AddDetector(Pos);
- m_DetectorPosition.push_back(Pos);
+ // Compute the TVector3 corresponding
+ TVector3 Pos(R*sin(Theta)*cos(Phi),R*sin(Theta)*sin(Phi),R*cos(Theta));
+ // Call the cartesian method
+ AddDetector(Pos);
+ m_DetectorPosition.push_back(Pos);
}
///////////////////////////////////////////////////////////////////////////
void TVendetaPhysics::BuildSimplePhysicalEvent() {
- BuildPhysicalEvent();
+ BuildPhysicalEvent();
}
@@ -82,91 +82,91 @@ void TVendetaPhysics::BuildSimplePhysicalEvent() {
///////////////////////////////////////////////////////////////////////////
void TVendetaPhysics::BuildPhysicalEvent() {
- // Treat Event, only if Fission Chamber has triggered
- if(m_AnodeNumber==-1)
- return;
+ // Treat Event, only if Fission Chamber has triggered
+ if(m_AnodeNumber==-1)
+ return;
- // apply thresholds and calibration
- //if(m_AnodeNumber==0)
- // return;
+ // apply thresholds and calibration
+ //if(m_AnodeNumber==0)
+ // return;
- PreTreat();
+ PreTreat();
- // match energy and time together
- unsigned int mysizeLGE = m_PreTreatedData->GetLGMultEnergy();
- unsigned int mysizeHGE = m_PreTreatedData->GetHGMultEnergy();
+ // match energy and time together
+ unsigned int mysizeLGE = m_PreTreatedData->GetLGMultEnergy();
+ unsigned int mysizeHGE = m_PreTreatedData->GetHGMultEnergy();
- for (UShort_t e = 0; e < mysizeLGE ; e++) {
+ for (UShort_t e = 0; e < mysizeLGE ; e++) {
- LG_DetectorNumber.push_back(m_PreTreatedData->GetLGDetectorNbr(e));
- LG_Q1.push_back(m_PreTreatedData->GetLGQ1(e));
- LG_Q2.push_back(m_PreTreatedData->GetLGQ2(e));
- LG_Time.push_back(m_PreTreatedData->GetLGTime(e));
- LG_Qmax.push_back(m_PreTreatedData->GetLGQmax(e));
+ LG_DetectorNumber.push_back(m_PreTreatedData->GetLGDetectorNbr(e));
+ LG_Q1.push_back(m_PreTreatedData->GetLGQ1(e));
+ LG_Q2.push_back(m_PreTreatedData->GetLGQ2(e));
+ LG_Time.push_back(m_PreTreatedData->GetLGTime(e));
+ LG_Qmax.push_back(m_PreTreatedData->GetLGQmax(e));
- }
+ }
+
+ for (UShort_t e = 0; e < mysizeHGE ; e++) {
+ HG_DetectorNumber.push_back(m_PreTreatedData->GetHGDetectorNbr(e));
+ HG_Q1.push_back(m_PreTreatedData->GetHGQ1(e));
+ HG_Q2.push_back(m_PreTreatedData->GetHGQ2(e));
+ HG_Time.push_back(m_PreTreatedData->GetHGTime(e));
+ HG_Qmax.push_back(m_PreTreatedData->GetHGQmax(e));
+ }
- for (UShort_t e = 0; e < mysizeHGE ; e++) {
- HG_DetectorNumber.push_back(m_PreTreatedData->GetHGDetectorNbr(e));
- HG_Q1.push_back(m_PreTreatedData->GetHGQ1(e));
- HG_Q2.push_back(m_PreTreatedData->GetHGQ2(e));
- HG_Time.push_back(m_PreTreatedData->GetHGTime(e));
- HG_Qmax.push_back(m_PreTreatedData->GetHGQmax(e));
- }
-
- m_AnodeNumber=-1;
+ m_AnodeNumber=-1;
}
///////////////////////////////////////////////////////////////////////////
void TVendetaPhysics::PreTreat() {
- // This method typically applies thresholds and calibrations
- // Might test for disabled channels for more complex detector
-
- // clear pre-treated object
- ClearPreTreatedData();
- // instantiate CalibrationManager
- static CalibrationManager* Cal = CalibrationManager::getInstance();
- unsigned int mysizeLG = m_EventData->GetLGMultEnergy();
- unsigned int mysizeHG = m_EventData->GetHGMultEnergy();
-
- // LG pretreat
- for (UShort_t i = 0; i < mysizeLG ; ++i){
- int det = m_EventData->GetLGDetectorNbr(i);
- double Qmax = m_EventData->GetLGQmax(i);
- double TimeOffset=0;
- TimeOffset = Cal->GetValue("Vendeta/DET"+NPL::itoa(det)+"_LG_ANODE"+NPL::itoa(m_AnodeNumber)+"_TIMEOFFSET",0);
- if(m_AnodeNumber==0){
- // Apply calibration from Anode 6 in case of Pulser
- TimeOffset = Cal->GetValue("Vendeta/DET"+NPL::itoa(det)+"_LG_ANODE"+NPL::itoa(6)+"_TIMEOFFSET",0);
- }
- double Time = m_EventData->GetLGTime(i) + TimeOffset;
- m_PreTreatedData->SetLGDetectorNbr(det);
- m_PreTreatedData->SetLGQ1(m_EventData->GetLGQ1(i));
- m_PreTreatedData->SetLGQ2(m_EventData->GetLGQ2(i));
- m_PreTreatedData->SetLGTime(Time);
- m_PreTreatedData->SetLGQmax(Qmax);
+ // This method typically applies thresholds and calibrations
+ // Might test for disabled channels for more complex detector
+
+ // clear pre-treated object
+ ClearPreTreatedData();
+ // instantiate CalibrationManager
+ static CalibrationManager* Cal = CalibrationManager::getInstance();
+ unsigned int mysizeLG = m_EventData->GetLGMultEnergy();
+ unsigned int mysizeHG = m_EventData->GetHGMultEnergy();
+
+ // LG pretreat
+ for (UShort_t i = 0; i < mysizeLG ; ++i){
+ int det = m_EventData->GetLGDetectorNbr(i);
+ double Qmax = m_EventData->GetLGQmax(i);
+ double TimeOffset=0;
+ TimeOffset = Cal->GetValue("Vendeta/DET"+NPL::itoa(det)+"_LG_ANODE"+NPL::itoa(m_AnodeNumber)+"_TIMEOFFSET",0);
+ if(m_AnodeNumber==0){
+ // Apply calibration from Anode 6 in case of Pulser
+ TimeOffset = Cal->GetValue("Vendeta/DET"+NPL::itoa(det)+"_LG_ANODE"+NPL::itoa(6)+"_TIMEOFFSET",0);
}
+ double Time = m_EventData->GetLGTime(i) + TimeOffset;
+ m_PreTreatedData->SetLGDetectorNbr(det);
+ m_PreTreatedData->SetLGQ1(m_EventData->GetLGQ1(i));
+ m_PreTreatedData->SetLGQ2(m_EventData->GetLGQ2(i));
+ m_PreTreatedData->SetLGTime(Time);
+ m_PreTreatedData->SetLGQmax(Qmax);
+ }
- // HG pretreat
- for (UShort_t i = 0; i < mysizeHG ; ++i){
- int det = m_EventData->GetHGDetectorNbr(i);
- double Qmax = m_EventData->GetHGQmax(i);
- double TimeOffset=0;
- TimeOffset = Cal->GetValue("Vendeta/DET"+NPL::itoa(det)+"_HG_ANODE"+NPL::itoa(m_AnodeNumber)+"_TIMEOFFSET",0);
- if(m_AnodeNumber==0){
- // Apply calibration from Anode 6 in case of Pulser
- TimeOffset = Cal->GetValue("Vendeta/DET"+NPL::itoa(det)+"_HG_ANODE"+NPL::itoa(6)+"_TIMEOFFSET",0);
- }
-
- double Time = m_EventData->GetHGTime(i) + TimeOffset;
- m_PreTreatedData->SetHGDetectorNbr(det);
- m_PreTreatedData->SetHGQ1(m_EventData->GetHGQ1(i));
- m_PreTreatedData->SetHGQ2(m_EventData->GetHGQ2(i));
- m_PreTreatedData->SetHGTime(Time);
- m_PreTreatedData->SetHGQmax(Qmax);
+ // HG pretreat
+ for (UShort_t i = 0; i < mysizeHG ; ++i){
+ int det = m_EventData->GetHGDetectorNbr(i);
+ double Qmax = m_EventData->GetHGQmax(i);
+ double TimeOffset=0;
+ TimeOffset = Cal->GetValue("Vendeta/DET"+NPL::itoa(det)+"_HG_ANODE"+NPL::itoa(m_AnodeNumber)+"_TIMEOFFSET",0);
+ if(m_AnodeNumber==0){
+ // Apply calibration from Anode 6 in case of Pulser
+ TimeOffset = Cal->GetValue("Vendeta/DET"+NPL::itoa(det)+"_HG_ANODE"+NPL::itoa(6)+"_TIMEOFFSET",0);
}
+ double Time = m_EventData->GetHGTime(i) + TimeOffset;
+ m_PreTreatedData->SetHGDetectorNbr(det);
+ m_PreTreatedData->SetHGQ1(m_EventData->GetHGQ1(i));
+ m_PreTreatedData->SetHGQ2(m_EventData->GetHGQ2(i));
+ m_PreTreatedData->SetHGTime(Time);
+ m_PreTreatedData->SetHGQmax(Qmax);
+ }
+
}
diff --git a/NPLib/Physics/NPReaction.cxx b/NPLib/Physics/NPReaction.cxx
index 40f07bbecfd13e2260b743ac66edb63d0def98af..b548dc168d24c69a70cafa4742fadda43c165b97 100644
--- a/NPLib/Physics/NPReaction.cxx
+++ b/NPLib/Physics/NPReaction.cxx
@@ -173,8 +173,10 @@ Reaction::~Reaction() {}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
bool Reaction::CheckKinematic() {
double theta = fThetaCM;
- if (m1 > m2)
+ /*if (m1 > m2){
theta = M_PI - fThetaCM;
+ fThetaCM = M_PI - fThetaCM;
+ }*/
fEnergyImpulsionCM_3 = TLorentzVector(pCM_3 * sin(theta), 0, pCM_3 * cos(theta), ECM_3);
fEnergyImpulsionCM_4 = fTotalEnergyImpulsionCM - fEnergyImpulsionCM_3;
@@ -259,8 +261,10 @@ void Reaction::KineRelativistic(double& ThetaLab3, double& KineticEnergyLab3, do
// case of inverse kinematics
double theta = fThetaCM;
- if (m1 > m2)
+ /*if (m1 > m2){
theta = M_PI - fThetaCM;
+ //fThetaCM = M_PI - fThetaCM;
+ }*/
fEnergyImpulsionCM_3 = TLorentzVector(pCM_3 * sin(theta), 0, pCM_3 * cos(theta), ECM_3);
fEnergyImpulsionCM_4 = fTotalEnergyImpulsionCM - fEnergyImpulsionCM_3;
@@ -646,7 +650,7 @@ TGraph* Reaction::GetKinematicLine3(double AngleStep_CM) {
vector<double> vy;
double theta3, E3, theta4, E4;
- for (double angle = 0; angle < 360; angle += AngleStep_CM) {
+ for (double angle = 0; angle < 181; angle += AngleStep_CM) {
SetThetaCM(angle * deg);
KineRelativistic(theta3, E3, theta4, E4);
fParticle3.SetKineticEnergy(E3);
@@ -668,7 +672,7 @@ TGraph* Reaction::GetKinematicLine4(double AngleStep_CM) {
vector<double> vy;
double theta3, E3, theta4, E4;
- for (double angle = 0; angle < 360; angle += AngleStep_CM) {
+ for (double angle = 0; angle < 181; angle += AngleStep_CM) {
SetThetaCM(angle * deg);
KineRelativistic(theta3, E3, theta4, E4);
fParticle4.SetKineticEnergy(E4);
@@ -689,7 +693,7 @@ TGraph* Reaction::GetTheta3VsTheta4(double AngleStep_CM) {
vector<double> vy;
double theta3, E3, theta4, E4;
- for (double angle = 0; angle < 360; angle += AngleStep_CM) {
+ for (double angle = 0; angle < 181; angle += AngleStep_CM) {
SetThetaCM(angle * deg);
KineRelativistic(theta3, E3, theta4, E4);
@@ -708,7 +712,7 @@ TGraph* Reaction::GetBrhoLine3(double AngleStep_CM) {
double theta3, E3, theta4, E4;
double Brho;
- for (double angle = 0; angle < 360; angle += AngleStep_CM) {
+ for (double angle = 0; angle < 181; angle += AngleStep_CM) {
SetThetaCM(angle * deg);
KineRelativistic(theta3, E3, theta4, E4);
fParticle3.SetKineticEnergy(E3);
@@ -728,7 +732,7 @@ TGraph* Reaction::GetThetaLabVersusThetaCM(double AngleStep_CM) {
vector<double> vy;
double theta3, E3, theta4, E4;
- for (double angle = 0; angle < 360; angle += AngleStep_CM) {
+ for (double angle = 0; angle < 181; angle += AngleStep_CM) {
SetThetaCM(angle * deg);
KineRelativistic(theta3, E3, theta4, E4);
@@ -739,6 +743,46 @@ TGraph* Reaction::GetThetaLabVersusThetaCM(double AngleStep_CM) {
fAngleLine = new TGraph(vx.size(), &vx[0], &vy[0]);
return (fAngleLine);
}
+////////////////////////////////////////////////////////////////////////////////////////////
+TGraph* Reaction::GetJacobian(double AngleStep_CM){
+ vector<double> vx;
+ vector<double> vy;
+ double theta3, E3, theta4, E4;
+
+ double E3_CM = GetE_CM_3();
+ double P3_CM = GetP_CM_3();
+ double B = P3_CM/E3_CM;
+ double Mass1 = fParticle1.Mass();
+ double Mass2 = fParticle2.Mass();
+ double BeamEnergy = GetBeamEnergy();
+ double beta = sqrt(BeamEnergy*BeamEnergy+2*Mass1*BeamEnergy)/(BeamEnergy+Mass1+Mass2);
+ double r = beta/B;
+ double gamma = 1./sqrt(1-beta*beta);
+
+ double Jacobian_num;
+ double Jacobian_denum;
+ double Jacobian;
+ for (double angle = 0; angle < 181; angle += AngleStep_CM) {
+ SetThetaCM(angle * deg);
+ KineRelativistic(theta3, E3, theta4, E4);
+
+ Jacobian_num = abs(gamma*(1.+r*cos(fThetaCM)));
+
+ Jacobian_denum = pow(gamma,2)*pow(r+cos(fThetaCM),2) + pow(sin(fThetaCM),2);
+ Jacobian_denum = pow(Jacobian_denum,3./2.);
+
+ Jacobian = Jacobian_num / Jacobian_denum;
+
+ vx.push_back(fThetaCM / deg);
+ vy.push_back(Jacobian);
+ }
+
+ TGraph* gJaco = new TGraph(vx.size(), &vx[0], &vy[0]);
+
+ return gJaco;
+
+}
+
////////////////////////////////////////////////////////////////////////////////////////////
TGraph* Reaction::GetELabVersusThetaCM(double AngleStep_CM) {
@@ -746,7 +790,7 @@ TGraph* Reaction::GetELabVersusThetaCM(double AngleStep_CM) {
vector<double> vy;
double theta3, E3, theta4, E4;
- for (double angle = 0; angle < 360; angle += AngleStep_CM) {
+ for (double angle = 0; angle < 180; angle += AngleStep_CM) {
SetThetaCM(angle * deg);
KineRelativistic(theta3, E3, theta4, E4);
diff --git a/NPLib/Physics/NPReaction.h b/NPLib/Physics/NPReaction.h
index 0f61637780f971794e2a42a0cfd1ca40a81f167f..455aadcc6b2408a346c6acc7ee4dd16d05e64419 100644
--- a/NPLib/Physics/NPReaction.h
+++ b/NPLib/Physics/NPReaction.h
@@ -264,6 +264,7 @@ namespace NPL {
TGraph* GetThetaLabVersusThetaCM(double AngleStep_CM = 1);
TGraph* GetELabVersusThetaCM(double AngleStep_CM = 1);
TGraph* GetTheta3VsTheta4(double AngleStep_CM = 1);
+ TGraph* GetJacobian(double AngleStep_CM = 1);
void PrintKinematic();
double GetP_CM_1() { return pCM_1; }
double GetP_CM_2() { return pCM_2; }
diff --git a/NPLib/Utility/npcalibration b/NPLib/Utility/npcalibration
deleted file mode 100755
index 400a50a9cc0087aa1c26dc69a6ee84eed9eec2e2..0000000000000000000000000000000000000000
Binary files a/NPLib/Utility/npcalibration and /dev/null differ
diff --git a/NPSimulation/Detectors/PISTA/PISTA.cc b/NPSimulation/Detectors/PISTA/PISTA.cc
index c5f7d993eb85cdb3330bee0e7f0510b47acb8d28..0407fa38f37aa10910780654f92cabf565afa00b 100644
--- a/NPSimulation/Detectors/PISTA/PISTA.cc
+++ b/NPSimulation/Detectors/PISTA/PISTA.cc
@@ -69,7 +69,7 @@ namespace PISTA_NS{
const double TrapezoidHeight = 57.7*mm;
const double TrapezoidLength = 1*cm;
const double FirstStageThickness = 100*um;
- const double SecondStageThickness = 1*mm;
+ const double SecondStageThickness = 1.0*mm;
const double DistanceBetweenSi = 4*mm;
const double FirstStageNbrOfStrips = 91;
const double SecondStageNbrOfStrips = 57;
@@ -319,34 +319,89 @@ void PISTA::ReadSensitive(const G4Event* ){
// DE
DSSDScorers::PS_Rectangle* FirstStageScorer= (DSSDScorers::PS_Rectangle*) m_FirstStageScorer->GetPrimitive(0);
- unsigned int sizeDE = FirstStageScorer->GetLengthMult();
- for(unsigned int i = 0 ; i < sizeDE ; i++){
- double Energy = RandGauss::shoot(FirstStageScorer->GetEnergyLength(i), DE_ResoEnergy);
- if(Energy>EnergyThreshold){
+ unsigned int sizeDEFront = FirstStageScorer->GetWidthMult();
+ for(unsigned int i = 0 ; i < sizeDEFront ; i++){
+ double EnergyFront = RandGauss::shoot(FirstStageScorer->GetEnergyWidth(i), DE_ResoEnergy);
+ if(EnergyFront>EnergyThreshold){
+ double Time = RandGauss::shoot(FirstStageScorer->GetTimeWidth(i), ResoTime);
+ int DetNbr = FirstStageScorer->GetDetectorWidth(i);
+ int StripFront = 92-FirstStageScorer->GetStripWidth(i);
+ m_Event->SetPISTA_DE_DetectorNbr(DetNbr);
+ m_Event->SetPISTA_DE_StripNbr(StripFront);
+ m_Event->SetPISTA_DE_StripEnergy(EnergyFront);
+ m_Event->SetPISTA_DE_StripTime(Time);
+ }
+ }
+
+ unsigned int sizeDEBack = FirstStageScorer->GetLengthMult();
+ for(unsigned int i = 0 ; i < sizeDEBack ; i++){
+ double EnergyBack = RandGauss::shoot(FirstStageScorer->GetEnergyLength(i), DE_ResoEnergy);
+ if(EnergyBack>EnergyThreshold){
double Time = RandGauss::shoot(FirstStageScorer->GetTimeLength(i), ResoTime);
int DetNbr = FirstStageScorer->GetDetectorLength(i);
- int StripFront = 92-FirstStageScorer->GetStripWidth(i);
- m_Event->SetPISTA_DE(DetNbr, StripFront, Energy, Energy, Time, Time);
m_Event->SetPISTA_DE_BackDetector(DetNbr);
+ m_Event->SetPISTA_DE_BackEnergy(EnergyBack);
+ m_Event->SetPISTA_DE_BackTime(Time);
+
}
}
+
+ /*for(unsigned int i = 0 ; i < sizeDEFront ; i++){
+ double EnergyFront = RandGauss::shoot(FirstStageScorer->GetEnergyWidth(i), DE_ResoEnergy);
+ double EnergyBack = RandGauss::shoot(FirstStageScorer->GetEnergyLength(i), DE_ResoEnergy);
+ if(EnergyFront>EnergyThreshold){
+ double Time = RandGauss::shoot(FirstStageScorer->GetTimeLength(i), ResoTime);
+ int DetNbr = FirstStageScorer->GetDetectorWidth(i);
+ int StripFront = 92-FirstStageScorer->GetStripWidth(i);
+ m_Event->SetPISTA_DE(DetNbr, StripFront, EnergyFront, EnergyBack, Time, Time);
+ m_Event->SetPISTA_DE_BackDetector(DetNbr);
+ }
+ }*/
FirstStageScorer->clear();
///////////
// E
DSSDScorers::PS_Rectangle* SecondStageScorer= (DSSDScorers::PS_Rectangle*) m_SecondStageScorer->GetPrimitive(0);
- unsigned int sizeE = SecondStageScorer->GetLengthMult();
- for(unsigned int i = 0 ; i < sizeE ; i++){
- double Energy = RandGauss::shoot(SecondStageScorer->GetEnergyLength(i), E_ResoEnergy);
- if(Energy>EnergyThreshold){
+ unsigned int sizeEFront = SecondStageScorer->GetLengthMult();
+ for(unsigned int i = 0 ; i < sizeEFront ; i++){
+ double EnergyFront = RandGauss::shoot(SecondStageScorer->GetEnergyLength(i), E_ResoEnergy);
+ if(EnergyFront>EnergyThreshold){
double Time = RandGauss::shoot(SecondStageScorer->GetTimeLength(i), ResoTime);
int DetNbr = SecondStageScorer->GetDetectorLength(i);
int StripFront = SecondStageScorer->GetStripLength(i);
- m_Event->SetPISTA_E(DetNbr, StripFront, Energy, Energy, Time, Time);
+ m_Event->SetPISTA_E_DetectorNbr(DetNbr);
+ m_Event->SetPISTA_E_StripNbr(StripFront);
+ m_Event->SetPISTA_E_StripEnergy(EnergyFront);
+ m_Event->SetPISTA_E_StripTime(Time);
+ }
+ }
+
+ unsigned int sizeEBack = SecondStageScorer->GetWidthMult();
+ for(unsigned int i = 0 ; i < sizeEBack ; i++){
+ double EnergyBack = RandGauss::shoot(SecondStageScorer->GetEnergyWidth(i), E_ResoEnergy);
+ if(EnergyBack>EnergyThreshold){
+ double Time = RandGauss::shoot(SecondStageScorer->GetTimeWidth(i), ResoTime);
+ int DetNbr = SecondStageScorer->GetDetectorWidth(i);
m_Event->SetPISTA_E_BackDetector(DetNbr);
+ m_Event->SetPISTA_E_BackEnergy(EnergyBack);
+ m_Event->SetPISTA_E_BackTime(Time);
}
}
+
+
+
+ /*for(unsigned int i = 0 ; i < sizeEFront ; i++){
+ double EnergyFront = RandGauss::shoot(SecondStageScorer->GetEnergyLength(i), E_ResoEnergy);
+ double EnergyBack = RandGauss::shoot(SecondStageScorer->GetEnergyWidth(i), E_ResoEnergy);
+ if(EnergyFront>EnergyThreshold){
+ double Time = RandGauss::shoot(SecondStageScorer->GetTimeLength(i), ResoTime);
+ int DetNbr = SecondStageScorer->GetDetectorLength(i);
+ int StripFront = SecondStageScorer->GetStripLength(i);
+ m_Event->SetPISTA_E(DetNbr, StripFront, EnergyFront, EnergyBack, Time, Time);
+ m_Event->SetPISTA_E_BackDetector(DetNbr);
+ }
+ }*/
SecondStageScorer->clear();
diff --git a/NPSimulation/Detectors/Vendeta/Vendeta.cc b/NPSimulation/Detectors/Vendeta/Vendeta.cc
index af97c3a19098f87a52aab1b17a654449d8f5bf2c..c0ca4ccc14e5b5ea390c6e4d309dcd66933ebc0b 100644
--- a/NPSimulation/Detectors/Vendeta/Vendeta.cc
+++ b/NPSimulation/Detectors/Vendeta/Vendeta.cc
@@ -61,7 +61,7 @@ using namespace CLHEP;
namespace Vendeta_NS{
// Energy and time Resolution
const double EnergyThreshold = 0.1*MeV;
- const double ResoTime = 0.6*ns ;
+ const double ResoTime = 0.8/2.35*ns ;
const double ResoEnergy = 0.1*MeV ;
const double Thickness = 51.*mm ;
const double Radius = 127./2*mm ;
diff --git a/NPSimulation/Scorers/PlasticBar.cc b/NPSimulation/Scorers/PlasticBar.cc
index 9b80e7d832a2afcee216d6fd1f76a1c0bc03127d..1ea8e9ad5aff2a936240a92efb354976e35cbb6a 100644
--- a/NPSimulation/Scorers/PlasticBar.cc
+++ b/NPSimulation/Scorers/PlasticBar.cc
@@ -27,11 +27,10 @@ using namespace PlasticBar;
using namespace std;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-unsigned int
-PlasticBarData::CalculateIndex(const vector<unsigned int>& level) {
+unsigned int PlasticBarData::CalculateIndex(const vector<unsigned int>& level) {
- unsigned int size = level.size();
- unsigned int result = 0;
+ unsigned int size = level.size();
+ unsigned int result = 0;
unsigned int multiplier = 1;
for (unsigned int i = 0; i < size; i++) {
result += level[i] * multiplier;
@@ -41,57 +40,52 @@ PlasticBarData::CalculateIndex(const vector<unsigned int>& level) {
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-double
-PS_PlasticBar::EnergyToLight(double Energy, int Z){
- //Get Energy to Light conversion factor depending on the particle doing the ionisation
+double PS_PlasticBar::EnergyToLight(double Energy, int Z) {
+ // Get Energy to Light conversion factor depending on the particle doing the ionisation
double a0 = t_Convertor[Z][0];
double a1 = t_Convertor[Z][1];
double a2 = t_Convertor[Z][2];
double a3 = t_Convertor[Z][3];
- double Light = a0*Energy - a1*( 1. - std::exp(-a2*std::pow(Energy,a3)) );
- return( std::max(Light ,0.) );
+ double Light = a0 * Energy - a1 * (1. - std::exp(-a2 * std::pow(Energy, a3)));
+ return (std::max(Light, 0.));
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void
-PS_PlasticBar::Compute_Light(){
-
+void PS_PlasticBar::Compute_Light() {
+
std::vector<PlasticBarData>::iterator it = m_Data.begin();
- std::vector<std::array<double,7>>::iterator it_Energy = t_Energy_by_ChargeNumber.begin();
- std::vector<std::array<double,7>>::iterator it_Light = t_Light_by_ChargeNumber.begin();
+ std::vector<std::array<double, 7>>::iterator it_Energy = t_Energy_by_ChargeNumber.begin();
+ std::vector<std::array<double, 7>>::iterator it_Light = t_Light_by_ChargeNumber.begin();
- for(; it_Energy < t_Energy_by_ChargeNumber.end() ; it_Energy++, it_Light++, it++){
- //cout << "Barre ID : " << (*it).GetIndex() << endl;
+ for (; it_Energy < t_Energy_by_ChargeNumber.end(); it_Energy++, it_Light++, it++) {
+ // cout << "Barre ID : " << (*it).GetIndex() << endl;
double Light = 0, Energy = 0;
double ConvertedLight = 0;
- //cout << " Energy Light Z " << endl;
- for(int i=0 ; i<7 ; i++){
+ // cout << " Energy Light Z " << endl;
+ for (int i = 0; i < 7; i++) {
Energy += (*it_Energy)[i];
- if(t_DoConversion){
- ConvertedLight = EnergyToLight( (*it_Energy)[i], i);
+ if (t_DoConversion) {
+ ConvertedLight = EnergyToLight((*it_Energy)[i], i);
(*it_Light)[i] = ConvertedLight;
Light += ConvertedLight;
- t_TotalLight_by_Z[i] += ConvertedLight;
+ t_TotalLight_by_Z[i] += ConvertedLight;
}
- //cout << " " << Energy;
- //cout << " " << Light;
- //cout << " " << i << endl;
+ // cout << " " << Energy;
+ // cout << " " << Light;
+ // cout << " " << i << endl;
}
- //cout << endl;
+ // cout << endl;
(*it).SetEnergy(Energy);
- if(t_DoConversion)
+ if (t_DoConversion)
(*it).SetLight(Light);
else
(*it).SetLight(Energy);
-
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-vector<PlasticBarData>::iterator
-PlasticBarDataVector::find(const unsigned int& index) {
- for (vector<PlasticBarData>::iterator it = m_Data.begin();
- it != m_Data.end(); it++) {
+vector<PlasticBarData>::iterator PlasticBarDataVector::find(const unsigned int& index) {
+ for (vector<PlasticBarData>::iterator it = m_Data.begin(); it != m_Data.end(); it++) {
if ((*it).GetIndex() == index)
return it;
}
@@ -99,14 +93,15 @@ PlasticBarDataVector::find(const unsigned int& index) {
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-PS_PlasticBar::PS_PlasticBar(G4String name, vector<G4int> NestingLevel, G4int depth)
- : G4VPrimitiveScorer(name, depth) {
+PS_PlasticBar::PS_PlasticBar(G4String name, vector<G4int> NestingLevel, G4int depth) : G4VPrimitiveScorer(name, depth) {
m_NestingLevel = NestingLevel;
auto tree = RootOutput::getInstance()->GetTree();
tree->Branch("PlasticBar_Time", &t_Time);
tree->Branch("PlasticBar_Position", &t_Position);
tree->Branch("PlasticBar_TotalEnergy_by_Z", &t_TotalEnergy_by_Z);
tree->Branch("PlasticBar_TotalLight_by_Z", &t_TotalLight_by_Z);
+ t_TotalEnergy_by_Z = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
+ t_TotalLight_by_Z = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -118,35 +113,36 @@ G4bool PS_PlasticBar::ProcessHits(G4Step* aStep, G4TouchableHistory*) {
// Contain Energy, Time + as many copy number as nested volume
unsigned int mysize = m_NestingLevel.size();
- // first lets consider a step going inside the bar, the interaction point can be associated to the position of the end of the step (which helps avoiding having all the steps at the entrance of the bar)
+ // first lets consider a step going inside the bar, the interaction point can be associated to the position of the end
+ // of the step (which helps avoiding having all the steps at the entrance of the bar)
G4String processname = aStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName();
- if(processname!="Transportation"){
+ if (processname != "Transportation") {
G4String particlename = aStep->GetTrack()->GetParticleDefinition()->GetParticleName();
int Z = aStep->GetTrack()->GetParticleDefinition()->GetAtomicNumber();
- //under the energy cut, neutrons can deposit energy via hadIon. But the conversion cant be done since we dont know the recoil particle. We chose to neglect their effect. (validated because deposited energy is less than 1MeV, which is the PM threshold)
- if(particlename!="neutron" && Z<7){
- t_Energy = aStep->GetTotalEnergyDeposit();
- t_Time = aStep->GetPostStepPoint()->GetGlobalTime();
- t_Position = aStep->GetPostStepPoint()->GetPosition()[1];
+ // under the energy cut, neutrons can deposit energy via hadIon. But the conversion cant be done since we dont know
+ // the recoil particle. We chose to neglect their effect. (validated because deposited energy is less than 1MeV,
+ // which is the PM threshold)
+ if (particlename != "neutron" && Z < 7) {
+ t_Energy = aStep->GetTotalEnergyDeposit();
+ t_Time = aStep->GetPostStepPoint()->GetGlobalTime();
+ t_Position = aStep->GetPostStepPoint()->GetPosition()[1];
t_Level.clear();
-
+
for (unsigned int i = 0; i < mysize; i++) {
- t_Level.push_back(
- aStep->GetPostStepPoint()->GetTouchableHandle()->GetCopyNumber(m_NestingLevel[i])
- );
+ t_Level.push_back(aStep->GetPostStepPoint()->GetTouchableHandle()->GetCopyNumber(m_NestingLevel[i]));
}
- if(t_Time < 150){
+ if (t_Time < 150) {
// Check if the bar has been hit before, if not, set the first interaction point.
std::vector<PlasticBarData>::iterator it;
it = m_Data.find(PlasticBarData::CalculateIndex(t_Level));
if (it == m_Data.end()) {
m_Data.Set(0.0, 0.0, t_Position, t_Time, t_Level);
- AddEntry(); //Adds null entry to t_Energy_by_ChargeNumber
- auto lastEnergyEntry = t_Energy_by_ChargeNumber.end() - 1; //iterator
+ AddEntry(); // Adds null entry to t_Energy_by_ChargeNumber
+ auto lastEnergyEntry = t_Energy_by_ChargeNumber.end() - 1; // iterator
(*lastEnergyEntry)[Z] += t_Energy;
t_TotalEnergy_by_Z[Z] += t_Energy;
}
@@ -155,16 +151,18 @@ G4bool PS_PlasticBar::ProcessHits(G4Step* aStep, G4TouchableHistory*) {
t_TotalEnergy_by_Z[Z] += t_Energy;
}
}
- }
- //if(it->GetLight()>160){
- // cout << "{!!!} Deposited Light : " << it->GetLight() << endl;
- //}
- //it->Add_Light(t_Light);
- //it->Set_Light(EnergyToLight(it->GetEnergy(),Z));
- //cout << particlename << " of Z = " << Z << " has left :" << endl;
- //cout << t_Energy << " MeV" << endl;
- //cout << EnergyToLight(it->GetEnergy(),Z)*MeV << " MeVEE" << endl;
- //now, if the next interaction is a Transportation, and there is an energy deposition, it means that we are observing the last step. To avoid having all the positions on the exit of the bar, lets consider the PreStepPosition instead.
+ }
+ // if(it->GetLight()>160){
+ // cout << "{!!!} Deposited Light : " << it->GetLight() << endl;
+ // }
+ // it->Add_Light(t_Light);
+ // it->Set_Light(EnergyToLight(it->GetEnergy(),Z));
+ // cout << particlename << " of Z = " << Z << " has left :" << endl;
+ // cout << t_Energy << " MeV" << endl;
+ // cout << EnergyToLight(it->GetEnergy(),Z)*MeV << " MeVEE" << endl;
+ // now, if the next interaction is a Transportation, and there is an energy deposition, it means that we are
+ // observing the last step. To avoid having all the positions on the exit of the bar, lets consider the
+ // PreStepPosition instead.
}
return TRUE;
}
@@ -173,9 +171,7 @@ G4bool PS_PlasticBar::ProcessHits(G4Step* aStep, G4TouchableHistory*) {
void PS_PlasticBar::Initialize(G4HCofThisEvent*) { clear(); }
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void PS_PlasticBar::EndOfEvent(G4HCofThisEvent*) {
- Compute_Light();
-}
+void PS_PlasticBar::EndOfEvent(G4HCofThisEvent*) { Compute_Light(); }
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void PS_PlasticBar::clear() {
diff --git a/NPSimulation/Scorers/PlasticBar.hh b/NPSimulation/Scorers/PlasticBar.hh
index baf6c1a8f908184fc342707504056096b1754fe3..691743be1066079af7f324eb75351d3c24e1ac71 100644
--- a/NPSimulation/Scorers/PlasticBar.hh
+++ b/NPSimulation/Scorers/PlasticBar.hh
@@ -26,152 +26,142 @@
*****************************************************************************/
#include "G4VPrimitiveScorer.hh"
#include "NPSHitsMap.hh"
-//#include "NPSecondaries.hh"
+// #include "NPSecondaries.hh"
+#include <array>
#include <map>
using namespace CLHEP;
namespace PlasticBar {
-
-
-
-
// Hold One hit info
- class PlasticBarData{
- public:
- PlasticBarData(const double& Energy, const double& Light, const double& Position, const double& Time,const std::vector<unsigned int>& Nesting){
- m_Index=CalculateIndex(Nesting);
- m_Level=Nesting;
- m_Energy=Energy;
- m_Light=Light;
- m_Position=Position;
- m_Time=Time;
- };
- ~PlasticBarData(){};
-
- private:
- unsigned int m_Index;
- std::vector<unsigned int> m_Level;
- double m_Energy;
- double m_Light;
- double m_Position;
- double m_Time;
-
- public:
- static unsigned int CalculateIndex(const std::vector<unsigned int>& Nesting);
-
- public:
- inline unsigned int GetIndex() const {return m_Index;};
- inline std::vector<unsigned int> GetLevel() const {return m_Level;};
- inline double GetEnergy() const {return m_Energy;};
- inline double GetLight() const {return m_Light;};
- inline double GetPosition() const {return m_Position;};
- inline double GetTime() const {return m_Time;};
-
- public:
- inline void SetEnergy(const double& Energy) {m_Energy = Energy;};
- inline void SetLight(const double& Light) {m_Light = Light;};
- };
-
-
-
-
+ class PlasticBarData {
+ public:
+ PlasticBarData(const double& Energy, const double& Light, const double& Position, const double& Time,
+ const std::vector<unsigned int>& Nesting) {
+ m_Index = CalculateIndex(Nesting);
+ m_Level = Nesting;
+ m_Energy = Energy;
+ m_Light = Light;
+ m_Position = Position;
+ m_Time = Time;
+ };
+ ~PlasticBarData(){};
+
+ private:
+ unsigned int m_Index;
+ std::vector<unsigned int> m_Level;
+ double m_Energy;
+ double m_Light;
+ double m_Position;
+ double m_Time;
+
+ public:
+ static unsigned int CalculateIndex(const std::vector<unsigned int>& Nesting);
+
+ public:
+ inline unsigned int GetIndex() const { return m_Index; };
+ inline std::vector<unsigned int> GetLevel() const { return m_Level; };
+ inline double GetEnergy() const { return m_Energy; };
+ inline double GetLight() const { return m_Light; };
+ inline double GetPosition() const { return m_Position; };
+ inline double GetTime() const { return m_Time; };
+
+ public:
+ inline void SetEnergy(const double& Energy) { m_Energy = Energy; };
+ inline void SetLight(const double& Light) { m_Light = Light; };
+ };
// Manage a vector of DSSD hit
- class PlasticBarDataVector{
- public:
- PlasticBarDataVector(){};
- ~PlasticBarDataVector(){};
-
- private:
- std::vector<PlasticBarData> m_Data;
-
- public:
- std::vector<PlasticBarData>::iterator find(const unsigned int& index);
- inline void clear(){m_Data.clear();} ;
- inline std::vector<PlasticBarData>::iterator end() {return m_Data.end();};
- inline std::vector<PlasticBarData>::iterator begin() {return m_Data.begin();};
- inline unsigned int size() {return m_Data.size();};
- PlasticBarData* operator[](const unsigned int& i){return &m_Data[i];};
-
- inline void Set(const double& Energy, const double& Light, const double& Position, const double& Time, const std::vector<unsigned int>& Nesting) {
- m_Data.push_back(PlasticBarData(Energy, Light, Position, Time, Nesting));
- };
+ class PlasticBarDataVector {
+ public:
+ PlasticBarDataVector(){};
+ ~PlasticBarDataVector(){};
+
+ private:
+ std::vector<PlasticBarData> m_Data;
+
+ public:
+ std::vector<PlasticBarData>::iterator find(const unsigned int& index);
+ inline void clear() { m_Data.clear(); };
+ inline std::vector<PlasticBarData>::iterator end() { return m_Data.end(); };
+ inline std::vector<PlasticBarData>::iterator begin() { return m_Data.begin(); };
+ inline unsigned int size() { return m_Data.size(); };
+ PlasticBarData* operator[](const unsigned int& i) { return &m_Data[i]; };
+
+ inline void Set(const double& Energy, const double& Light, const double& Position, const double& Time,
+ const std::vector<unsigned int>& Nesting) {
+ m_Data.push_back(PlasticBarData(Energy, Light, Position, Time, Nesting));
+ };
};
-
-
-
-
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
- class PS_PlasticBar : public G4VPrimitiveScorer{
-
- public: // with description
- PS_PlasticBar(G4String name, std::vector<G4int>NestingLevel, G4int depth=0);
- ~PS_PlasticBar();
-
- protected: // with description
- G4bool ProcessHits(G4Step*, G4TouchableHistory*);
- double EnergyToLight(double Energy, int Z);
- void Compute_Light();
-
- public:
- void Initialize(G4HCofThisEvent*);
- void EndOfEvent(G4HCofThisEvent*);
- void clear();
- void DrawAll();
- void PrintAll();
-
- private:
- // How much level of volume nesting should be considered
- // Give the list of the nesting level at which the copy number should be return.
- // 0 is the lowest level possible (the actual volume copy number in which the interaction happen)
- std::vector<G4int> m_NestingLevel;
-
- private:
- PlasticBarDataVector m_Data;
-
- private:
- std::vector< std::array<double,7> > t_Energy_by_ChargeNumber;
- std::vector< std::array<double,7> > t_Light_by_ChargeNumber;
- inline void AddEntry() {
- t_Energy_by_ChargeNumber.push_back({0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0});
- t_Light_by_ChargeNumber.push_back({0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0});
- };
- std::array<double,7> t_TotalEnergy_by_Z = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
- std::array<double,7> t_TotalLight_by_Z = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
-
- private:
- double t_Energy;
- double t_Light;
- double t_Position;
- double t_Time;
- std::vector<unsigned int> t_Level;
-
- private:
- bool t_DoConversion = 1;
-
- const double t_Convertor[7][4] = {
- {1.0 , 0.0 , 0.0 , 0.0 }, //electron
- {0.90 , 7.55 , 0.099 , 0.74}, //proton
- {0.78 , 39.3 , 0.022 , 0.91}, //alpha
- {0.61 , 57.1 , 0.012 , 0.95}, //7Li
- {0.44 , 45.5 , 0.010 , 0.92}, //9Be
- {0.35 , 34.5 , 0.011 , 0.91}, //11B
- {0.30 , 25.6 , 0.013 , 0.91}, //12C
- };
-
- public:
- inline unsigned int GetMult() {return m_Data.size();};
- inline double GetEnergy(const unsigned int& i) {return m_Data[i]->GetEnergy();};
- inline double GetLight(const unsigned int& i) {return m_Data[i]->GetLight();};
- inline double GetPosition(const unsigned int& i) {return m_Data[i]->GetPosition();};
- inline double GetTime(const unsigned int& i) {return m_Data[i]->GetTime();};
- inline std::vector<unsigned int> GetLevel(const unsigned int& i) {
- return m_Data[i]->GetLevel();
- };
+ class PS_PlasticBar : public G4VPrimitiveScorer {
+
+ public: // with description
+ PS_PlasticBar(G4String name, std::vector<G4int> NestingLevel, G4int depth = 0);
+ ~PS_PlasticBar();
+
+ protected: // with description
+ G4bool ProcessHits(G4Step*, G4TouchableHistory*);
+ double EnergyToLight(double Energy, int Z);
+ void Compute_Light();
+
+ public:
+ void Initialize(G4HCofThisEvent*);
+ void EndOfEvent(G4HCofThisEvent*);
+ void clear();
+ void DrawAll();
+ void PrintAll();
+
+ private:
+ // How much level of volume nesting should be considered
+ // Give the list of the nesting level at which the copy number should be return.
+ // 0 is the lowest level possible (the actual volume copy number in which the interaction happen)
+ std::vector<G4int> m_NestingLevel;
+
+ private:
+ PlasticBarDataVector m_Data;
+
+ private:
+ std::vector<std::array<double, 7>> t_Energy_by_ChargeNumber;
+ std::vector<std::array<double, 7>> t_Light_by_ChargeNumber;
+ inline void AddEntry() {
+ std::array<double, 7> nil_array = {{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
+ t_Energy_by_ChargeNumber.push_back(nil_array);
+ t_Light_by_ChargeNumber.push_back(nil_array);
+ };
+ std::array<double, 7> t_TotalEnergy_by_Z = {{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
+ std::array<double, 7> t_TotalLight_by_Z = {{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}};
+
+ private:
+ double t_Energy;
+ double t_Light;
+ double t_Position;
+ double t_Time;
+ std::vector<unsigned int> t_Level;
+
+ private:
+ bool t_DoConversion = 1;
+
+ const double t_Convertor[7][4] = {
+ {1.0, 0.0, 0.0, 0.0}, // electron
+ {0.90, 7.55, 0.099, 0.74}, // proton
+ {0.78, 39.3, 0.022, 0.91}, // alpha
+ {0.61, 57.1, 0.012, 0.95}, // 7Li
+ {0.44, 45.5, 0.010, 0.92}, // 9Be
+ {0.35, 34.5, 0.011, 0.91}, // 11B
+ {0.30, 25.6, 0.013, 0.91}, // 12C
+ };
+
+ public:
+ inline unsigned int GetMult() { return m_Data.size(); };
+ inline double GetEnergy(const unsigned int& i) { return m_Data[i]->GetEnergy(); };
+ inline double GetLight(const unsigned int& i) { return m_Data[i]->GetLight(); };
+ inline double GetPosition(const unsigned int& i) { return m_Data[i]->GetPosition(); };
+ inline double GetTime(const unsigned int& i) { return m_Data[i]->GetTime(); };
+ inline std::vector<unsigned int> GetLevel(const unsigned int& i) { return m_Data[i]->GetLevel(); };
};
-}
+} // namespace PlasticBar
#endif
diff --git a/Projects/Inelastic/238Uel.reaction b/Projects/Inelastic/238Uel.reaction
new file mode 100644
index 0000000000000000000000000000000000000000..8bb56b1d059533b59fd99eeb4baf60005abe2aaa
--- /dev/null
+++ b/Projects/Inelastic/238Uel.reaction
@@ -0,0 +1,28 @@
+Beam
+ Particle= n
+ Energy= 2.1
+ SigmaEnergy= 0.010
+ SigmaThetaX= 0.
+ SigmaPhiY= 0.
+ SigmaX= 10 mm
+ SigmaY= 10 mm
+ MeanThetaX= 0
+ MeanPhiY= 0
+ MeanX= 0
+ MeanY= 0
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+TwoBodyReaction
+ Beam= n
+ Target= 238U
+ Light= n
+ Heavy= 238U
+ ExcitationEnergyLight= 0.0
+ ExcitationEnergyHeavy= 0.0
+ CrossSectionPath= flat.txt CS
+ ShootLight= 1
+ ShootHeavy= 1
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
diff --git a/Projects/Inelastic/Analysis.cxx b/Projects/Inelastic/Analysis.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..452e686e36385936e9c8d8b30f8a545a8dd13004
--- /dev/null
+++ b/Projects/Inelastic/Analysis.cxx
@@ -0,0 +1,156 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2014 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: Adrien MATTA contact address: a.matta@surrey.ac.uk *
+ * *
+ * Creation Date : march 2025 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * Class describing the property of an Analysis object *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+#include<iostream>
+using namespace std;
+#include"Analysis.h"
+#include"NPAnalysisFactory.h"
+#include"NPDetectorManager.h"
+#include"NPOptionManager.h"
+#include"RootOutput.h"
+#include"RootInput.h"
+////////////////////////////////////////////////////////////////////////////////
+Analysis::Analysis(){
+}
+////////////////////////////////////////////////////////////////////////////////
+Analysis::~Analysis(){
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::Init(){
+ m_Vendeta= (TVendetaPhysics*) m_DetectorManager->GetDetector("Vendeta");
+ InitialConditions = new TInitialConditions();
+ InteractionCoordinates = new TInteractionCoordinates();
+ ReactionConditions = new TReactionConditions();
+
+ InitInputBranch();
+ InitOutputBranch();
+
+ my_Reaction = new NPL::Reaction("1n(238U,1n)238U@2.1");
+
+ neutron = new NPL::Nucleus("1n");
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::TreatEvent(){
+ ReInitValue();
+ Einit = InitialConditions->GetKineticEnergy(0);
+ double init_ThetaLab = ReactionConditions->GetTheta(0)*deg;
+ double init_BeamEnergy = ReactionConditions->GetBeamEnergy();
+ //my_Reaction->SetBeamEnergy(init_BeamEnergy);
+
+ neutron->SetKineticEnergy(Einit);
+ double beam_TOF = neutron->GetTimeOfFlight();
+
+ double Xtarget = InitialConditions->GetIncidentPositionX();
+ double Ytarget = InitialConditions->GetIncidentPositionY();
+ double Ztarget = 0;//InitialConditions->GetIncidentPositionZ();
+ TVector3 TargetPos = TVector3(Xtarget,Ytarget,Ztarget);
+
+ m_Vendeta->SetAnodeNumber(0);
+ m_Vendeta->BuildPhysicalEvent();
+ unsigned int size = m_Vendeta->LG_DetectorNumber.size();
+ for(int i=0; i<size; i++){
+ double Rdet, R;
+ Rdet = m_Vendeta->GetDistanceFromTarget(m_Vendeta->LG_DetectorNumber[i]);
+ TVector3 DetPos = m_Vendeta->GetVectorDetectorPosition(m_Vendeta->LG_DetectorNumber[i]);
+ //TVector3 HitPos = DetPos-TargetPos;
+
+ R= Rdet*mm;
+ Distance.push_back(R);
+ Det.push_back(m_Vendeta->LG_DetectorNumber[i]);
+ T.push_back(m_Vendeta->LG_Time[i]);
+ double T_stop = (m_Vendeta->LG_Time[i])*1e-9;
+ neutron->SetTimeOfFlight((T_stop-beam_TOF)/(Rdet*1e-3));
+ //neutron->SetTimeOfFlight((T_stop)/(Rdet*1e-3-8e-3));
+ Elab.push_back(neutron->GetEnergy());
+
+
+ double DeltaTheta = atan(63.5/Rdet);
+ double exp_ThetaLab = m_Vendeta->GetVectorDetectorPosition(m_Vendeta->LG_DetectorNumber[i]).Theta();
+ double random_ThetaLab = ra.Uniform(exp_ThetaLab-DeltaTheta, exp_ThetaLab+DeltaTheta);
+ double dEx = my_Reaction->ReconstructRelativistic(Elab[i], random_ThetaLab);
+
+ ThetaLab.push_back(random_ThetaLab/deg);
+ //ThetaLab.push_back(exp_ThetaLab/deg);
+ Ex.push_back(dEx);
+ }
+}
+
+///////////////////////////////////////////////////////////////////////////////
+void Analysis::InitOutputBranch() {
+ RootOutput::getInstance()->GetTree()->Branch("Einit",&Einit,"Einit/D");
+ RootOutput::getInstance()->GetTree()->Branch("ThetaLab",&ThetaLab);
+ RootOutput::getInstance()->GetTree()->Branch("Elab",&Elab);
+ RootOutput::getInstance()->GetTree()->Branch("Ex",&Ex);
+ RootOutput::getInstance()->GetTree()->Branch("T",&T);
+ RootOutput::getInstance()->GetTree()->Branch("Distance",&Distance);
+ RootOutput::getInstance()->GetTree()->Branch("Det",&Det);
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::InitInputBranch(){
+ RootInput:: getInstance()->GetChain()->SetBranchStatus("InitialConditions",true );
+ RootInput:: getInstance()->GetChain()->SetBranchStatus("fIC_*",true );
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("InitialConditions",&InitialConditions);
+
+ RootInput:: getInstance()->GetChain()->SetBranchStatus("ReactionConditions",true );
+ RootInput:: getInstance()->GetChain()->SetBranchStatus("fRC_*",true );
+ RootInput:: getInstance()->GetChain()->SetBranchAddress("ReactionConditions",&ReactionConditions);
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::ReInitValue(){
+ Einit = -100;
+ Ex.clear();
+ ThetaLab.clear();
+ Elab.clear();
+ T.clear();
+ Distance.clear();
+ Det.clear();
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void Analysis::End(){
+}
+
+
+////////////////////////////////////////////////////////////////////////////////
+// Construct Method to be pass to the DetectorFactory //
+////////////////////////////////////////////////////////////////////////////////
+NPL::VAnalysis* Analysis::Construct(){
+ return (NPL::VAnalysis*) new Analysis();
+}
+
+////////////////////////////////////////////////////////////////////////////////
+// Registering the construct method to the factory //
+////////////////////////////////////////////////////////////////////////////////
+extern "C"{
+ class proxy{
+ public:
+ proxy(){
+ NPL::AnalysisFactory::getInstance()->SetConstructor(Analysis::Construct);
+ }
+ };
+
+ proxy p;
+}
+
diff --git a/Projects/Inelastic/Analysis.h b/Projects/Inelastic/Analysis.h
new file mode 100644
index 0000000000000000000000000000000000000000..0d321075e8fb7868e0a79501c3720d3e4b87b6f3
--- /dev/null
+++ b/Projects/Inelastic/Analysis.h
@@ -0,0 +1,67 @@
+#ifndef Analysis_h
+#define Analysis_h
+/*****************************************************************************
+ * Copyright (C) 2009-2014 this file is part of the NPTool Project *
+ * *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
+ * For the list of contributors see $NPTOOL/Licence/Contributors *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: Adrien MATTA contact address: a.matta@surrey.ac.uk *
+ * *
+ * Creation Date : march 2025 *
+ * Last update : *
+ *---------------------------------------------------------------------------*
+ * Decription: *
+ * Class describing the property of an Analysis object *
+ * *
+ *---------------------------------------------------------------------------*
+ * Comment: *
+ * *
+ * *
+ *****************************************************************************/
+#include "NPVAnalysis.h"
+#include "TVendetaPhysics.h"
+#include "TInitialConditions.h"
+#include "TInteractionCoordinates.h"
+#include "TReactionConditions.h"
+#include "NPNucleus.h"
+#include "NPReaction.h"
+#include "TRandom3.h"
+
+class Analysis: public NPL::VAnalysis{
+ public:
+ Analysis();
+ ~Analysis();
+
+ public:
+ void Init();
+ void TreatEvent();
+ void End();
+ void InitOutputBranch();
+ void InitInputBranch();
+ void ReInitValue();
+ static NPL::VAnalysis* Construct();
+
+private:
+ double Einit;
+ vector<double> ThetaLab;
+ vector<double> Ex;
+ vector<double> Elab;
+ vector<double> T;
+ vector<double> Distance;
+ vector<int> Det;
+
+ private:
+ TVendetaPhysics* m_Vendeta;
+ TInitialConditions* InitialConditions;
+ TInteractionCoordinates* InteractionCoordinates;
+ TReactionConditions* ReactionConditions;
+
+ NPL::Reaction* my_Reaction;
+ NPL::Nucleus* neutron;
+
+ TRandom3 ra;
+};
+#endif
diff --git a/Projects/Inelastic/CMakeLists.txt b/Projects/Inelastic/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..22c74affdfc45019bdda2594f8439c52d4ab97ec
--- /dev/null
+++ b/Projects/Inelastic/CMakeLists.txt
@@ -0,0 +1,5 @@
+cmake_minimum_required (VERSION 2.8)
+# Setting the policy to match Cmake version
+cmake_policy(VERSION ${CMAKE_MAJOR_VERSION}.${CMAKE_MINOR_VERSION})
+# include the default NPAnalysis cmake file
+include("../../NPLib/ressources/CMake/NPAnalysis.cmake")
diff --git a/Projects/Inelastic/Geometry/GenerateGeometry.C b/Projects/Inelastic/Geometry/GenerateGeometry.C
new file mode 100644
index 0000000000000000000000000000000000000000..49a31993415710c40a2588050aba4864bce8bade
--- /dev/null
+++ b/Projects/Inelastic/Geometry/GenerateGeometry.C
@@ -0,0 +1,47 @@
+double R = 2500; // mm
+
+void GenerateGeometry()
+{
+
+ ofstream ofile;
+ string ofilename = "Vendeta.detector";
+ ofile.open(ofilename.c_str());
+ ofile << "\%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << endl;
+
+ ofile << "Target " << endl;
+ ofile << " THICKNESS= 10 micrometer" << endl;
+ ofile << " RADIUS= 20 mm" << endl;
+ ofile << " MATERIAL= 238U " << endl;
+ ofile << " ANGLE= 0 deg " << endl;
+ ofile << " X= 0 mm " << endl;
+ ofile << " Y= 0 mm " << endl;
+ ofile << " Z= 0 mm " << endl;
+
+ double Theta = 20;
+ double Phi[6] = {0, 20, 40, 140, 160, 180};
+ for(int p=0; p<6; p++){
+ if(p%2==0) Theta = 20;
+ else Theta = 15;
+ cout << p << " " << Theta << endl;
+ for(int i = 0; i <15; i++){
+ ofile << "\%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << endl;
+ ofile << "Vendeta" << endl;
+ ofile << " R= " << R << " mm" << endl;
+ ofile << " THETA= " << Theta << " deg" << endl;
+ ofile << " PHI= " << Phi[p] << " deg" << endl;
+
+ //ofile << "\%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << endl;
+ //ofile << "Vendeta" << endl;
+ //ofile << " R= " << R << " mm" << endl;
+ //ofile << " THETA= " << -Theta << " deg" << endl;
+ //ofile << " PHI= " << Phi << " deg" << endl;
+
+
+ Theta += 10;
+ }
+ }
+
+ ofile.close();
+
+
+}
diff --git a/Projects/Inelastic/PhysicsListOption.txt b/Projects/Inelastic/PhysicsListOption.txt
new file mode 100644
index 0000000000000000000000000000000000000000..ec1a9581f1eaf7645d75dab4f7a5ce99feb20aad
--- /dev/null
+++ b/Projects/Inelastic/PhysicsListOption.txt
@@ -0,0 +1,15 @@
+EmPhysicsList Option4
+DefaultCutOff 1
+DriftElectronPhysics 0
+IonBinaryCascadePhysics 0
+NPIonInelasticPhysics 0
+EmExtraPhysics 0
+HadronElasticPhysics 0
+StoppingPhysics 0
+OpticalPhysics 0
+HadronPhysicsINCLXX 0
+HadronPhysicsQGSP_BIC_HP 0
+HadronPhysicsQGSP_BERT_HP 0
+Decay 0
+Menate_R 0
+NeutronHP 1
diff --git a/Projects/Inelastic/RunToTreat.txt b/Projects/Inelastic/RunToTreat.txt
new file mode 100644
index 0000000000000000000000000000000000000000..d1ac9e045b470346e5cc1e4ba58670e7ff74a5b4
--- /dev/null
+++ b/Projects/Inelastic/RunToTreat.txt
@@ -0,0 +1,4 @@
+TTreeName
+ SimulatedTree
+RootFileName
+ ../../Outputs/Simulation/vendeta_el_*.root
diff --git a/Projects/Inelastic/Vendeta_inelastic.detector b/Projects/Inelastic/Vendeta_inelastic.detector
new file mode 100644
index 0000000000000000000000000000000000000000..777bad378275b19c72f7b11de4c6a9f665d8473a
--- /dev/null
+++ b/Projects/Inelastic/Vendeta_inelastic.detector
@@ -0,0 +1,459 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Target
+ THICKNESS= 10 micrometer
+ RADIUS= 20 mm
+ MATERIAL= 238U
+ ANGLE= 0 deg
+ X= 0 mm
+ Y= 0 mm
+ Z= 0 mm
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 20 deg
+ PHI= 0 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 30 deg
+ PHI= 0 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 40 deg
+ PHI= 0 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 50 deg
+ PHI= 0 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 60 deg
+ PHI= 0 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 70 deg
+ PHI= 0 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 80 deg
+ PHI= 0 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 90 deg
+ PHI= 0 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 100 deg
+ PHI= 0 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 110 deg
+ PHI= 0 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 120 deg
+ PHI= 0 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 130 deg
+ PHI= 0 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 140 deg
+ PHI= 0 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 150 deg
+ PHI= 0 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 160 deg
+ PHI= 0 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 15 deg
+ PHI= 20 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 25 deg
+ PHI= 20 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 35 deg
+ PHI= 20 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 45 deg
+ PHI= 20 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 55 deg
+ PHI= 20 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 65 deg
+ PHI= 20 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 75 deg
+ PHI= 20 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 85 deg
+ PHI= 20 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 95 deg
+ PHI= 20 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 105 deg
+ PHI= 20 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 115 deg
+ PHI= 20 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 125 deg
+ PHI= 20 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 135 deg
+ PHI= 20 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 145 deg
+ PHI= 20 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 155 deg
+ PHI= 20 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 20 deg
+ PHI= 40 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 30 deg
+ PHI= 40 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 40 deg
+ PHI= 40 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 50 deg
+ PHI= 40 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 60 deg
+ PHI= 40 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 70 deg
+ PHI= 40 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 80 deg
+ PHI= 40 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 90 deg
+ PHI= 40 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 100 deg
+ PHI= 40 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 110 deg
+ PHI= 40 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 120 deg
+ PHI= 40 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 130 deg
+ PHI= 40 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 140 deg
+ PHI= 40 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 150 deg
+ PHI= 40 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 160 deg
+ PHI= 40 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 15 deg
+ PHI= 140 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 25 deg
+ PHI= 140 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 35 deg
+ PHI= 140 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 45 deg
+ PHI= 140 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 55 deg
+ PHI= 140 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 65 deg
+ PHI= 140 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 75 deg
+ PHI= 140 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 85 deg
+ PHI= 140 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 95 deg
+ PHI= 140 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 105 deg
+ PHI= 140 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 115 deg
+ PHI= 140 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 125 deg
+ PHI= 140 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 135 deg
+ PHI= 140 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 145 deg
+ PHI= 140 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 155 deg
+ PHI= 140 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 20 deg
+ PHI= 160 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 30 deg
+ PHI= 160 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 40 deg
+ PHI= 160 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 50 deg
+ PHI= 160 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 60 deg
+ PHI= 160 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 70 deg
+ PHI= 160 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 80 deg
+ PHI= 160 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 90 deg
+ PHI= 160 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 100 deg
+ PHI= 160 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 110 deg
+ PHI= 160 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 120 deg
+ PHI= 160 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 130 deg
+ PHI= 160 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 140 deg
+ PHI= 160 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 150 deg
+ PHI= 160 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 160 deg
+ PHI= 160 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 15 deg
+ PHI= 180 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 25 deg
+ PHI= 180 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 35 deg
+ PHI= 180 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 45 deg
+ PHI= 180 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 55 deg
+ PHI= 180 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 65 deg
+ PHI= 180 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 75 deg
+ PHI= 180 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 85 deg
+ PHI= 180 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 95 deg
+ PHI= 180 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 105 deg
+ PHI= 180 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 115 deg
+ PHI= 180 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 125 deg
+ PHI= 180 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 135 deg
+ PHI= 180 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 145 deg
+ PHI= 180 deg
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Vendeta
+ R= 2500 mm
+ THETA= 155 deg
+ PHI= 180 deg
diff --git a/Projects/Inelastic/run.mac b/Projects/Inelastic/run.mac
new file mode 100644
index 0000000000000000000000000000000000000000..506556ed2882f22f4f6013ba592297aba63dec25
--- /dev/null
+++ b/Projects/Inelastic/run.mac
@@ -0,0 +1 @@
+/run/beamOn 100000
diff --git a/Projects/Inelastic/sim.sh b/Projects/Inelastic/sim.sh
new file mode 100755
index 0000000000000000000000000000000000000000..9fc763e80f426cd97d338787190ca3a183ce1cda
--- /dev/null
+++ b/Projects/Inelastic/sim.sh
@@ -0,0 +1,7 @@
+#!/bin/bash
+
+for i in {1..10}
+do
+ npsimulation -D Vendeta_inelastic.detector -E 238Uel.reaction --random-seed $i -O vendeta_el_$i -B run.mac &
+done
+
diff --git a/Projects/PISTA/Analysis.cxx b/Projects/PISTA/Analysis.cxx
index 5943a650cea221a4b76ccea7f493dbd830b3327c..a9c20e246d5bfeaca90cca35cbf5ac1aae480d73 100644
--- a/Projects/PISTA/Analysis.cxx
+++ b/Projects/PISTA/Analysis.cxx
@@ -78,7 +78,8 @@ void Analysis::TreatEvent(){
TVector3 BeamDirection = InitialConditions->GetBeamDirection();
TVector3 BeamPosition(XTarget,YTarget,ZTarget);
- TVector3 PositionOnTarget(0,0,0);
+ //TVector3 PositionOnTarget(0,0,0);
+ TVector3 PositionOnTarget(-1,0.5,0);
//TVector3 PositionOnTarget(Rand.Gaus(XTarget, 0.6/2.35), Rand.Gaus(YTarget, 0.6/2.35), 0);
//TVector3 PositionOnTarget(XTarget, YTarget, 0);
BeamEnergy = 1417.;//InitialConditions->GetIncidentInitialKineticEnergy();
@@ -113,7 +114,8 @@ void Analysis::TreatEvent(){
ThetaNormalTarget = HitDirection.Angle(TVector3(0,0,1));
//Elab = Be10C.EvaluateInitialEnergy(Energy,TargetThickness*0.5,ThetaNormalTarget);
- Elab = C12C.EvaluateInitialEnergy(Energy,TargetThickness*0.5,ThetaNormalTarget);
+ Elab = Energy;
+ //Elab = C12C.EvaluateInitialEnergy(Energy,TargetThickness*0.5,ThetaNormalTarget);
OptimumEx = Transfer->ReconstructRelativistic(OriginalElab, OriginalThetaLab*deg);
Ex = Transfer->ReconstructRelativistic(Elab, ThetaLab);
ThetaCM = Transfer->EnergyLabToThetaCM(Elab, ThetaLab)/deg;
diff --git a/Projects/Vendeta_sim/sim.sh b/Projects/Vendeta_sim/sim.sh
index 054ff32b4a29d126ce0c6be8a29057193d584c7f..3d7a2f27bb174949d7b773728f6dbbd3b81b55db 100755
--- a/Projects/Vendeta_sim/sim.sh
+++ b/Projects/Vendeta_sim/sim.sh
@@ -2,6 +2,6 @@
for i in {1..2}
do
- npsimulation -D Vendeta.detector -E neutron.source --random-seed $i -O vendeta_sim_$i -B run.mac &
+ npsimulation -D Vendeta_inelastic.detector -E 238Uel.reaction --random-seed $i -O vendeta_sim_$i -B run.mac &
done
diff --git a/Projects/e850/pista_e850_part1.detector b/Projects/e850/pista_e850_part1.detector
new file mode 100644
index 0000000000000000000000000000000000000000..539972b62974dec5dde6aced06caa869f78a2f24
--- /dev/null
+++ b/Projects/e850/pista_e850_part1.detector
@@ -0,0 +1,65 @@
+%%%%%%%%%%%%%%%%%%%%
+Target
+ THICKNESS= 0.44 micrometer
+ RADIUS= 5 mm
+ MATERIAL= C
+ ANGLE= 0 deg
+ X= 0 mm
+ Y= 0 mm
+ Z= 0 mm
+
+%%%%%%%%%%%%%%%%%%%%
+PISTA
+ POS_A= 36.71 96.1 50.7
+ POS_B= -35.92 95.8 51.3
+ POS_C= -19.35 56.5 93.4
+ POS_D= 20.8 56.6 93.1
+
+%%%%%%%%%%%%%%%%%%%%
+PISTA
+ POS_A= -42.15 93.1 51.1
+ POS_B= -93.62 41.8 51.6
+ POS_C= -54.12 25.3 93.6
+ POS_D= -25.66 53.6 93.4
+
+%%%%%%%%%%%%%%%%%%%%
+PISTA
+ POS_A= -95.98 35.6 51.3
+ POS_B= -95.63 -37.1 51.5
+ POS_C= -56.22 -20.7 93.7
+ POS_D= -56.39 19.5 93.6
+
+%%%%%%%%%%%%%%%%%%%%
+PISTA
+ POS_A= -93.34 -43.5 50.9
+ POS_B= -41.88 -94.8 51.3
+ POS_C= -25.41 -55.2 93.2
+ POS_D= -53.84 -26.8 93.0
+
+%%%%%%%%%%%%%%%%%%%%
+PISTA
+ POS_A= -35.66 -97.1 51.3
+ POS_B= 36.98 -96.9 51.1
+ POS_C= 20.9 -57.4 93.4
+ POS_D= -19.24 -57.6 93.5
+
+%%%%%%%%%%%%%%%%%%%%
+PISTA
+ POS_A= 43.66 -94.4 50.9
+ POS_B= 95.19 -43.3 50.8
+ POS_C= 56.03 -26.6 93.1
+ POS_D= 27.54 -54.9 93.2
+
+%%%%%%%%%%%%%%%%%%%%
+PISTA
+ POS_A= 96.87 -36.8 50.3
+ POS_B= 97.41 35.8 50.4
+ POS_C= 57.83 20.0 92.6
+ POS_D= 57.52 -20.1 92.6
+
+%%%%%%%%%%%%%%%%%%%%
+PISTA
+ POS_A= 94.39 42.0 50.8
+ POS_B= 42.94 93.3 50.7
+ POS_C= 26.47 54.1 93.0
+ POS_D= 54.9 25.7 93.1
diff --git a/Projects/e850/pista_e850_part2.detector b/Projects/e850/pista_e850_part2.detector
new file mode 100644
index 0000000000000000000000000000000000000000..86a3ec0571201e5a91ab7b1c0943f97259c632b3
--- /dev/null
+++ b/Projects/e850/pista_e850_part2.detector
@@ -0,0 +1,65 @@
+%%%%%%%%%%%%%%%%%%%%
+Target
+ THICKNESS= 0.44 micrometer
+ RADIUS= 5 mm
+ MATERIAL= C
+ ANGLE= 0 deg
+ X= 0 mm
+ Y= 0 mm
+ Z= 0 mm
+
+%%%%%%%%%%%%%%%%%%%%
+PISTA
+ POS_A= 36.58 96.2 50.7
+ POS_B= -36.06 96.1 51.1
+ POS_C= -19.69 56.8 93.4
+ POS_D= 20.46 56.8 93.1
+
+%%%%%%%%%%%%%%%%%%%%
+PISTA
+ POS_A= -42.32 93.5 51.1
+ POS_B= -93.68 42.2 51.5
+ POS_C= -54.14 25.7 93.5
+ POS_D= -25.74 54.1 93.3
+
+%%%%%%%%%%%%%%%%%%%%
+PISTA
+ POS_A= -96.11 35.8 51.3
+ POS_B= -95.9 -36.9 51.5
+ POS_C= -56.48 -20.6 93.7
+ POS_D= -56.58 19.6 93.6
+
+%%%%%%%%%%%%%%%%%%%%
+PISTA
+ POS_A= -93.58 -43.6 51.0
+ POS_B= -42.0 -94.7 51.3
+ POS_C= -25.55 -55.1 93.2
+ POS_D= -54.04 -26.8 93.1
+
+%%%%%%%%%%%%%%%%%%%%
+PISTA
+ POS_A= -35.96 -97.3 51.4
+ POS_B= 36.67 -97.1 50.9
+ POS_C= 20.73 -57.2 92.8
+ POS_D= -19.41 -57.3 93.1
+
+%%%%%%%%%%%%%%%%%%%%
+PISTA
+ POS_A= 43.74 -94.7 50.8
+ POS_B= 95.19 -43.5 50.8
+ POS_C= 55.98 -26.8 93.0
+ POS_D= 27.53 -55.1 93.0
+
+%%%%%%%%%%%%%%%%%%%%
+PISTA
+ POS_A= 97.15 -37.1 50.1
+ POS_B= 97.49 35.6 49.9
+ POS_C= 57.64 19.8 91.9
+ POS_D= 57.43 -20.3 92.0
+
+%%%%%%%%%%%%%%%%%%%%
+PISTA
+ POS_A= 94.08 42.0 51.0
+ POS_B= 42.78 93.4 50.9
+ POS_C= 26.27 54.4 93.4
+ POS_D= 54.6 25.9 93.4