Skip to content
Snippets Groups Projects
Commit 4b03f46b authored by Adrien Matta's avatar Adrien Matta :skull_crossbones:
Browse files

* Adding new project MUGAST_LISE

parent b35e1187
No related branches found
No related tags found
No related merge requests found
......@@ -396,9 +396,9 @@ void GaspardTrackerTrapezoid::ConstructDetector(G4LogicalVolume* world)
// w perpendicular to (u,v) plan and pointing ThirdStage
// Phi is angle between X axis and projection in (X,Y) plan
// Theta is angle between position vector and z axis
G4double wX = m_R[i] * sin(Theta / rad) * cos(Phi / rad);
G4double wY = m_R[i] * sin(Theta / rad) * sin(Phi / rad);
G4double wZ = m_R[i] * cos(Theta / rad);
G4double wX = m_R[i] * sin(Theta) * cos(Phi);
G4double wY = m_R[i] * sin(Theta) * sin(Phi);
G4double wZ = m_R[i] * cos(Theta);
MMw = G4ThreeVector(wX, wY, wZ);
// vector corresponding to the center of the module
......@@ -406,9 +406,9 @@ void GaspardTrackerTrapezoid::ConstructDetector(G4LogicalVolume* world)
// vector parallel to one axis of silicon plane
// in fact, this is vector u
G4double ii = cos(Theta / rad) * cos(Phi / rad);
G4double jj = cos(Theta / rad) * sin(Phi / rad);
G4double kk = -sin(Theta / rad);
G4double ii = cos(Theta) * cos(Phi);
G4double jj = cos(Theta) * sin(Phi);
G4double kk = -sin(Theta);
G4ThreeVector Y = G4ThreeVector(ii, jj, kk);
MMw = MMw.unit();
......
%%%%%%%%%%%%%%%%%%%%%% S1107 at Triumf %%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Beam
Particle= 68Ni
Energy= 1700
SigmaEnergy= 10
ExcitationEnergy= 0 MeV
SigmaThetaX= 3 mrad
SigmaPhiY= 3 mrad
SigmaX= 25 mm
SigmaY= 25 mm
MeanThetaX= 1 mrad
MeanPhiY= 1 mrad
MeanX= 0 mm
MeanY= 0 mm
%EnergyProfilePath= eLise.txt EL
%XThetaXProfilePath=
%YPhiYProfilePath=
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
TwoBodyReaction
Beam= 68Ni
Target= 2H
Light= 1H
Heavy= 69Ni
ExcitationEnergyLight= 0.0 MeV
ExcitationEnergyHeavy= 1.0 MeV
CrossSectionPath= fe61_f5_2_250.txt h
ShootLight= 1
ShootHeavy= 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
/*****************************************************************************
* 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 "NPFunction.h"
#include "NPAnalysisFactory.h"
#include "NPDetectorManager.h"
#include "NPOptionManager.h"
////////////////////////////////////////////////////////////////////////////////
Analysis::Analysis(){
}
////////////////////////////////////////////////////////////////////////////////
Analysis::~Analysis(){
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::Init() {
// initialize input and output branches
InitOutputBranch();
InitInputBranch();
myInit = new TInitialConditions();
// get MUST2 and Gaspard objects
M2 = (TMust2Physics*) m_DetectorManager -> GetDetector("M2Telescope");
GD = (GaspardTracker*) m_DetectorManager -> GetDetector("GaspardTracker");
// get reaction information
myReaction.ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
OriginalBeamEnergy = myReaction.GetBeamEnergy();
// target thickness
TargetThickness = m_DetectorManager->GetTargetThickness();
string TargetMaterial = m_DetectorManager->GetTargetMaterial();
// Cryo target case
WindowsThickness = m_DetectorManager->GetWindowsThickness();
string WindowsMaterial = m_DetectorManager->GetWindowsMaterial();
// energy losses
string light=NPL::ChangeNameToG4Standard(myReaction.GetNucleus3().GetName());
string beam=NPL::ChangeNameToG4Standard(myReaction.GetNucleus1().GetName());
LightTarget = NPL::EnergyLoss(light+"_"+TargetMaterial+".G4table","G4Table",100 );
LightAl = NPL::EnergyLoss(light+"_Al.G4table","G4Table",100);
LightSi = NPL::EnergyLoss(light+"_Si.G4table","G4Table",100);
BeamCD2 = NPL::EnergyLoss(beam+"_"+TargetMaterial+".G4table","G4Table",100);
if(WindowsThickness){
BeamWindow= new NPL::EnergyLoss(beam+"_"+WindowsMaterial+".G4table","G4Table",100);
LightWindow= new NPL::EnergyLoss(light+"_"+WindowsMaterial+".G4table","G4Table",100);
}
else{
BeamWindow= NULL;
LightWindow=NULL;
}
// initialize various parameters
Rand = TRandom3();
DetectorNumber = 0;
ThetaNormalTarget = 0;
ThetaM2Surface = 0;
Si_E_M2 = 0;
CsI_E_M2 = 0;
Energy = 0;
ThetaGDSurface = 0;
X = 0;
Y = 0;
Z = 0;
dE = 0;
dTheta=0;
BeamDirection = TVector3(0,0,1);
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::TreatEvent() {
// Reinitiate calculated variable
ReInitValue();
//double zImpact = Rand.Uniform(-TargetThickness*0.5,TargetThickness*0.5);
double zImpact = 0 ;
BeamImpact = TVector3(0,0,zImpact);
// determine beam energy for a randomized interaction point in target
// 1% FWHM randominastion (E/100)/2.35
myReaction.SetBeamEnergy(Rand.Gaus(myInit->GetIncidentFinalKineticEnergy(),myInit->GetIncidentFinalKineticEnergy()/235));
//////////////////////////// LOOP on MUST2 //////////////////
for(unsigned int countMust2 = 0 ; countMust2 < M2->Si_E.size() ; countMust2++){
/************************************************/
//Part 0 : Get the usefull Data
// MUST2
int TelescopeNumber = M2->TelescopeNumber[countMust2];
/************************************************/
// Part 1 : Impact Angle
ThetaM2Surface = 0;
ThetaNormalTarget = 0;
TVector3 HitDirection = M2 -> GetPositionOfInteraction(countMust2) - BeamImpact ;
ThetaLab = HitDirection.Angle( BeamDirection );
X = M2 -> GetPositionOfInteraction(countMust2).X();
Y = M2 -> GetPositionOfInteraction(countMust2).Y();
Z = M2 -> GetPositionOfInteraction(countMust2).Z();
ThetaM2Surface = HitDirection.Angle(- M2 -> GetTelescopeNormal(countMust2) );
ThetaNormalTarget = HitDirection.Angle( TVector3(0,0,1) ) ;
/************************************************/
/************************************************/
// Part 2 : Impact Energy
Energy = ELab = 0;
Si_E_M2 = M2->Si_E[countMust2];
CsI_E_M2= M2->CsI_E[countMust2];
// if CsI
if(CsI_E_M2>0 ){
// The energy in CsI is calculate form dE/dx Table because
Energy = CsI_E_M2;
Energy = LightAl.EvaluateInitialEnergy( Energy ,0.4*micrometer , ThetaM2Surface);
Energy+=Si_E_M2;
}
else
Energy = Si_E_M2;
// Evaluate energy using the thickness
ELab = LightAl.EvaluateInitialEnergy( Energy ,0.4*micrometer , ThetaM2Surface);
// Target Correction
ELab = LightTarget.EvaluateInitialEnergy( ELab ,TargetThickness/2.-zImpact, ThetaNormalTarget);
if(LightWindow)
ELab = LightWindow->EvaluateInitialEnergy( ELab ,WindowsThickness, ThetaNormalTarget);
/************************************************/
/************************************************/
// Part 3 : Excitation Energy Calculation
Ex = myReaction.ReconstructRelativistic( ELab , ThetaLab );
ThetaLab=ThetaLab/deg;
/************************************************/
/************************************************/
// Part 4 : Theta CM Calculation
ThetaCM = myReaction.EnergyLabToThetaCM( ELab , ThetaLab)/deg;
/************************************************/
}//end loop MUST2
////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////
//////////////////////////// LOOP on GASPARD //////////////////
if(GD->GetEnergyDeposit()>0){
/************************************************/
// Part 1 : Impact Angle
ThetaGDSurface = 0;
ThetaNormalTarget = 0;
TVector3 HitDirection = GD -> GetPositionOfInteraction() - BeamImpact ;
ThetaLab = HitDirection.Angle( BeamDirection );
X = GD -> GetPositionOfInteraction().X();
Y = GD -> GetPositionOfInteraction().Y();
Z = GD -> GetPositionOfInteraction().Z();
ThetaGDSurface = HitDirection.Angle( TVector3(0,0,1) ) ;
ThetaNormalTarget = HitDirection.Angle( TVector3(0,0,1) ) ;
/************************************************/
/************************************************/
// Part 2 : Impact Energy
Energy = ELab = 0;
Energy = GD->GetEnergyDeposit();
// Target Correction
ELab = LightTarget.EvaluateInitialEnergy( Energy ,TargetThickness*0.5-zImpact, ThetaNormalTarget);
if(LightWindow)
ELab = LightWindow->EvaluateInitialEnergy( ELab ,WindowsThickness, ThetaNormalTarget);
dE= myInit->GetKineticEnergy(0) - ELab ;
dTheta = (myInit->GetThetaLab_WorldFrame(0)*deg-ThetaLab)/deg ;
/************************************************/
/************************************************/
// Part 3 : Excitation Energy Calculation
Ex = myReaction.ReconstructRelativistic( ELab , ThetaLab );
/************************************************/
/************************************************/
// Part 4 : Theta CM Calculation
ThetaCM = myReaction.EnergyLabToThetaCM( ELab , ThetaLab)/deg;
ThetaLab=ThetaLab/deg;
/************************************************/
}//end loop GASPARD
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::End(){
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::InitOutputBranch() {
RootOutput::getInstance()->GetTree()->Branch("Ex",&Ex,"Ex/D");
RootOutput::getInstance()->GetTree()->Branch("ELab",&ELab,"ELab/D");
RootOutput::getInstance()->GetTree()->Branch("ThetaLab",&ThetaLab,"ThetaLab/D");
RootOutput::getInstance()->GetTree()->Branch("ThetaCM",&ThetaCM,"ThetaCM/D");
RootOutput::getInstance()->GetTree()->Branch("Run",&Run,"Run/I");
RootOutput::getInstance()->GetTree()->Branch("X",&X,"X/D");
RootOutput::getInstance()->GetTree()->Branch("Y",&Y,"Y/D");
RootOutput::getInstance()->GetTree()->Branch("Z",&Z,"Z/D");
RootOutput::getInstance()->GetTree()->Branch("dE",&dE,"dE/D");
RootOutput::getInstance()->GetTree()->Branch("dTheta",&dTheta,"dTheta/D");
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::InitInputBranch(){
RootInput::getInstance()->GetChain()->SetBranchStatus("Run",true);
RootInput::getInstance()->GetChain()->SetBranchAddress("Run",&Run);
RootInput::getInstance()->GetChain()->SetBranchStatus("InitialConditions",true);
RootInput::getInstance()->GetChain()->SetBranchAddress("InitialConditions",&myInit);
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::ReInitValue(){
Ex = -1000 ;
ELab = -1000;
ThetaLab = -1000;
ThetaCM = -1000;
X = -1000;
Y = -1000;
Z = -1000;
dE= -1000;
dTheta= -1000;
}
////////////////////////////////////////////////////////////////////////////////
// Construct Method to be pass to the AnalysisFactory //
////////////////////////////////////////////////////////////////////////////////
NPL::VAnalysis* Analysis::Construct(){
return (NPL::VAnalysis*) new Analysis();
}
////////////////////////////////////////////////////////////////////////////////
// Registering the construct method to the factory //
////////////////////////////////////////////////////////////////////////////////
extern "C"{
class proxy_analysis{
public:
proxy_analysis(){
NPL::AnalysisFactory::getInstance()->SetConstructor(Analysis::Construct);
}
};
proxy_analysis p_analysis;
}
#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"NPEnergyLoss.h"
#include"NPReaction.h"
#include"RootOutput.h"
#include"RootInput.h"
#include "TMust2Physics.h"
#include "GaspardTracker.h"
#include "TInitialConditions.h"
#include <TRandom3.h>
#include <TVector3.h>
#include <TMath.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 Ex;
double ELab;
double ThetaLab;
double ThetaCM;
NPL::Reaction myReaction;
// Energy loss table: the G4Table are generated by the simulation
NPL::EnergyLoss LightTarget;
NPL::EnergyLoss LightAl;
NPL::EnergyLoss LightSi;
NPL::EnergyLoss BeamCD2;
NPL::EnergyLoss* BeamWindow;
NPL::EnergyLoss* LightWindow;
double TargetThickness ;
double WindowsThickness;
// Beam Energy
double OriginalBeamEnergy ; // AMEV
double FinalBeamEnergy;
// intermediate variable
TVector3 BeamDirection;
TVector3 BeamImpact;
TRandom3 Rand ;
int Run;
int DetectorNumber ;
double ThetaNormalTarget;
double ThetaM2Surface ;
double Si_E_M2 ;
double CsI_E_M2 ;
double Energy ;
double ThetaGDSurface ;
double X ;
double Y ;
double Z ;
double dE;
double dTheta;
// Branches and detectors
TMust2Physics* M2;
GaspardTracker* GD;
TInitialConditions* myInit ;
};
#endif
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")
%%%%%%%%%%Detector%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1
Target
THICKNESS= 9.5
RADIUS= 7.5 mm
MATERIAL= CD2
Angle= 0 deg
X= 0 mm
Y= 0 mm
Z= 0 mm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
AnnularS1
Z= -100
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
GaspardTracker Trapezoid
R= -90 mm
THETA= 60 deg
PHI= 0 deg
BETA= -60 0 -90 deg
FIRSTSTAGE= 1
SECONDSTAGE= 0
THIRDSTAGE= 0
VIS= all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%GaspardTracker Trapezoid
R= -90 mm
THETA= 60 deg
PHI= 45 deg
BETA= -60 0 -90 deg
FIRSTSTAGE= 1
SECONDSTAGE= -1
THIRDSTAGE= 0
VIS= all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%GaspardTracker Trapezoid
R= -90 mm
THETA= 60 deg
PHI= 90 deg
BETA= -60 0 -90 deg
FIRSTSTAGE= 1
SECONDSTAGE= 0
THIRDSTAGE= 0
VIS= all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%GaspardTracker Trapezoid
R= -90 mm
THETA= 60 deg
PHI= 135 deg
BETA= -60 0 -90 deg
FIRSTSTAGE= 1
SECONDSTAGE= 0
THIRDSTAGE= 0
VIS= all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%GaspardTracker Trapezoid
R= -90 mm
THETA= 60 deg
PHI= 180 deg
BETA= -60 0 -90 deg
FIRSTSTAGE= 1
SECONDSTAGE= 0
THIRDSTAGE= 0
VIS= all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%GaspardTracker Trapezoid
R= -90 mm
THETA= 60 deg
PHI= 225 deg
BETA= -60 0 -90 deg
FIRSTSTAGE= 1
SECONDSTAGE= 0
THIRDSTAGE= 0
VIS= all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%GaspardTracker Trapezoid
R= -90 mm
THETA= 60 deg
PHI= 270 deg
BETA= -60 0 -90 deg
FIRSTSTAGE= 1
SECONDSTAGE= 0
THIRDSTAGE= 0
VIS= all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%GaspardTracker Trapezoid
R= -90 mm
THETA= 60 deg
PHI= 315 deg
BETA= -60 0 -90 deg
FIRSTSTAGE= 1
SECONDSTAGE= 0
THIRDSTAGE= 0
VIS= all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Tiara Barrel
InnerBarrel= 0
OuterBarrel= 0
Chamber= 0
X= 0 mm
Y= 0 mm
Z= 0 mm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%% An Isotropic Source to be used as EventGenerator %%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Energy are given in MeV , Position in mm %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Isotropic
EnergyLow= 5 MeV
EnergyHigh= 5 MeV
HalfOpenAngleMin= 90 deg
HalfOpenAngleMax= 180 deg
x0= 0 mm
y0= 0 mm
z0= 0 mm
Particle= alpha
ExcitationEnergy= 0 MeV
% Supported particle type: proton, neutron, deuton, triton, He3 , alpha
......@@ -88,19 +88,20 @@ void GeometricalEfficiency(const char * fname = "myResult.root"){
}
hEmittTheta->Sumw2();
hEmittThetaCM->Sumw2();
hDetecTheta->Sumw2();
hDetecThetaCM->Sumw2();
hEmittThetaCM->Sumw2();
hDetecTheta->Sumw2();
hDetecThetaCM->Sumw2();
TCanvas *c0 = new TCanvas("c0", "Distrib",800,800);
hEmittTheta->Draw("");
TCanvas *c0 = new TCanvas("c0", "Distrib",800,800);
hEmittTheta->Draw("");
hDetecTheta->SetMarkerColor(kAzure+7);
hDetecTheta->Draw("same");
// efficiency in lab frame in %
TCanvas *c = new TCanvas("c", "efficiency",800,800);
c->SetTopMargin(0.01);
c->SetRightMargin(0.03);
TH1F *hEfficiency = new TH1F("hEfficiency", "Efficiency", 180, 0, 90);
TH1F *hEfficiency = new TH1F("hEfficiency", "Efficiency", 180, 0, 180);
hEfficiency->Divide(hDetecTheta, hEmittTheta, 1, 1);
hEfficiency->GetXaxis()->SetTitle("#Theta (deg)");
hEfficiency->GetYaxis()->SetTitle("#epsilon");
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment