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

* progress on strasse analysis

parent 47d712fd
No related branches found
No related tags found
No related merge requests found
...@@ -59,7 +59,6 @@ void TReactionConditions::Clear(){ ...@@ -59,7 +59,6 @@ void TReactionConditions::Clear(){
fRC_Momentum_Direction_X.clear(); fRC_Momentum_Direction_X.clear();
fRC_Momentum_Direction_Y.clear(); fRC_Momentum_Direction_Y.clear();
fRC_Momentum_Direction_Z.clear(); fRC_Momentum_Direction_Z.clear();
fRC_Momentum.clear();
} }
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
void TReactionConditions::Dump() const{ void TReactionConditions::Dump() const{
......
...@@ -53,6 +53,7 @@ private: ...@@ -53,6 +53,7 @@ private:
double fRC_Beam_Emittance_Theta; double fRC_Beam_Emittance_Theta;
double fRC_Beam_Emittance_Phi; double fRC_Beam_Emittance_Phi;
double fRC_Beam_Reaction_Energy; double fRC_Beam_Reaction_Energy;
double fRC_Vertex_Position_X; double fRC_Vertex_Position_X;
double fRC_Vertex_Position_Y; double fRC_Vertex_Position_Y;
double fRC_Vertex_Position_Z; double fRC_Vertex_Position_Z;
...@@ -68,8 +69,6 @@ private: ...@@ -68,8 +69,6 @@ private:
vector<double> fRC_Momentum_Direction_X; vector<double> fRC_Momentum_Direction_X;
vector<double> fRC_Momentum_Direction_Y; vector<double> fRC_Momentum_Direction_Y;
vector<double> fRC_Momentum_Direction_Z; vector<double> fRC_Momentum_Direction_Z;
vector<TVector3> fRC_Momentum;
public: public:
TReactionConditions(); TReactionConditions();
...@@ -82,7 +81,7 @@ public: ...@@ -82,7 +81,7 @@ public:
///////////////////// SETTERS //////////////////////// ///////////////////// SETTERS ////////////////////////
// Beam parameter // Beam parameter
void SetBeamParticleName (const string & Beam_Particle_Name) {fRC_Beam_Particle_Name = Beam_Particle_Name;}//! 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 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 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;}//! void SetBeamEmittanceThetaX (const double & Beam_Emittance_ThetaX) {fRC_Beam_Emittance_ThetaX = Beam_Emittance_ThetaX;}//!
...@@ -124,6 +123,8 @@ public: ...@@ -124,6 +123,8 @@ public:
double GetVertexPositionX () const {return fRC_Vertex_Position_X ;}//! double GetVertexPositionX () const {return fRC_Vertex_Position_X ;}//!
double GetVertexPositionY () const {return fRC_Vertex_Position_Y ;}//! double GetVertexPositionY () const {return fRC_Vertex_Position_Y ;}//!
double GetVertexPositionZ () const {return fRC_Vertex_Position_Z ;}//! 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 GetExcitation3 () const {return fRC_ExcitationEnergy3 ;}//!
double GetExcitation4 () const {return fRC_ExcitationEnergy4 ;}//! double GetExcitation4 () const {return fRC_ExcitationEnergy4 ;}//!
double GetThetaCM () const {return fRC_ThetaCM;}//! double GetThetaCM () const {return fRC_ThetaCM;}//!
...@@ -138,7 +139,8 @@ public: ...@@ -138,7 +139,8 @@ public:
double GetMomentumDirectionX (const int &i) const {return fRC_Momentum_Direction_X[i];}//! 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 GetMomentumDirectionY (const int &i) const {return fRC_Momentum_Direction_Y[i];}//!
double GetMomentumDirectionZ (const int &i) const {return fRC_Momentum_Direction_Z[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 GetBeamDirection () const ;
TVector3 GetParticleDirection (const int i) const ; TVector3 GetParticleDirection (const int i) const ;
......
...@@ -75,10 +75,12 @@ void Analysis::TreatEvent(){ ...@@ -75,10 +75,12 @@ void Analysis::TreatEvent(){
TVector3 InnerPos1 = Strasse->GetInnerPositionOfInteraction(0); TVector3 InnerPos1 = Strasse->GetInnerPositionOfInteraction(0);
TVector3 OuterPos1 = Strasse->GetOuterPositionOfInteraction(0); TVector3 OuterPos1 = Strasse->GetOuterPositionOfInteraction(0);
TVector3 Proton1 = OuterPos1-InnerPos1; TVector3 Proton1 = OuterPos1-InnerPos1;
Proton1=Proton1.Unit();
// Proton 2 // Proton 2
TVector3 InnerPos2 = Strasse->GetInnerPositionOfInteraction(1); TVector3 InnerPos2 = Strasse->GetInnerPositionOfInteraction(1);
TVector3 OuterPos2 = Strasse->GetOuterPositionOfInteraction(1); TVector3 OuterPos2 = Strasse->GetOuterPositionOfInteraction(1);
TVector3 Proton2 = OuterPos2-InnerPos2; TVector3 Proton2 = OuterPos2-InnerPos2;
Proton2=Proton2.Unit();
deltaPhi = abs(Proton1.Phi()/deg-Proton2.Phi()/deg); deltaPhi = abs(Proton1.Phi()/deg-Proton2.Phi()/deg);
sumTheta = Proton1.Theta()/deg+Proton2.Theta()/deg; sumTheta = Proton1.Theta()/deg+Proton2.Theta()/deg;
...@@ -105,34 +107,75 @@ void Analysis::TreatEvent(){ ...@@ -105,34 +107,75 @@ void Analysis::TreatEvent(){
unsigned int i1,i2; unsigned int i1,i2;
i1 = Catana->FindClosestHitToLine(InnerPos1,OuterPos1,d1); i1 = Catana->FindClosestHitToLine(InnerPos1,OuterPos1,d1);
i2 = Catana->FindClosestHitToLine(InnerPos2,OuterPos2,d2); 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]); E1 = ReconstructProtonEnergy(Vertex,Proton1,Catana->Energy[i1]);
E2 = ReconstructProtonEnergy(Vertex,Proton2,Catana->Energy[i2]); 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 // Vertex = RC->GetVertexPosition();
TVector3 PA(0,0,sqrt(TA*(TA+2*m_QFS->GetNucleusA()->Mass()))); // for like there is no BDC //TA = RC->GetBeamEnergy();
Proton1=E1*Proton1.Unit(); ////////////////////////////////////
Proton2=E2*Proton2.Unit();
// 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()); LV_A.SetVectM(PA,m_QFS->GetNucleusA()->Mass());
double P1= sqrt(E1*(E1+2*NPUNITS::proton_mass_c2)); double P1= sqrt(E1*(E1+2*NPUNITS::proton_mass_c2));
double P2= sqrt(E2*(E2+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_p1.SetVectM(Proton1.Unit()*P1,NPUNITS::proton_mass_c2);
LV_p2.SetVectM(Proton2.Unit()*P2,NPUNITS::proton_mass_c2); LV_p2.SetVectM(Proton2.Unit()*P2,NPUNITS::proton_mass_c2);
// computing Ex from Missing Mass // computing Ex from Missing Mass
LV_B = LV_A + LV_T - LV_p1 - LV_p2; LV_B = LV_A + LV_T - LV_p1 - LV_p2;
//LV_B = RC->GetParticleMomentum(2); //LV_B = RC->GetParticleMomentum(2);
Ex = LV_B.M() - m_QFS->GetNucleusB()->Mass(); Ex = LV_B.M() - m_QFS->GetNucleusB()->Mass();
} }
} }
} }
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
double Analysis::ReconstructProtonEnergy(const TVector3& x0, const TVector3& dir,const double& Ecatana){ 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()); Normal.SetPhi(dir.Phi());
double Theta = dir.Angle(Normal); double Theta = dir.Angle(Normal);
// Catana Al housing // Catana Al housing
...@@ -146,12 +189,11 @@ double Analysis::ReconstructProtonEnergy(const TVector3& x0, const TVector3& dir ...@@ -146,12 +189,11 @@ double Analysis::ReconstructProtonEnergy(const TVector3& x0, const TVector3& dir
// LH2 target // LH2 target
static TVector3 x1; static TVector3 x1;
x1= x0+dir; x1= x0+dir;
TVector3 T1(0,30,0); TVector3 T1(0,15,0);
TVector3 T2(0,30,1); TVector3 T2(0,15,1);
T1.SetPhi(dir.Phi()); T1.SetPhi(dir.Phi());
T2.SetPhi(dir.Phi()); T2.SetPhi(dir.Phi());
TVector3 Vertex,delta; double d = NPL::MinimumDistancePointLine(T1,T2,x0);
double d = NPL::MinimumDistanceTwoLines(x0,x1,T1,T2,Vertex,delta);
E = protonTarget.EvaluateInitialEnergy(E,d,Theta); E = protonTarget.EvaluateInitialEnergy(E,d,Theta);
} }
...@@ -185,6 +227,17 @@ void Analysis::InitOutputBranch() { ...@@ -185,6 +227,17 @@ void Analysis::InitOutputBranch() {
RootOutput::getInstance()->GetTree()->Branch("deltaPhi",&deltaPhi,"deltaPhi/D"); RootOutput::getInstance()->GetTree()->Branch("deltaPhi",&deltaPhi,"deltaPhi/D");
RootOutput::getInstance()->GetTree()->Branch("sumTheta",&sumTheta,"sumTheta/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("Distance",&Distance,"Distance/D");
RootOutput::getInstance()->GetTree()->Branch("InteractionCoordinates","TInteractionCoordinates",&DC); RootOutput::getInstance()->GetTree()->Branch("InteractionCoordinates","TInteractionCoordinates",&DC);
...@@ -212,6 +265,15 @@ void Analysis::ReInitValue(){ ...@@ -212,6 +265,15 @@ void Analysis::ReInitValue(){
Distance=-1000; Distance=-1000;
sumTheta=-1000; sumTheta=-1000;
deltaPhi=-1000; deltaPhi=-1000;
alpha=-1000;
Theta1=-1000;
Phi1=-1000;
Theta2=-1000;
Phi2=-1000;
Theta1s=-1000;
Phi1s=-1000;
Theta2s=-1000;
Phi2s=-1000;
} }
......
...@@ -50,7 +50,7 @@ class Analysis: public NPL::VAnalysis{ ...@@ -50,7 +50,7 @@ class Analysis: public NPL::VAnalysis{
void ReInitValue(); void ReInitValue();
static NPL::VAnalysis* Construct(); static NPL::VAnalysis* Construct();
TVector3 InterpolateInPlaneZ(TVector3,TVector3,double); TVector3 InterpolateInPlaneZ(TVector3,TVector3,double);
private: private:
double Ex; double Ex;
double E1; double E1;
...@@ -66,6 +66,17 @@ class Analysis: public NPL::VAnalysis{ ...@@ -66,6 +66,17 @@ class Analysis: public NPL::VAnalysis{
double Distance; double Distance;
double deltaPhi; double deltaPhi;
double sumTheta; 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_A;
TLorentzVector LV_T; TLorentzVector LV_T;
TLorentzVector LV_B; TLorentzVector LV_B;
......
...@@ -3,7 +3,7 @@ Target ...@@ -3,7 +3,7 @@ Target
THICKNESS= 150 mm THICKNESS= 150 mm
ANGLE= 0 deg ANGLE= 0 deg
RADIUS= 15 mm RADIUS= 15 mm
MATERIAL= LH2 MATERIAL= LH2
X= 0 mm X= 0 mm
Y= 0 mm Y= 0 mm
Z= 0 mm Z= 0 mm
......
void Shift(){ void Shift(){
TFile* file_ok = TFile::Open("../../Outputs/Analysis/strasse_ok.root"); TFile* file_ok = TFile::Open("../../Outputs/Analysis/strasse_ok.root");
TFile* file_shifted = TFile::Open("../../Outputs/Analysis/strasse_shifted.root"); TFile* file_shifted = TFile::Open("../../Outputs/Analysis/strasse_shifted.root");
TTree* ok= (TTree*) file_ok->FindObjectAny("PhysicsTree"); TTree* ok= (TTree*) file_ok->FindObjectAny("PhysicsTree");
TTree* shifted= (TTree*) file_shifted->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); TCanvas* ctheta= new TCanvas("ControlTheta","ControlTheta",2000,1000);
ctheta->Divide(2,1); ctheta->Divide(2,1);
ctheta->cd(1); ctheta->cd(1);
string cond = "Theta12!=-1000"; ok->Draw("Theta12>>ht(2000)",cond.c_str(),"") ;
ok->Draw("Theta12>>ht(5000)",cond.c_str(),"") ; shifted->Draw("Theta12>>hts(2000)",cond.c_str(),"same") ;
shifted->Draw("Theta12>>hts(5000)",cond.c_str(),"same") ;
TH1* hts= (TH1*) gDirectory->FindObjectAny("hts"); TH1* hts= (TH1*) gDirectory->FindObjectAny("hts");
hts->SetFillColor(kOrange-3); hts->SetFillColor(kOrange-3);
hts->SetLineColor(kOrange-3); hts->SetLineColor(kOrange-3);
ctheta->cd(2); ctheta->cd(2);
cond = "deltaPhi!=-1000"; cond = "deltaPhi!=-1000";
ok->Draw("deltaPhi>>hp(5000)",cond.c_str(),"") ; ok->Draw("deltaPhi>>hp(2000)",cond.c_str(),"") ;
shifted->Draw("deltaPhi>>hps(5000)",cond.c_str(),"same") ; shifted->Draw("deltaPhi>>hps(2000)",cond.c_str(),"same") ;
TH1* hps= (TH1*) gDirectory->FindObjectAny("hps"); TH1* hps= (TH1*) gDirectory->FindObjectAny("hps");
hps->SetFillColor(kOrange-3); hps->SetFillColor(kOrange-3);
hps->SetLineColor(kOrange-3); hps->SetLineColor(kOrange-3);
......
%%%%%%%%%%%%%%%%%% QFS Example in progress %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Beam Beam
Particle= 12C Particle= 12C
......
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