From 0768500eee200378ab8565d0cde8932d91bf7322 Mon Sep 17 00:00:00 2001
From: adrien matta <matta@lpccaen.in2p3.fr>
Date: Wed, 29 Jul 2020 17:14:13 +0200
Subject: [PATCH] * progress on strasse analysis
---
NPLib/Physics/TReactionConditions.cxx | 1 -
NPLib/Physics/TReactionConditions.h | 10 ++-
Projects/Strasse/Analysis.cxx | 88 ++++++++++++++++---
Projects/Strasse/Analysis.h | 13 ++-
.../geometry/strasse_optimized.detector | 2 +-
Projects/Strasse/macro/Shift.cxx | 62 +++++++++++--
Projects/Strasse/reaction/C12_p2p.reaction | 1 -
7 files changed, 151 insertions(+), 26 deletions(-)
diff --git a/NPLib/Physics/TReactionConditions.cxx b/NPLib/Physics/TReactionConditions.cxx
index f4400ba39..460075917 100644
--- a/NPLib/Physics/TReactionConditions.cxx
+++ b/NPLib/Physics/TReactionConditions.cxx
@@ -59,7 +59,6 @@ void TReactionConditions::Clear(){
fRC_Momentum_Direction_X.clear();
fRC_Momentum_Direction_Y.clear();
fRC_Momentum_Direction_Z.clear();
- fRC_Momentum.clear();
}
////////////////////////////////////////////////////////////////////////////////
void TReactionConditions::Dump() const{
diff --git a/NPLib/Physics/TReactionConditions.h b/NPLib/Physics/TReactionConditions.h
index 93c70537e..1a85cda8c 100644
--- a/NPLib/Physics/TReactionConditions.h
+++ b/NPLib/Physics/TReactionConditions.h
@@ -53,6 +53,7 @@ private:
double fRC_Beam_Emittance_Theta;
double fRC_Beam_Emittance_Phi;
double fRC_Beam_Reaction_Energy;
+
double fRC_Vertex_Position_X;
double fRC_Vertex_Position_Y;
double fRC_Vertex_Position_Z;
@@ -68,8 +69,6 @@ private:
vector<double> fRC_Momentum_Direction_X;
vector<double> fRC_Momentum_Direction_Y;
vector<double> fRC_Momentum_Direction_Z;
- vector<TVector3> fRC_Momentum;
-
public:
TReactionConditions();
@@ -82,7 +81,7 @@ public:
///////////////////// SETTERS ////////////////////////
// Beam parameter
void SetBeamParticleName (const string & Beam_Particle_Name) {fRC_Beam_Particle_Name = Beam_Particle_Name;}//!
- void SetBeamReactionEnergy (const double & Beam_Reaction_Energy) {fRC_Beam_Reaction_Energy = Beam_Reaction_Energy;}//!
+ void SetBeamReactionEnergy (const double & Beam_Reaction_Energy) {fRC_Beam_Reaction_Energy = Beam_Reaction_Energy;}//!
void SetBeamEmittanceTheta (const double & Beam_Emittance_Theta) {fRC_Beam_Emittance_Theta = Beam_Emittance_Theta;}//!
void SetBeamEmittancePhi (const double & Beam_Emittance_Phi) {fRC_Beam_Emittance_Phi = Beam_Emittance_Phi;}//!
void SetBeamEmittanceThetaX (const double & Beam_Emittance_ThetaX) {fRC_Beam_Emittance_ThetaX = Beam_Emittance_ThetaX;}//!
@@ -124,6 +123,8 @@ public:
double GetVertexPositionX () const {return fRC_Vertex_Position_X ;}//!
double GetVertexPositionY () const {return fRC_Vertex_Position_Y ;}//!
double GetVertexPositionZ () const {return fRC_Vertex_Position_Z ;}//!
+ TVector3 GetVertexPosition () const {
+ return TVector3(fRC_Vertex_Position_X,fRC_Vertex_Position_Y,fRC_Vertex_Position_Z) ;}//!
double GetExcitation3 () const {return fRC_ExcitationEnergy3 ;}//!
double GetExcitation4 () const {return fRC_ExcitationEnergy4 ;}//!
double GetThetaCM () const {return fRC_ThetaCM;}//!
@@ -138,7 +139,8 @@ public:
double GetMomentumDirectionX (const int &i) const {return fRC_Momentum_Direction_X[i];}//!
double GetMomentumDirectionY (const int &i) const {return fRC_Momentum_Direction_Y[i];}//!
double GetMomentumDirectionZ (const int &i) const {return fRC_Momentum_Direction_Z[i];}//!
- TVector3 GetParticleMomentum (const int &i) const {return TVector3(fRC_Momentum_Direction_X[i],fRC_Momentum_Direction_Y[i],fRC_Momentum_Direction_Z[i]).Unit();}//!
+ TVector3 GetParticleMomentum (const int &i) const {
+ return TVector3(fRC_Momentum_Direction_X[i],fRC_Momentum_Direction_Y[i],fRC_Momentum_Direction_Z[i]).Unit();}//!
TVector3 GetBeamDirection () const ;
TVector3 GetParticleDirection (const int i) const ;
diff --git a/Projects/Strasse/Analysis.cxx b/Projects/Strasse/Analysis.cxx
index 1f5acb933..442d6c0d4 100644
--- a/Projects/Strasse/Analysis.cxx
+++ b/Projects/Strasse/Analysis.cxx
@@ -75,10 +75,12 @@ void Analysis::TreatEvent(){
TVector3 InnerPos1 = Strasse->GetInnerPositionOfInteraction(0);
TVector3 OuterPos1 = Strasse->GetOuterPositionOfInteraction(0);
TVector3 Proton1 = OuterPos1-InnerPos1;
+ Proton1=Proton1.Unit();
// Proton 2
TVector3 InnerPos2 = Strasse->GetInnerPositionOfInteraction(1);
TVector3 OuterPos2 = Strasse->GetOuterPositionOfInteraction(1);
TVector3 Proton2 = OuterPos2-InnerPos2;
+ Proton2=Proton2.Unit();
deltaPhi = abs(Proton1.Phi()/deg-Proton2.Phi()/deg);
sumTheta = Proton1.Theta()/deg+Proton2.Theta()/deg;
@@ -105,34 +107,75 @@ void Analysis::TreatEvent(){
unsigned int i1,i2;
i1 = Catana->FindClosestHitToLine(InnerPos1,OuterPos1,d1);
i2 = Catana->FindClosestHitToLine(InnerPos2,OuterPos2,d2);
- if(i1!=i2){
+ //if(i1!=i2){
+ if(true){
+ double TA = BeamTarget.Slow(InitialBeamEnergy,abs(VertexZ-75),RC->GetBeamDirection().Angle(TVector3(0,0,1)));
+
+ ////////////////////////////////////
+ // From Reaction Conditions
+ double E1s = RC->GetKineticEnergy(0);
+ double E2s = RC->GetKineticEnergy(1);
+ TVector3 Proton1s=RC->GetParticleMomentum(0).Unit();
+ TVector3 Proton2s=RC->GetParticleMomentum(1).Unit();
+ // Matching the right energy with the right proton
+ if((Proton1s-Proton1).Mag()<(Proton1s-Proton2).Mag()){
+ E1=E1s;
+ E2=E2s;
+ alpha=Proton1s.Angle(Proton1)/deg;
+ Theta1=Proton1.Theta();
+ Phi1=Proton1.Phi();
+ Theta2=Proton2.Theta();
+ Phi2=Proton2.Phi();
+ Theta1s=Proton1s.Theta();
+ Phi1s=Proton1s.Phi();
+ Theta2s=Proton2s.Theta();
+ Phi2s=Proton2s.Phi();
+
+ }
+ else{
+ E2=E1s;
+ E1=E2s;
+ alpha=Proton2s.Angle(Proton1)/deg;
+ Theta1=Proton1.Theta()/deg;
+ Phi1=Proton1.Phi()/deg;
+ Theta2=Proton2.Theta()/deg;
+ Phi2=Proton2.Phi()/deg;
+ Theta1s=Proton2s.Theta()/deg;
+ Phi1s=Proton2s.Phi()/deg;
+ Theta2s=Proton1s.Theta()/deg;
+ Phi2s=Proton1s.Phi()/deg;
+ }
+ // From detectors
E1 = ReconstructProtonEnergy(Vertex,Proton1,Catana->Energy[i1]);
E2 = ReconstructProtonEnergy(Vertex,Proton2,Catana->Energy[i2]);
- double TA = BeamTarget.Slow(InitialBeamEnergy,abs(VertexZ-75),0);
- // setting up Lorentz Vector from measured trajectories and energies
- TVector3 PA(0,0,sqrt(TA*(TA+2*m_QFS->GetNucleusA()->Mass()))); // for like there is no BDC
- Proton1=E1*Proton1.Unit();
- Proton2=E2*Proton2.Unit();
+
+ // Vertex = RC->GetVertexPosition();
+ //TA = RC->GetBeamEnergy();
+ ////////////////////////////////////
+ // setting up Lorentz Vector from measured trajectories and energies
+ // TVector3 PA(0,0,sqrt(TA*(TA+2*m_QFS->GetNucleusA()->Mass()))); // for like there is no BDC
+ double beam_mom=sqrt(TA*(TA+2*m_QFS->GetNucleusA()->Mass()));
+ // TVector3 PA(0,0,beam_mom); // for like there is no BDC
+ TVector3 PA=beam_mom*RC->GetBeamDirection().Unit();
+
LV_A.SetVectM(PA,m_QFS->GetNucleusA()->Mass());
double P1= sqrt(E1*(E1+2*NPUNITS::proton_mass_c2));
double P2= sqrt(E2*(E2+2*NPUNITS::proton_mass_c2));
LV_p1.SetVectM(Proton1.Unit()*P1,NPUNITS::proton_mass_c2);
LV_p2.SetVectM(Proton2.Unit()*P2,NPUNITS::proton_mass_c2);
-
// computing Ex from Missing Mass
LV_B = LV_A + LV_T - LV_p1 - LV_p2;
//LV_B = RC->GetParticleMomentum(2);
Ex = LV_B.M() - m_QFS->GetNucleusB()->Mass();
}
-
}
}
////////////////////////////////////////////////////////////////////////////////
double Analysis::ReconstructProtonEnergy(const TVector3& x0, const TVector3& dir,const double& Ecatana){
- TVector3 Normal = TVector3(0,0,1);
+ TVector3 Normal = TVector3(0,1,0);
Normal.SetPhi(dir.Phi());
double Theta = dir.Angle(Normal);
// Catana Al housing
@@ -146,12 +189,11 @@ double Analysis::ReconstructProtonEnergy(const TVector3& x0, const TVector3& dir
// LH2 target
static TVector3 x1;
x1= x0+dir;
- TVector3 T1(0,30,0);
- TVector3 T2(0,30,1);
+ TVector3 T1(0,15,0);
+ TVector3 T2(0,15,1);
T1.SetPhi(dir.Phi());
T2.SetPhi(dir.Phi());
- TVector3 Vertex,delta;
- double d = NPL::MinimumDistanceTwoLines(x0,x1,T1,T2,Vertex,delta);
+ double d = NPL::MinimumDistancePointLine(T1,T2,x0);
E = protonTarget.EvaluateInitialEnergy(E,d,Theta);
}
@@ -185,6 +227,17 @@ void Analysis::InitOutputBranch() {
RootOutput::getInstance()->GetTree()->Branch("deltaPhi",&deltaPhi,"deltaPhi/D");
RootOutput::getInstance()->GetTree()->Branch("sumTheta",&sumTheta,"sumTheta/D");
+ RootOutput::getInstance()->GetTree()->Branch("alpha",&alpha,"alpha/D");
+
+ RootOutput::getInstance()->GetTree()->Branch("Theta1",&Theta1,"Theta1/D");
+ RootOutput::getInstance()->GetTree()->Branch("Phi1",&Phi1,"Phi1/D");
+ RootOutput::getInstance()->GetTree()->Branch("Theta2",&Theta2,"Theta2/D");
+ RootOutput::getInstance()->GetTree()->Branch("Phi2",&Phi2,"Phi2/D");
+ RootOutput::getInstance()->GetTree()->Branch("Theta1s",&Theta1s,"Theta1s/D");
+ RootOutput::getInstance()->GetTree()->Branch("Phi1s",&Phi1s,"Phi1s/D");
+ RootOutput::getInstance()->GetTree()->Branch("Theta2s",&Theta2s,"Theta2s/D");
+ RootOutput::getInstance()->GetTree()->Branch("Phi2s",&Phi2s,"Phi2s/D");
+
RootOutput::getInstance()->GetTree()->Branch("Distance",&Distance,"Distance/D");
RootOutput::getInstance()->GetTree()->Branch("InteractionCoordinates","TInteractionCoordinates",&DC);
@@ -212,6 +265,15 @@ void Analysis::ReInitValue(){
Distance=-1000;
sumTheta=-1000;
deltaPhi=-1000;
+ alpha=-1000;
+ Theta1=-1000;
+ Phi1=-1000;
+ Theta2=-1000;
+ Phi2=-1000;
+ Theta1s=-1000;
+ Phi1s=-1000;
+ Theta2s=-1000;
+ Phi2s=-1000;
}
diff --git a/Projects/Strasse/Analysis.h b/Projects/Strasse/Analysis.h
index 52cf9bea4..9754678af 100644
--- a/Projects/Strasse/Analysis.h
+++ b/Projects/Strasse/Analysis.h
@@ -50,7 +50,7 @@ class Analysis: public NPL::VAnalysis{
void ReInitValue();
static NPL::VAnalysis* Construct();
TVector3 InterpolateInPlaneZ(TVector3,TVector3,double);
-
+
private:
double Ex;
double E1;
@@ -66,6 +66,17 @@ class Analysis: public NPL::VAnalysis{
double Distance;
double deltaPhi;
double sumTheta;
+ double alpha;
+ double Theta1=-1000;
+ double Phi1=-1000;
+ double Theta2=-1000;
+ double Phi2=-1000;
+ double Theta1s=-1000;
+ double Phi1s=-1000;
+ double Theta2s=-1000;
+ double Phi2s=-1000;
+
+
TLorentzVector LV_A;
TLorentzVector LV_T;
TLorentzVector LV_B;
diff --git a/Projects/Strasse/geometry/strasse_optimized.detector b/Projects/Strasse/geometry/strasse_optimized.detector
index 25d406002..f15b2d807 100644
--- a/Projects/Strasse/geometry/strasse_optimized.detector
+++ b/Projects/Strasse/geometry/strasse_optimized.detector
@@ -3,7 +3,7 @@ Target
THICKNESS= 150 mm
ANGLE= 0 deg
RADIUS= 15 mm
- MATERIAL= LH2
+ MATERIAL= LH2
X= 0 mm
Y= 0 mm
Z= 0 mm
diff --git a/Projects/Strasse/macro/Shift.cxx b/Projects/Strasse/macro/Shift.cxx
index 235948643..6a593fa26 100644
--- a/Projects/Strasse/macro/Shift.cxx
+++ b/Projects/Strasse/macro/Shift.cxx
@@ -1,21 +1,73 @@
void Shift(){
+
TFile* file_ok = TFile::Open("../../Outputs/Analysis/strasse_ok.root");
TFile* file_shifted = TFile::Open("../../Outputs/Analysis/strasse_shifted.root");
TTree* ok= (TTree*) file_ok->FindObjectAny("PhysicsTree");
TTree* shifted= (TTree*) file_shifted->FindObjectAny("PhysicsTree");
+
+ string cond = "sqrt(VertexX*VertexX+VertexY*VertexY)<20 && VertexZ>-80 && VertexZ<80";
+ // Vertex
+ TCanvas* cVertex = new TCanvas("ControlVertex","ControlVertex",1000,1000);
+ cVertex->Divide(2,2);
+ cond = "VertexX!=-1000";
+ cVertex->cd(1);
+ ok->Draw("VertexX-fRC_Vertex_Position_X>>hxok(500,-2,2)",cond.c_str(),"") ;
+ shifted->Draw("VertexX-fRC_Vertex_Position_X>>hx(500,-2,2)",cond.c_str(),"same") ;
+ TH1* hx= (TH1*) gDirectory->FindObjectAny("hx");
+ hx->SetLineColor(kOrange-3);
+ hx->SetFillColor(kOrange-3);
+
+ TH1* hxok= (TH1*) gDirectory->FindObjectAny("hxok");
+ hxok->GetXaxis()->SetTitle("X_{reconstructed}-X_{real} (mm)");
+
+ cVertex->cd(2);
+
+ ok->Draw("VertexY-fRC_Vertex_Position_Y>>hyok(500,-2,2)",cond.c_str(),"") ;
+ shifted->Draw("VertexY-fRC_Vertex_Position_Y>>hy(500,-2,2)",cond.c_str(),"same") ;
+
+ TH1* hy= (TH1*) gDirectory->FindObjectAny("hy");
+ hy->SetLineColor(kOrange-3);
+ hy->SetFillColor(kOrange-3);
+ TH1* hyok= (TH1*) gDirectory->FindObjectAny("hyok");
+ hyok->GetXaxis()->SetTitle("Y_{reconstructed}-Y_{real} (mm)");
+
+
+ cVertex->cd(3);
+
+ ok->Draw("VertexZ-fRC_Vertex_Position_Z>>hzok(2000,-20,20)",cond.c_str(),"") ;
+ shifted->Draw("VertexZ-fRC_Vertex_Position_Z>>hz(2000,-20,20)",cond.c_str(),"same") ;
+
+ TH1* hz= (TH1*) gDirectory->FindObjectAny("hz");
+ hz->SetLineColor(kOrange-3);
+ hz->SetFillColor(kOrange-3);
+ TH1* hzok= (TH1*) gDirectory->FindObjectAny("hzok");
+ hzok->GetXaxis()->SetTitle("Z_{reconstructed}-Z_{real} (mm)");
+
+
+ cVertex->cd(4);
+ ok->Draw("Distance>>hdok(500,0,2)", cond.c_str(),"");
+ shifted->Draw("Distance>>hd(500,0,2)", cond.c_str(),"same");
+ TH1* hd= (TH1*) gDirectory->FindObjectAny("hd");
+ hd->SetLineColor(kOrange-3);
+ hd->SetFillColor(kOrange-3);
+ TH1* hdok= (TH1*) gDirectory->FindObjectAny("hdok");
+ hdok->GetXaxis()->SetTitle("Minimum track distance (mm)");
+
+
+
+ // Angle
TCanvas* ctheta= new TCanvas("ControlTheta","ControlTheta",2000,1000);
ctheta->Divide(2,1);
ctheta->cd(1);
- string cond = "Theta12!=-1000";
- ok->Draw("Theta12>>ht(5000)",cond.c_str(),"") ;
- shifted->Draw("Theta12>>hts(5000)",cond.c_str(),"same") ;
+ ok->Draw("Theta12>>ht(2000)",cond.c_str(),"") ;
+ shifted->Draw("Theta12>>hts(2000)",cond.c_str(),"same") ;
TH1* hts= (TH1*) gDirectory->FindObjectAny("hts");
hts->SetFillColor(kOrange-3);
hts->SetLineColor(kOrange-3);
ctheta->cd(2);
cond = "deltaPhi!=-1000";
- ok->Draw("deltaPhi>>hp(5000)",cond.c_str(),"") ;
- shifted->Draw("deltaPhi>>hps(5000)",cond.c_str(),"same") ;
+ ok->Draw("deltaPhi>>hp(2000)",cond.c_str(),"") ;
+ shifted->Draw("deltaPhi>>hps(2000)",cond.c_str(),"same") ;
TH1* hps= (TH1*) gDirectory->FindObjectAny("hps");
hps->SetFillColor(kOrange-3);
hps->SetLineColor(kOrange-3);
diff --git a/Projects/Strasse/reaction/C12_p2p.reaction b/Projects/Strasse/reaction/C12_p2p.reaction
index 47d81151a..6e9e34a8c 100755
--- a/Projects/Strasse/reaction/C12_p2p.reaction
+++ b/Projects/Strasse/reaction/C12_p2p.reaction
@@ -1,4 +1,3 @@
-%%%%%%%%%%%%%%%%%% QFS Example in progress %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Beam
Particle= 12C
--
GitLab