Skip to content
Snippets Groups Projects
Commit 3ca6f6c0 authored by deserevi's avatar deserevi
Browse files

* Fix bug in TW1Physics.cxx

   + GetPositionOfInteraction() method is now implemented

* Add analysis support to W1
   + To be tested
parent 3c98d8a3
No related branches found
No related tags found
No related merge requests found
Analyse:
make -C ./src
clean:
make clean -C ./src
distclean:
make clean -C ./src
rm Analysis
TTreeName
SimulatedTree
RootFileName
../../Outputs/Simulation/myResult.root
% ../../Outputs/Simulation/134Snpt_1h9_10MeVA_T0_B0_E0_S2mm.root
% ../../Outputs/Simulation/132Sndp_3p3_10MeVA_T0_B1_E0_S05mm.root
% ../../Outputs/Simulation/134Snpt_1h9_10MeVA_T1_B1_E0_S05mm.root
// 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 "DetectorManager.h"
#include "NPOptionManager.h"
#include "TW1Physics.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 <TRandom3.h>
// NPL
#include "NPReaction.h"
#include "RootInput.h"
#include "RootOutput.h"
#include "TInteractionCoordinates.h"
#include "TInitialConditions.h"
// Use CLHEP System of unit and Physical Constant
#include "CLHEP/Units/GlobalSystemOfUnits.h"
#include "CLHEP/Units/PhysicalConstants.h"
// ----------------------------------------------------------------------------------------------
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
{
// Declare your Energy loss here
// EnergyLoss LightTargetCD2 = EnergyLoss("proton_cd2.txt", 100, 1, 1); // LISE++
// For 132Sn(d,p)
// CD2
EnergyLoss LightTarget = EnergyLoss("proton_CD2.G4table", "G4Table", 1000); // G4
EnergyLoss BeamTarget = EnergyLoss("Sn132[0.0]_CD2.G4table", "G4Table", 1000); // G4
// solid D2
// EnergyLoss LightTarget = EnergyLoss("proton_D2_solid.G4table", "G4Table", 1000); // G4
// EnergyLoss BeamTarget = EnergyLoss("Sn132[0.0]_D2_solid.G4table", "G4Table", 1000); // G4
// For 134Sn(p,t)
// CH2
// EnergyLoss LightTarget = EnergyLoss("triton_CH2.G4table", "G4Table", 1000); // G4
// EnergyLoss BeamTarget = EnergyLoss("Sn134[0.0]_CH2.G4table", "G4Table", 1000); // G4
// solid H2
// EnergyLoss LightTarget = EnergyLoss("triton_H2_solid.G4table", "G4Table", 1000); // G4
// EnergyLoss BeamTarget = EnergyLoss("Sn134[0.0]_H2_solid.G4table", "G4Table", 1000); // G4
// For 132Sn(d,t)
// EnergyLoss LightTarget = EnergyLoss("triton_CD2.G4table", "G4Table", 1000); // G4
// EnergyLoss BeamTarget = EnergyLoss("Sn132[0.0]_CD2.G4table", "G4Table", 1000); // G4
}
using namespace ENERGYLOSS ;
// ----------------------------------------------------------------------------------------------
/////////////////////////////////////////////////////////////////////////////////////////////////
void DisplayInputCrossSection()
{
// Path to cross-section files
TString path = gSystem->Getenv("NPTOOL");
path += "/Inputs/CrossSection/";
// Read cross-sections 132Sn(d,p)
/* TGraph *gr1 = new TGraph(path + "132Sndp_5A_MeV_3p3_ZR_FRC.lis");
TGraph *gr2 = new TGraph(path + "132Sndp_5A_MeV_2f7_ZR_FRC.lis");
TGraph *gr3 = new TGraph(path + "132Sndp_10A_MeV_3p3_ZR_FRC.lis");
TGraph *gr4 = new TGraph(path + "132Sndp_10A_MeV_2f7_ZR_FRC.lis");*/
// Jacques
TGraph *gr1 = ReadCrossSection("132Sndp_5A_MeV_3p3_ZR_FRC.lis");
TGraph *gr2 = ReadCrossSection("132Sndp_5A_MeV_2f7_ZR_FRC.lis");
TGraph *gr3 = ReadCrossSection("132Sndp_10A_MeV_3p3_ZR_FRC.lis");
TGraph *gr4 = ReadCrossSection("132Sndp_10A_MeV_2f7_ZR_FRC.lis");
// Didier
TGraph *gr5 = new TGraph(path + "sn132dp_gs_10AMeV.txt");
// Read cross-section 134Sn(p,t)
// Didier
TGraph *grpt1 = ReadCrossSectionPT("CS_Ep10MeV_sn134pt_gs_1h9demi.dat");
TGraph *grpt2 = ReadCrossSectionPT("CS_Ep15MeV_sn134pt_gs_1h9demi.dat");
TGraph *grpt3 = ReadCrossSectionPT("CS_Ep20MeV_sn134pt_gs_1h9demi.dat");
// Read cross-section 132Sn(d,d)
// Angel
TGraph *grdd = new TGraph(path + "132Sndd_10A_MeV_ruth.dat");
TGraph *grpp = new TGraph(path + "132Snpp_10A_MeV_ruth.dat");
// Draw cross-sections
TCanvas *can = new TCanvas("can");
can->SetLogy();
can->Draw();
// TH2F *hframe = new TH2F("hframe", "^{2}H(^{132}Sn,p)^{133}Sn", 180, 0, 180, 100, 1e-2, 100);
// TH2F *hframe = new TH2F("hframe", "", 180, 0, 180, 100, 1e-2, 100);
// TH2F *hframe = new TH2F("hframe", "^{1}H(^{134}Sn,t)^{132}Sn_{g.s.}", 180, 0, 180, 100, 1e-8, 1e-5);
// TH2F *hframe = new TH2F("hframe", "", 180, 0, 180, 100, 1e-8, 1e-5);
TH2F *hframe = new TH2F("hframe", "", 180, 0, 180, 100, 1e-3, 10);
hframe->Draw();
hframe->SetXTitle("#Theta_{c.m.} [deg]");
// hframe->SetYTitle("d#sigma/d#Omega [mb/sr]");
hframe->SetYTitle("d#sigma/d#Omega / (d#sigma/d#Omega)_{R}");
// hframe->SetYTitle("#propto d#sigma/d#Omega");
/* grpt1->SetLineColor(kRed); grpt1->Draw("l");
grpt2->SetLineColor(kMagenta); grpt2->Draw("l");
grpt3->SetLineColor(kBlue); grpt3->Draw("l");*/
/* gr1->SetLineColor(kRed); gr1->SetLineStyle(2); gr1->Draw("l");
gr2->SetLineColor(kRed); gr2->Draw("l");
gr3->SetLineColor(kBlue); gr3->SetLineStyle(2); gr3->Draw("l");
gr4->SetLineColor(kBlue); gr4->Draw("l");*/
// gr5->Draw("l");
// grdd->Draw("l");
grpp->Draw("l");
// TLegend
TLegend *leg = new TLegend(0.50, 0.64, 0.82, 0.84);
// leg->AddEntry(grdd, "10 MeV/u", "l");
leg->AddEntry(grpp, "10 MeV/u", "l");
/* TLegend *leg = new TLegend(0.50, 0.64, 0.82, 0.84);
leg->AddEntry(grpt1, "1h9/2 10 MeV/u", "l");
leg->AddEntry(grpt2, "1h9/2 15 MeV/u", "l");
leg->AddEntry(grpt3, "1h9/2 20 MeV/u", "l");*/
/* TLegend *leg = new TLegend(0.16, 0.17, 0.48, 0.37);
leg->AddEntry(gr1, "3p3/2 5 MeV/u", "l");
leg->AddEntry(gr2, "2f7/2 5 MeV/u", "l");
leg->AddEntry(gr3, "3p3/2 10 MeV/u", "l");
leg->AddEntry(gr4, "2f7/2 10 MeV/u", "l");*/
leg->SetBorderSize(1);
leg->Draw();
/* TMultiGraph *mgr = new TMultiGraph();
mgr->Add(gr1, "lp");
mgr->Add(gr2, "lp");
mgr->Add(gr3, "lp");
mgr->Add(gr4, "lp");
mgr->Draw("a*");*/
// gr1->Draw("alp");
}
TGraph* ReadCrossSection(const char* fname)
{
// Path to cross-section files
TString path = gSystem->Getenv("NPTOOL");
path += "/Inputs/CrossSection/";
// Open file
ifstream fich;
fich.open(path + fname);
if (!fich) cout << "Probleme d'ouverture dans le fichier " << fname << endl;
// Read file
Double_t angle, sigma;
Int_t nlines = 0;
TGraph *gr = new TGraph();
while (fich >> angle >> sigma) {
gr->SetPoint(nlines++, angle, sigma * 15); // 15: fm^2 -> mb + D0^2
}
// Close file
fich.close();
// TGraph name
gr->SetTitle(fname);
// test pour savoir si on a bien rempli le TGraph
if (gr->GetN() == 0)
cout << "Mauvaise lecture du fichier --> probablement mauvais format" << endl;
return gr;
}
TGraph* ReadCrossSectionPT(const char* fname)
{
// Path to cross-section files
TString path = gSystem->Getenv("NPTOOL");
path += "/Inputs/CrossSection/";
// Open file
ifstream fich;
fich.open(path + fname);
if (!fich) cout << "Probleme d'ouverture dans le fichier " << fname << endl;
// Read file
Double_t angle, sigma, dum;
Int_t nlines = 0;
TGraph *gr = new TGraph();
// while (fich >> angle >> sigma >> dum >> dum >> dum >> dum >> dum >> dum >> dum >> dum >> dum) {
while (fich >> angle >> sigma) {
gr->SetPoint(nlines++, angle, sigma);
}
// Close file
fich.close();
// TGraph name
gr->SetTitle(fname);
// test pour savoir si on a bien rempli le TGraph
if (gr->GetN() == 0)
cout << "Mauvaise lecture du fichier --> probablement mauvais format" << endl;
return gr;
}
#include "ObjectManager.hh"
using namespace std;
int main(int argc,char** argv)
{
// command line parsing
NPOptionManager* myOptionManager = NPOptionManager::getInstance(argc,argv);
string reactionfileName = myOptionManager->GetReactionFilePath();
string detectorfileName = myOptionManager->GetDetectorFilePath();
string calibrationfileName = myOptionManager->GetCalibrationFilePath();
string runToReadfileName = myOptionManager->GetRunToReadFilePath();
string OutputfileName = myOptionManager->GetOutputFilePath();
// Instantiate RootInput and RootOutput singleton classes
RootInput:: getInstance(runToReadfileName);
// RootOutput::getInstance("Analysis/"+OutputfileName, "AnalyzedTree") ;
RootOutput::getInstance("Analysis/W1_AnalyzedData", "AnalyzedTree");
// Initialize the reaction
NPL::Reaction* myReaction = new Reaction();
myReaction->ReadConfigurationFile(reactionfileName);
// Initialize the detector
NPA::DetectorManager* myDetector = new DetectorManager;
myDetector->ReadConfigurationFile(detectorfileName);
// Calculate beam energy at target middle
// Target informations
cout << "/////////// Target information ///////////" << endl;
cout << "Thickness (um): " << myDetector->GetTargetThickness() << endl;
// Get nominal beam energy
Double_t BeamEnergyNominal = myReaction->GetBeamEnergy() * MeV;
cout << "Nominal beam energy (MeV): " << BeamEnergyNominal << endl;
// Slow beam at target middle
Double_t BeamEnergy = BeamTarget.Slow(BeamEnergyNominal, myDetector->GetTargetThickness()/2 * micrometer, 0);
cout << "Middle-target beam energy (MeV): " << BeamEnergy << endl << endl;
// Set energy beam at target middle
myReaction->SetBeamEnergy(BeamEnergy);
// Attach more branch to the output
double Ex = 0 ; double ExNoStrips = 0 ; double EE = 0 ; double TT = 0 ; double X = 0 ; double Y = 0 ; int det ;
RootOutput::getInstance()->GetTree()->Branch("ExcitationEnergy",&Ex,"Ex/D") ;
RootOutput::getInstance()->GetTree()->Branch("ExcitationEnergyNoStrips",&ExNoStrips,"ExNoStrips/D") ;
RootOutput::getInstance()->GetTree()->Branch("E",&EE,"EE/D") ;
RootOutput::getInstance()->GetTree()->Branch("A",&TT,"TT/D") ;
RootOutput::getInstance()->GetTree()->Branch("X",&X,"X/D") ;
RootOutput::getInstance()->GetTree()->Branch("Y",&Y,"Y/D") ;
// Get GaspardTracker pointer
TW1Physics* W1 = (TW1Physics*) myDetector->m_Detector["W1"];
// Get the input TChain and treat it
TChain* chain = RootInput:: getInstance() -> GetChain();
// Connect TInitialConditions branch
TInitialConditions *initCond = 0;
chain->SetBranchAddress("InitialConditions", &initCond);
chain->SetBranchStatus("InitialConditions", 1);
// Connect TInteractionCoordinates branch
TInteractionCoordinates *interCoord = 0;
chain->SetBranchAddress("InteractionCoordinates", &interCoord);
chain->SetBranchStatus("InteractionCoordinates", 0);
// Analysis is here!
int nentries = chain->GetEntries();
cout << "/////////// Loop information ///////////" << endl;
cout << "Number of entries to be analysed: " << nentries << endl;
// Default initialization
double XTarget = 0;
double YTarget = 0;
double BeamTheta = 0;
double BeamPhi = 0;
// random generator
TRandom3 *gene = new TRandom3();
// Loop on all events
for (int i = 0; i < nentries; i ++) {
if (i%10000 == 0 && i!=0) cout << "\r" << i << " analyzed events" << flush;
chain -> GetEntry(i);
// Treat Gaspard event
myDetector->ClearEventPhysics();
myDetector->BuildPhysicalEvent();
// Get total energy
double E = W1->GetEnergy(0);
// if there is a hit in the detector array, treat it.
double Theta, ThetaStrip, angle, ThetaCM;
double DetecX, DetecY, DetecZ;
double r;
TVector3 A;
if (E > -1000) {
// Get c.m. angle
ThetaCM = initCond->GetICEmittedAngleThetaCM(0) * deg;
// Get exact scattering angle from TInteractionCoordinates object
// interCoord->Dump();
// cout << i << " mult: " << interCoord->GetDetectedMultiplicity() << endl;
DetecX = interCoord->GetDetectedPositionX(0);
DetecY = interCoord->GetDetectedPositionY(0);
DetecZ = interCoord->GetDetectedPositionZ(0);
TVector3 Detec(DetecX, DetecY, DetecZ);
// Get interaction position in detector
// This takes into account the strips
A = W1->GetPositionOfInteraction(0);
// Get beam interaction coordinates on target (from initial condition)
XTarget = initCond->GetICPositionX(0);
YTarget = initCond->GetICPositionY(0);
BeamTheta = initCond->GetICIncidentAngleTheta(0)*deg;
BeamPhi = initCond->GetICIncidentAnglePhi(0)*deg;
TVector3 BeamDirection = TVector3(cos(BeamPhi)*sin(BeamTheta), sin(BeamPhi)*sin(BeamTheta), cos(BeamTheta));
// Hit direction taking into account beam position on target
TVector3 HitDirection = A - TVector3(XTarget, YTarget, 0);
// cout << "A: " << A.X() << " " << A.Y() << " " << A.Z() << endl;
// cout << "HitDirection: " << HitDirection.X() << " " << HitDirection.Y() << " " << HitDirection.Z() << endl;
// Calculate scattering angle w.r.t. optical beam axis (do not take into account beam position on target)
ThetaStrip = ThetaCalculation(A, TVector3(0,0,1));
Theta = ThetaCalculation(Detec, TVector3(0, 0, 1));
// Calculate scattering angle w.r.t. beam (ideal case)
// ThetaStrip = ThetaCalculation(HitDirection, BeamDirection);
// Theta = ThetaCalculation(Detec - TVector3(XTarget, YTarget, 0), BeamDirection);
// Calculate scattering angle w.r.t. beam (finite spatial resolution)
// double resol = 800; // in micrometer
// angle = gene->Rndm() * 2*3.14;
// r = fabs(gene->Gaus(0, resol)) * micrometer;
// ThetaStrip = ThetaCalculation(A - TVector3(XTarget + r*cos(angle), YTarget + r*sin(angle), 0), BeamDirection);
// Theta = ThetaCalculation(Detec - TVector3(XTarget + r*cos(angle), YTarget + r*sin(angle), 0), BeamDirection);
//
// Correct for energy loss in the target
E = LightTarget.EvaluateInitialEnergy(E, myDetector->GetTargetThickness()/2 * micrometer, ThetaStrip);
// Calculate excitation energy
// if (Theta/deg > 150 && Theta/deg < 180) {
// if (Theta/deg < 60 && ThetaCM/deg < 90) {
// if (Theta/deg > 35 && Theta/deg < 45 && E/MeV < 17) {
// if (Theta/deg < 45) {
// if (E/MeV < 38) { // for (p,t) reaction
if (Theta/deg > 30) { // for (d,p) reaction
ExNoStrips = myReaction->ReconstructRelativistic(E, Theta / rad);
Ex = myReaction->ReconstructRelativistic(E, ThetaStrip);
}
else {
Ex = -200;
ExNoStrips = -200;
}
}
else {
Ex = -100;
ExNoStrips = -100;
}
EE = E ; TT = ThetaStrip/deg;
if (E>-1000) {
X = A . X();
Y = A . Y();
}
else {
X = -1000 ; Y = -1000;
}
// Fill output tree
RootOutput::getInstance()->GetTree()->Fill();
}
// delete singleton classes
RootOutput::getInstance()->Destroy();
RootInput::getInstance()->Destroy();
NPOptionManager::getInstance()->Destroy();
return 0;
}
double ThetaCalculation (TVector3 A , TVector3 B)
{
double Theta = acos( (A.Dot(B)) / (A.Mag()*B.Mag()) );
return Theta ;
}
###Make file made by Adrien MATTA/ Institut de Physique Nucleaire d'Orsay IPNO###
# Made to compile the ROOT Analyser for MUST2 experiment
CPP=g++
EXEC=Analysis
# local includes
NPAINCLUDES = ../include
# ROOT includes
CXXFLAGS += `root-config --cflags`
# CLHEP includes
CXXFLAGS += -I$(CLHEP_INCLUDE_DIR)
CXXFLAGS += -I$(NPAINCLUDES)
CXXFLAGS += -I$(NPLIB)/include
LDFLAGS = `root-config --libs` -lMathMore
LDFLAGS+= `$(NPLIB)/liblist`
LDFLAGS+= -L$(CLHEP_LIB_DIR) -l$(CLHEP_LIB)
SRC= $(wildcard *.cc)
INC= $(wildcard $(NPAINCLUDES)/*.hh)
OBJ=$(SRC:.cc=.o)
#all:$(EXEC)
# @$(CPP) -o $@ -c $< $(CXXFLAGS)
Analysis: Analysis.o $(INC)
@$(CPP) $(CXXFLAGS) -o $@ $^ $(LDFLAGS)
mv Analysis ../Analysis
%.o: %.cc
@$(CPP) $(CXXFLAGS) -o $@ -c $<
.PHONY: clean mrproper
clean:
rm -rf *.o
mrproper: clean
rm -rf $(EXEC)
...@@ -445,6 +445,17 @@ void TW1Physics::AddDetector(double theta, double phi, double distance, ...@@ -445,6 +445,17 @@ void TW1Physics::AddDetector(double theta, double phi, double distance,
TVector3 TW1Physics::GetPositionOfInteraction(int i)
{
TVector3 Position = TVector3(GetStripPositionX(fDetectorNumber[i], fFrontStrip[i], fBackStrip[i]),
GetStripPositionY(fDetectorNumber[i], fFrontStrip[i], fBackStrip[i]),
GetStripPositionZ(fDetectorNumber[i], fFrontStrip[i], fBackStrip[i]));
return Position;
}
TVector3 TW1Physics::GetDetectorNormal(int i) TVector3 TW1Physics::GetDetectorNormal(int i)
{ {
TVector3 U = TVector3(GetStripPositionX(fDetectorNumber[i], m_NumberOfStrips, 1), TVector3 U = TVector3(GetStripPositionX(fDetectorNumber[i], m_NumberOfStrips, 1),
......
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