Skip to content
Snippets Groups Projects
Commit d137e717 authored by Morfouace's avatar Morfouace
Browse files

*removing obsolete Hira folder in the analysis area

parent fbcc5715
No related branches found
No related tags found
No related merge requests found
#include "Analysis.h"
using namespace std;
int main(int argc, char** argv)
{
// command line parsing
NPOptionManager* myOptionManager = NPOptionManager::getInstance(argc,argv);
// Instantiate RootInput
string runToReadfileName = myOptionManager->GetRunToReadFile();
RootInput:: getInstance(runToReadfileName);
// if input files are not given, use those from TAsciiFile
if (myOptionManager->IsDefault("DetectorConfiguration")) {
string name = RootInput::getInstance()->DumpAsciiFile("DetectorConfiguration");
myOptionManager->SetDetectorFile(name);
}
// get input files from NPOptionManager
string detectorfileName = myOptionManager->GetDetectorFile();
string calibrationfileName = myOptionManager->GetCalibrationFile();
string OutputfileName = myOptionManager->GetOutputFile();
// Instantiate RootOutput
RootOutput::getInstance("Analysis/"+OutputfileName, "AnalysedTree");
// Instantiate the detector using a file
NPL::DetectorManager* myDetector = new DetectorManager();
myDetector->ReadConfigurationFile(detectorfileName);
// Get the formed Chained Tree and Treat it
TChain* Chain = RootInput:: getInstance() -> GetChain();
// Instantiate some Reaction
NPL::Reaction* TransfertReaction = new Reaction ;
TransfertReaction -> ReadConfigurationFile("34Ar_pd.reaction") ;
//TransfertReaction -> ReadConfigurationFile("46Ar_pd.reaction") ;
//TransfertReaction -> ReadConfigurationFile("11Be_d3He.reaction") ;
//Get Detector pointer :
THiraPhysics* Hira = (THiraPhysics*) myDetector -> GetDetector("HIRAArray") ;
//TInteractionCoordinate
TInteractionCoordinates* InteractionCoordinates = new TInteractionCoordinates();
Chain->SetBranchStatus("InteractionCoordinates",true);
Chain->SetBranchStatus("fDetected_*",true);
Chain->SetBranchAddress("InteractionCoordinates",&InteractionCoordinates);
RootOutput::getInstance()->GetTree()->Branch("InteractionCoordinates","TInteractionCoordinates",&InteractionCoordinates);
//TInitialConditions
TInitialConditions* InitialConditions = new TInitialConditions();
Chain->SetBranchStatus("InitialConditions",true);
Chain->SetBranchStatus("fIC_*",true);
Chain->SetBranchAddress("InitialConditions",&InitialConditions);
RootOutput::getInstance()->GetTree()->Branch("InitialConditions","TInitialConditions",&InitialConditions);
RootOutput::getInstance()->GetTree()->Branch( "ThetaLab" , &ThetaLab , "ThetaLab/D" );
RootOutput::getInstance()->GetTree()->Branch( "PhiLab" , &PhiLab , "PhiLab/D" ) ;
RootOutput::getInstance()->GetTree()->Branch("ThetaCM", &ThetaCM,"ThetaCM/D") ;
RootOutput::getInstance()->GetTree()->Branch( "E_ThinSi" , &E_ThinSi , "E_ThinSi/D" ) ;
RootOutput::getInstance()->GetTree()->Branch( "E_ThickSi" , &E_ThickSi , "E_ThickSi/D" ) ;
RootOutput::getInstance()->GetTree()->Branch( "E_CsI" , &E_CsI , "E_CsI/D" ) ;
RootOutput::getInstance()->GetTree()->Branch( "ELab" , &ELab , "ELab/D" ) ;
RootOutput::getInstance()->GetTree()->Branch("ExcitationEnergy", &ExcitationEnergy,"ExcitationEnergy/D") ;
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( "TelescopeNumber" , &TelescopeNumber , "TelescopeNumber/D" ) ;
// Get number of events to treat
cout << endl << "///////// Starting Analysis ///////// "<< endl;
int nentries = Chain->GetEntries();
cout << " Number of Event to be treated : " << nentries << endl;
clock_t begin = clock();
clock_t end = begin;
double EnergyThreshold = 0.;
// main loop on entries
for (int i = 0; i < nentries; i++) {
if (i%10000 == 0 && i!=0) {
cout.precision(5);
end = clock();
double TimeElapsed = (end-begin) / CLOCKS_PER_SEC;
double percent = (double)i/nentries;
double TimeToWait = (TimeElapsed/percent) - TimeElapsed;
cout << " "<< flush;
cout << "\r Progression:" << percent*100 << " % \t | \t Remaining time : ~" << TimeToWait <<"s"<< flush;
}
else if (i == nentries-1) cout << "\r Progression:" << " 100% " <<endl;
// get data
Chain -> GetEntry(i);
ReInitOuputValue();
myDetector->ClearEventPhysics();
myDetector->BuildPhysicalEvent();
if(Hira -> ThickSi_E.size() == 1){
for(int countHira = 0 ; countHira < Hira->ThickSi_E.size() ; countHira++){
TelescopeNumber = Hira->TelescopeNumber[countHira];
TargetThickness = myDetector->GetTargetThickness();
X = Hira->GetPositionOfInteraction(0).X();
Y = Hira->GetPositionOfInteraction(0).Y();
Z = Hira->GetPositionOfInteraction(0).Z();
TVector3 PositionOnHira = TVector3(X,Y,Z);
TVector3 ZUnit = TVector3(0,0,1);
//TVector3 TelescopeNormal = -Hira->GetTelescopeNormal(countHira);
TVector3 TargetNormal = TVector3(0,0,1);//TVector3(-sin(TargetAngle*deg),0,cos(TargetAngle*deg));
double X_target = InitialConditions->GetIncidentPositionX();
double Y_target = InitialConditions->GetIncidentPositionY();
double Z_target = InitialConditions->GetIncidentPositionZ();
TVector3 PositionOnTarget = TVector3(X_target,Y_target,Z_target);
//TVector3 PositionOnTarget = TVector3(0,0,0);
TVector3 HitDirection = PositionOnHira - PositionOnTarget;
TVector3 HitDirectionUnit = HitDirection.Unit();
//double ThetaNormalHira = HitDirection.Angle(TelescopeNormal);
double ThetaNormalTarget = HitDirection.Angle(TargetNormal);
//double ThetaBeam = InitialConditions->GetICIncidentAngleTheta(countHira);
//double PhiBeam = InitialConditions->GetICIncidentAnglePhi(countHira);
TVector3 BeamDirection = InitialConditions->GetBeamDirection();
//TVector3(sin(ThetaBeam*deg)*cos(PhiBeam*deg),sin(ThetaBeam*deg)*sin(PhiBeam*deg),cos(ThetaBeam*deg));
double XBeam = BeamDirection.X();
double YBeam = BeamDirection.Y();
double ZBeam = BeamDirection.Z();
ThetaLab = BeamDirection.Angle(HitDirection);
ThetaLab = ThetaLab/deg;
PhiLab = HitDirection.Phi()/deg;
E_ThickSi = Hira->ThickSi_E[countHira];
E_ThinSi = Hira->ThinSi_E[countHira];
ELab = E_ThinSi + E_ThickSi;
if(Hira->CsI_E.size() == 1){
for(int countCsI =0; countCsI<Hira->CsI_E.size(); countCsI++){
E_CsI = Hira->CsI_E[countCsI];
if(E_CsI>EnergyThreshold)ELab += E_CsI;
}
}
else E_CsI = -1000;
ELab = EL_deuteron_CH2.EvaluateInitialEnergy(ELab,
(TargetThickness/2)*micrometer,
ThetaNormalTarget);
// ********************** Angle in the CM frame *****************************
TransfertReaction -> SetNuclei3(ELab, ThetaLab*deg);
ThetaCM = TransfertReaction->GetThetaCM()/deg;
ExcitationEnergy = TransfertReaction->GetExcitation4();
RootOutput::getInstance()->GetTree()->Fill();
}
}
}
cout << "A total of " << nentries << " event has been annalysed " << endl ;
RootOutput::getInstance()->Destroy();
RootInput::getInstance()->Destroy();
NPOptionManager::getInstance()->Destroy();
return 0 ;
}
/////////////////////////////////////////////////////////////////////////////////////////
// On reinitialise les valeurs des branches de sorties à chaque entrée dans la boucle :
void ReInitOuputValue()
{
E_ThinSi = -1000;
E_ThickSi = -1000;
E_CsI = -1000;
ELab = -1000;
ThetaLab = -1000;
PhiLab = -1000;
ThetaCM = -1000;
ExcitationEnergy = -1000;
X = -1000;
Y = -1000;
Z = -1000;
TelescopeNumber = -1;
}
// You can use this file to declare your spectra, file, energy loss , ... and whatever you want.
// This way you can remove all unnecessary declaration in the main programm.
// In order to help debugging and organizing we use Name Space.
/////////////////////////////////////////////////////////////////////////////////////////////////
// -------------------------------------- VARIOUS INCLUDE ---------------------------------------
// NPA
#include "NPDetectorManager.h"
#include "NPOptionManager.h"
// STL C++
#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <cmath>
#include <cstdlib>
// ROOT
#include <TROOT.h>
#include <TChain.h>
#include <TFile.h>
#include <TLeaf.h>
#include <TVector3.h>
#include <TRandom.h>
// NPL
#include "RootInput.h"
#include "RootOutput.h"
#include "NPReaction.h"
#include "TInitialConditions.h"
#include "TInteractionCoordinates.h"
#include "THiraData.h"
#include "THiraPhysics.h"
// Use CLHEP System of unit and Physical Constant
#include "NPGlobalSystemOfUnits.h"
#include "NPPhysicalConstants.h"
Double_t TargetThickness;
Double_t E_ThinSi = 0;
Double_t E_ThickSi = 0;
Double_t E_CsI = 0;
Double_t ELab = 0;
Double_t ThetaLab = 0;
Double_t PhiLab = 0;
Double_t ThetaLab_simu = 0;
Double_t ThetaCM = 0;
Double_t ExcitationEnergy = 0;
Double_t X = 0;
Double_t Y = 0;
Double_t Z = 0;
Double_t TelescopeNumber = 0;
Double_t BeamEnergy = 0;
void ReInitOuputValue();
// ----------------------------------------------------------------------------------------------
double ThetaCalculation (TVector3 A , TVector3 B) ;
/////////////////////////////////////////////////////////////////////////////////////////////////
// ----------------------------------- DOUBLE, INT, BOOL AND MORE -------------------------------
namespace VARIABLE
{
// Declare your Variable here:
double X1,Y1,Z1 ;
int N1,N2 = 0 ;
bool check= false ;
// A Usefull Simple Random Generator
TRandom Rand;
}
using namespace VARIABLE ;
// ----------------------------------------------------------------------------------------------
/////////////////////////////////////////////////////////////////////////////////////////////////
// -----------------------------------GRAPH------------------------------------------------------
#include <TObject.h>
#include <TH1.h>
#include <TH1F.h>
#include <TH2.h>
#include <TH2F.h>
#include <TGraph2D.h>
namespace GRAPH
{
// Declare your Spectra here:
TH1F *myHist1D = new TH1F("Hist1D","Histogramm 1D ; x ; count", 1000 , -5 , 5 ) ;
TH2F *myHist2D = new TH2F("Hist2D","Histogramm 2D ; x ; y ", 128 , 1 , 128 , 128 , 1 , 128 ) ;
}
using namespace GRAPH ;
// --------------------------------------------------------------------------------------------
///////////////////////////////////////////////////////////////////////////////////////////////
// -----------------------------------CUT------------------------------------------------------
#include <TCutG.h>
namespace CUT
{
// Declare your Cut here:
}
using namespace CUT ;
// --------------------------------------------------------------------------------------------
////////////////////////////////////////////////////////////////////////////////////////////////
// -----------------------------------ENERGY LOSS----------------------------------------------
#include "NPEnergyLoss.h"
using namespace NPL ;
namespace ENERGYLOSS
{
//EnergyLoss EL_deuteron_CH2 = EnergyLoss("deuteron_CH2.G4table","G4Table",100 );
EnergyLoss EL_deuteron_CH2 = EnergyLoss("He3_CD2.G4table","G4Table",100 );
// Declare your Energy loss here :
/* EnergyLoss DeuerontTarget = EnergyLoss ( "CD2.txt" ,
100 ,
1,
1 );
*/
}
using namespace ENERGYLOSS ;
// ----------------------------------------------------------------------------------------------
/////////////////////////////////////////////////////////////////////////////////////////////////
# include same architecture file than for NPLib
# so that consistency is ensured
include $(NPTOOL)/NPLib/Makefile.arch
# additional libraries
LIBRARY = `$(NPTOOL)/NPLib/liblist`
PROGRAMS = Analysis
all: $(PROGRAMS)
Analysis: Analysis.o
$(LD) $(LDFLAGS) $^ $(LIBS) $(LIBRARY) $(OutPutOpt) $@
@echo "$@ done"
# rule for creating .o from .cxx
.SUFFIXES: .$(SrcSuf)
.$(SrcSuf).$(ObjSuf):
$(CXX) $(CXXFLAGS) $(INCLUDE) -c $<
# some cleaning
clean:
rm -rf *.o
distclean:
make clean; rm $(PROGRAMS)
# dependences
Analysis.o: Analysis.cxx Analysis.h
TTreeName
SimulatedTree
RootFileName
%/Users/pierremorfouace/Physics/NPTool/nptool/Outputs/Simulation/46Ar_pd_gs.root
%/Users/pierremorfouace/Physics/NPTool/nptool/Outputs/Simulation/46Ar_pd_1st.root
/Users/pierremorfouace/Physics/NPTool/nptool/Outputs/Simulation/test.root
void InitChain();
void LoadCut();
TChain *chain;
TGraph *kin;
int nentries;
TCutG *cut_kine;
void Show()
{
gStyle->SetOptStat(0);
InitChain();
LoadCut();
nentries = chain->GetEntries();
//**************** Condition ****************//
TString CondCut = "cut_kine";
TString CondTot = CondCut;//+"&&"+
//*******************************************//
TCanvas* mainC0 = new TCanvas("Hira XY", "Hira XY" , 800,600);
TCanvas* mainC1 = new TCanvas("PID", "PID" , 1200,600);
TCanvas* mainC2 = new TCanvas("Kinematics", "Kinematics" , 800,600);
TCanvas* mainC3 = new TCanvas("Phi-Theta", "Phi-Theta" , 800,600);
TCanvas* mainC4 = new TCanvas("BeamSpot", "BeamSpot" , 800,600);
TCanvas* mainC5 = new TCanvas("ThetaCM", "ThetaCM" , 800,600);
TCanvas* mainC6 = new TCanvas("Efficiency", "Efficiency" , 800,600);
TCanvas* mainC7 = new TCanvas("CrossSection","CrossSection",800,600);
mainC1->Divide(2,1);
mainC0->cd();
chain->Draw("Y:X>>h0(300,-300,300,200,-200,200)","","colz");
TH2F* h0 = (TH2F*)gDirectory->FindObjectAny("h0");
h0->GetXaxis()->SetTitle("X (mm)");
h0->GetYaxis()->SetTitle("Y (mm)");
h0->SetTitle("");
mainC1->cd(1);
chain->Draw("E_ThickSi:E_CsI>>h1(1000,0,600,1000,0,30)","","colz");
TH2F* h1 = (TH2F*)gDirectory->FindObjectAny("h1");
h1->GetXaxis()->SetTitle("E_{CsI} (MeV)");
h1->GetYaxis()->SetTitle("E_{Si} (MeV)");
h1->SetTitle("");
mainC1->cd(2);
chain->Draw("E_ThinSi:E_ThickSi>>h12(1000,0,30,1000,0,10)","","colz");
TH2F* h12 = (TH2F*)gDirectory->FindObjectAny("h12");
h12->GetXaxis()->SetTitle("E_{Si} (MeV)");
h12->GetYaxis()->SetTitle("#Delta E (MeV)");
h12->SetTitle("");
mainC2->cd();
chain->Draw("ELab:ThetaLab>>h2(1000,0,50,1000,0,50)","","colz");
h2->SetMinimum(1);
TH2F* h2 = (TH2F*)gDirectory->FindObjectAny("h2");
h2->GetXaxis()->SetTitle("#theta_{lab} (deg)");
h2->GetYaxis()->SetTitle("E (MeV)");
h2->SetTitle("");
NPL::Reaction *r = new NPL::Reaction("11Be(d,3He)10Li@770");
//NPL::Reaction *r = new NPL::Reaction("34Ar(p,d)33Ar@2380");
//NPL::Reaction *r = new NPL::Reaction("46Ar(p,d)45Ar@3220");
kin = r->GetKinematicLine3();
kin->SetLineColor(2);
kin->Draw("plsame");
mainC3->cd();
chain->Draw("PhiLab:ThetaLab>>h3(100,0,50,720,-180,180)","","colz");
TH2F* h3 = (TH2F*)gDirectory->FindObjectAny("h3");
h3->GetXaxis()->SetTitle("#theta_{lab} (deg)");
h3->GetYaxis()->SetTitle("#phi_{lab} (MeV)");
h3->SetTitle("");
mainC4->cd();
chain->Draw("fIC_Incident_Position_Y:fIC_Incident_Position_X>>h4(1000,-10,10,100,-10,10)","","colz");
TH2F* h4 = (TH2F*)gDirectory->FindObjectAny("h4");
h4->GetXaxis()->SetTitle("X_{beam} (mm)");
h4->GetYaxis()->SetTitle("Y_{beam} (mm)");
h4->SetTitle("");
mainC5->cd();
chain->Draw("ThetaCM>>h5(25,0,100)",CondTot,"E1");
TH1F* h5 = (TH1F*)gDirectory->FindObjectAny("h5");
h5->GetXaxis()->SetTitle("#theta_{CM} (deg)");
h5->SetTitle("");
//Efficiency calculation
mainC6->cd();
TH1F *Efficiency = new TH1F("Efficiency", "Efficiency", 25, 0, 100);
Efficiency->SetXTitle("#theta_{CM}");
Efficiency->SetTitle("Efficiency");
Efficiency->Add(h5,4*TMath::Pi()/nentries);
Efficiency->Draw();
// Cross-Section
mainC7->cd();
mainC7->SetLogy();
TH1F *h7 = new TH1F("h7", "h7", 25, 0, 100);
h7->SetXTitle("#theta_{CM}");
h7->SetYTitle("d#sigma/d#Omega (mb/sr)");
h7->SetTitle("");
h7->Divide(h5,Efficiency,1,nentries/(4*TMath::Pi()));//262000
h7->Draw();
return;
}
//////////////////////////////////////////////////
void InitChain()
{
if(chain!=NULL){delete chain;}
chain = new TChain("AnalysedTree","");
//chain->Add("/Users/pierremorfouace/Physics/NPTool/nptool/Outputs/Analysis/34Ar_pd.root");
//chain->Add("/Users/pierremorfouace/Physics/NPTool/nptool/Outputs/Analysis/46Ar_pd.root");
chain->Add("/Users/pierremorfouace/Physics/NPTool/nptool/Outputs/Analysis/11Be_d3He_1mg.root");
return;
}
///////////////////////////////////////////////////////////////////////////////
void LoadCut()
{
TString Cut_Name = "CUT/cut_kine.root";
TString Object_Name = "cut_kine";
TFile* f_cut_kine = new TFile(Cut_Name,"read");
cut_kine = (TCutG*) f_cut_kine->FindObjectAny(Object_Name);
return;
}
/*****************************************************************************
* 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(){
M2= (TMust2Physics*) m_DetectorManager->GetDetector("MUST2Array");
SSSD= (TSSSDPhysics*) m_DetectorManager->GetDetector("SSSD");
Initial=new TInitialConditions();
InitOutputBranch();
InitInputBranch();
Rand = TRandom3();
He10Reaction= new NPL::Reaction();
He10Reaction->ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
DetectorNumber = 0 ;
ThetaNormalTarget = 0 ;
ThetaM2Surface = 0;
X_M2 = 0 ;
Y_M2 = 0 ;
Z_M2 = 0 ;
Si_E_M2 = 0 ;
CsI_E_M2 = 0 ;
E_SSSD = 0 ;
Energy = 0;
E_M2 = 0;
Si_X_M2 = 0;
Si_Y_M2 = 0;
ZTarget = 0;
TargetThickness = m_DetectorManager->GetTargetThickness()*micrometer;
// Energy loss table: the G4Table are generated by the simulation
He3CD2 = EnergyLoss("He3_CD2.G4table","G4Table",100 );
He3Al = EnergyLoss("He3_Al.G4table","G4Table",10);
He3Si = EnergyLoss("He3_Si.G4table","G4Table",10);
//Li11CD2 = EnergyLoss("Li11[0.0]_CD2.G4table","G4Table",100);
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::TreatEvent(){
// Reinitiate calculated variable
ReInitValue();
// Get the Init information on beam position and energy
// and apply by hand the experimental resolution
// This is because the beam diagnosis are not simulated
// PPAC position resolution on target is assumed to be 1mm
double XTarget = Rand.Gaus(Initial->GetIncidentPositionX(),1);
double YTarget = Rand.Gaus(Initial->GetIncidentPositionY(),1);
TVector3 BeamDirection = Initial->GetBeamDirection();
// Beam energy is measured using F3 and F2 plastic TOF
double BeamEnergy = Rand.Gaus(Initial->GetIncidentInitialKineticEnergy(),4.5);
//BeamEnergy = Li11CD2.Slow(BeamEnergy,TargetThickness/2.,0);
He10Reaction->SetBeamEnergy(BeamEnergy);
//////////////////////////// LOOP on MUST2 + SSSD Hit //////////////////
for(unsigned int countSSSD = 0 ; countSSSD < SSSD->Energy.size() ; countSSSD++){
for(unsigned int countMust2 = 0 ; countMust2 < M2->Si_E.size() ; countMust2++){
/************************************************/
//Part 0 : Get the usefull Data
// MUST2
int X = M2->Si_X[countMust2];
int Y = M2->Si_Y[countMust2];
int TelescopeNumber = M2->TelescopeNumber[countMust2];
Si_X_M2 = X ;
Si_Y_M2 = Y ;
//SSSD
int SiNumber = SSSD->DetectorNumber[countSSSD];
/************************************************/
// Matching between Thin Si and MUST2, and Forward Telescope Only
if(TelescopeNumber==SiNumber && TelescopeNumber<5){
DetectorNumber = TelescopeNumber ;
/************************************************/
// Part 1 : Impact Angle
ThetaM2Surface = 0;
ThetaNormalTarget = 0;
if(XTarget>-1000 && YTarget>-1000){
TVector3 BeamImpact(XTarget,YTarget,0);
TVector3 HitDirection = M2 -> GetPositionOfInteraction(countMust2) - BeamImpact ;
ThetaLab = HitDirection.Angle( BeamDirection );
ThetaM2Surface = HitDirection.Angle(- M2 -> GetTelescopeNormal(countMust2) );
ThetaNormalTarget = HitDirection.Angle( TVector3(0,0,1) ) ;
X_M2 = M2 -> GetPositionOfInteraction(countMust2).X() ;
Y_M2 = M2 -> GetPositionOfInteraction(countMust2).Y() ;
Z_M2 = M2 -> GetPositionOfInteraction(countMust2).Z() ;
}
else{
BeamDirection = TVector3(-1000,-1000,-1000);
ThetaM2Surface = -1000 ;
ThetaNormalTarget = -1000 ;
}
/************************************************/
/************************************************/
// Part 2 : Impact Energy
Energy = ELab = 0;
Si_E_M2 = M2->Si_E[countMust2];
CsI_E_M2= M2->CsI_E[countMust2];
E_SSSD = SSSD->Energy[countSSSD];
// if CsI
if(CsI_E_M2>0 ){
// The energy in CsI is calculate form dE/dx Table because
// 20um resolution is poor
Energy =
He3Si.EvaluateEnergyFromDeltaE(Si_E_M2,300*micrometer,
ThetaM2Surface, 0.01*MeV,
450.*MeV,0.001*MeV ,1000);
E_M2=CsI_E_M2;
}
else
Energy = Si_E_M2;
E_M2 += Si_E_M2;
// Evaluate energy using the thickness
ELab = He3Al.EvaluateInitialEnergy( Energy ,2*0.4*micrometer , ThetaM2Surface);
ELab = He3Si.EvaluateInitialEnergy( ELab ,20*micrometer , ThetaM2Surface);
ELab = He3Al.EvaluateInitialEnergy( ELab ,0.4*micrometer , ThetaM2Surface);
// Target Correction
ELab = He3CD2.EvaluateInitialEnergy( ELab ,TargetThickness/2., ThetaNormalTarget);
/************************************************/
/************************************************/
// Part 3 : Excitation Energy Calculation
Ex = He10Reaction -> ReconstructRelativistic( ELab , ThetaLab );
/************************************************/
/************************************************/
// Part 4 : Theta CM Calculation
ThetaCM = He10Reaction -> EnergyLabToThetaCM( ELab , ThetaLab)/deg;
ThetaLab=ThetaLab/deg;
/************************************************/
}
} //end loop SSSD
}//end loop MUST2
}
////////////////////////////////////////////////////////////////////////////////
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");
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::InitInputBranch(){
RootInput:: getInstance()->GetChain()->SetBranchStatus("InitialConditions",true );
RootInput:: getInstance()->GetChain()->SetBranchStatus("fIC_*",true );
RootInput:: getInstance()->GetChain()->SetBranchAddress("InitialConditions",&Initial);
}
////////////////////////////////////////////////////////////////////////////////
void Analysis::ReInitValue(){
Ex = -1000 ;
ELab = -1000;
ThetaLab = -1000;
ThetaCM = -1000;
}
////////////////////////////////////////////////////////////////////////////////
// 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;
}
#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"THiraPhysics.h"
#include "TMust2Physics.h"
#include "TSSSDPhysics.h"
#include "TInitialConditions.h"
#include "NPEnergyLoss.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 Ex;
double ELab;
double ThetaLab;
double ThetaCM;
NPL::Reaction* He10Reaction;
// intermediate variable
TRandom3 Rand;
int DetectorNumber;
double ThetaNormalTarget;
double ThetaM2Surface;
double X_M2;
double Y_M2;
double Z_M2;
double Si_E_M2;
double CsI_E_M2;
double E_SSSD;
double Energy = 0;
double E_M2 = 0;
double Si_X_M2 = 0;
double Si_Y_M2 = 0;
double ZTarget = 0;
double TargetThickness;
NPL::EnergyLoss He3CD2 ;
NPL::EnergyLoss He3Al ;
NPL::EnergyLoss He3Si ;
NPL::EnergyLoss Li11CD2 ;
TMust2Physics* M2;
TSSSDPhysics* SSSD;
TInitialConditions* Initial;
};
#endif
cmake_minimum_required (VERSION 2.8)
include("../../NPLib/ressources/CMake/NPAnalysis.cmake")
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