Skip to content
Snippets Groups Projects
Commit 0ec7fa26 authored by flavigny's avatar flavigny
Browse files

* Update STRASSE analysis and last geometry

parent f80501ce
No related branches found
No related tags found
No related merge requests found
Pipeline #234311 passed
...@@ -19,69 +19,69 @@ ...@@ -19,69 +19,69 @@
* * * *
* * * *
*****************************************************************************/ *****************************************************************************/
#include<iostream> #include <iostream>
using namespace std; using namespace std;
#include"Analysis.h" #include "Analysis.h"
#include"NPAnalysisFactory.h" #include "NPAnalysisFactory.h"
#include"NPDetectorManager.h" #include "NPDetectorManager.h"
#include"NPOptionManager.h" #include "NPFunction.h"
#include"NPFunction.h" #include "NPOptionManager.h"
#include"NPTrackingUtility.h" #include "NPPhysicalConstants.h"
#include"NPPhysicalConstants.h" #include "NPTrackingUtility.h"
#include"TRandom3.h" #include "TRandom3.h"
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
Analysis::Analysis(){ Analysis::Analysis() {}
}
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
Analysis::~Analysis(){ Analysis::~Analysis() {}
}
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
void Analysis::Init(){ void Analysis::Init() {
IC= new TInitialConditions; IC = new TInitialConditions;
DC= new TInteractionCoordinates; DC = new TInteractionCoordinates;
RC= new TReactionConditions; RC = new TReactionConditions;
InitOutputBranch(); InitOutputBranch();
InitInputBranch(); InitInputBranch();
Strasse = (TStrassePhysics*) m_DetectorManager -> GetDetector("Strasse"); Strasse = (TStrassePhysics*)m_DetectorManager->GetDetector("Strasse");
// Catana = (TCatanaPhysics*) m_DetectorManager -> GetDetector("Catana"); // Catana = (TCatanaPhysics*) m_DetectorManager -> GetDetector("Catana");
// reaction properties // reaction properties
m_QFS = new NPL::QFS(); m_QFS = new NPL::QFS();
m_QFS->ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile()); m_QFS->ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
// reaction properties // reaction properties
myBeam = new NPL::Beam(); myBeam = new NPL::Beam();
myBeam->ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile()); myBeam->ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
InitialBeamEnergy = myBeam->GetEnergy()/* * myBeam->GetA()*/; InitialBeamEnergy = myBeam->GetEnergy() /* * myBeam->GetA()*/;
// target thickness // target thickness
TargetThickness = m_DetectorManager->GetTargetThickness(); // TargetThickness = m_DetectorManager->GetTargetThickness();
string TargetMaterial = m_DetectorManager->GetTargetMaterial(); // string TargetMaterial = m_DetectorManager->GetTargetMaterial();
TargetThickness = 150 * mm;
string TargetMaterial = "LH2";
// EnergyLoss Tables // EnergyLoss Tables
string BeamName = NPL::ChangeNameToG4Standard(myBeam->GetName()); string BeamName = NPL::ChangeNameToG4Standard(myBeam->GetName());
BeamTarget = NPL::EnergyLoss(BeamName+"_"+TargetMaterial+".G4table","G4Table",100); BeamTarget = NPL::EnergyLoss(BeamName + "_" + TargetMaterial + ".G4table", "G4Table", 100);
protonTarget = NPL::EnergyLoss("proton_"+TargetMaterial+".G4table","G4Table",100); protonTarget = NPL::EnergyLoss("proton_" + TargetMaterial + ".G4table", "G4Table", 100);
protonAl = NPL::EnergyLoss("proton_Al.G4table","G4Table",100); protonAl = NPL::EnergyLoss("proton_Al.G4table", "G4Table", 100);
protonSi = NPL::EnergyLoss("proton_Si.G4table","G4Table",100); protonSi = NPL::EnergyLoss("proton_Si.G4table", "G4Table", 100);
LV_T.SetVectM(TVector3(0,0,0),NPUNITS::proton_mass_c2); LV_T.SetVectM(TVector3(0, 0, 0), NPUNITS::proton_mass_c2);
rand= new TRandom3(); rand = new TRandom3();
} }
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
void Analysis::TreatEvent(){ void Analysis::TreatEvent() {
// Reinitiate calculated variable // Reinitiate calculated variable
ReInitValue(); ReInitValue();
bool perfectBeamEnergy=false; bool perfectBeamEnergy = false;
bool perfectBeamTracking=false; bool perfectBeamTracking = false;
bool perfectProtonTracking=false; bool perfectProtonTracking = false;
bool perfectProtonEnergy=false; bool perfectProtonEnergy = false;
bool CATANAEnergyReco=false; bool CATANAEnergyReco = false;
unsigned int size = Strasse->GetEventMultiplicity(); unsigned int size = Strasse->GetEventMultiplicity();
if(size==2){ // 2 proton detected if (size == 2) { // 2 proton detected
//////////////////////////////////////////////// ////////////////////////////////////////////////
////////////Protons (E,P) Reconstruction /////// ////////////Protons (E,P) Reconstruction ///////
...@@ -89,86 +89,99 @@ void Analysis::TreatEvent(){ ...@@ -89,86 +89,99 @@ void Analysis::TreatEvent(){
// Proton 1 // Proton 1
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(); 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(); 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;
Theta12 = Proton1.Angle(Proton2)/deg; Theta12 = Proton1.Angle(Proton2) / deg;
// computing minimum distance of the two lines // computing minimum distance of the two lines
TVector3 Vertex; TVector3 Vertex;
TVector3 delta; TVector3 delta;
Distance = NPL::MinimumDistanceTwoLines(InnerPos1,OuterPos1,InnerPos2,OuterPos2,Vertex,delta); Distance = NPL::MinimumDistanceTwoLines(InnerPos1, OuterPos1, InnerPos2, OuterPos2, Vertex, delta);
VertexX=Vertex.X(); VertexX = Vertex.X();
VertexY=Vertex.Y(); VertexY = Vertex.Y();
VertexZ=Vertex.Z(); VertexZ = Vertex.Z();
deltaX=delta.X(); deltaX = delta.X();
deltaY=delta.Y(); deltaY = delta.Y();
deltaZ=delta.Z(); deltaZ = delta.Z();
if(perfectProtonTracking){ if (perfectProtonTracking) {
// From Reaction Conditions // From Reaction Conditions
VertexX=RC->GetVertexPositionX(); VertexX = RC->GetVertexPositionX();
VertexY=RC->GetVertexPositionY(); VertexY = RC->GetVertexPositionY();
VertexZ=RC->GetVertexPositionZ(); VertexZ = RC->GetVertexPositionZ();
Proton2=RC->GetParticleMomentum(0).Unit(); Proton2 = RC->GetParticleDirection(0).Unit();
Proton1=RC->GetParticleMomentum(1).Unit(); Proton1 = RC->GetParticleDirection(1).Unit();
deltaX=0; deltaX = 0;
deltaY=0; deltaY = 0;
deltaZ=0; deltaZ = 0;
} }
if(perfectProtonEnergy){ if (perfectProtonEnergy) {
E1 = RC->GetKineticEnergy(1); E1 = RC->GetKineticEnergy(1);
E2 = RC->GetKineticEnergy(0); E2 = RC->GetKineticEnergy(0);
} }
else { else {
//General CATANA resolution for proton (1.4 % FWHM) // General CATANA resolution for proton (1.4 % FWHM)
// simply applied to E, crude approximation (no CATANA reco.) // simply applied to E, crude approximation (no CATANA reco.)
E1 = ApplyCATANAResoProton(RC->GetKineticEnergy(1)); E1 = ApplyCATANAResoProton(RC->GetKineticEnergy(1));
E2 = ApplyCATANAResoProton(RC->GetKineticEnergy(0)); E2 = ApplyCATANAResoProton(RC->GetKineticEnergy(0));
} }
if(CATANAEnergyReco){ if (CATANAEnergyReco) {
//if true, bypass previous proton energy calc. and use full CATANA Reco. // if true, bypass previous proton energy calc. and use full CATANA Reco.
/////////// ///////////
// TO DO // // TO DO //
/////////// ///////////
} }
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));
/////////////////////////////////////// ///////////////////////////////////////
//////Beam (E,P) Reconstruction /////// //////Beam (E,P) Reconstruction ///////
/////////////////////////////////////// ///////////////////////////////////////
TVector3 BeamDirection; TVector3 BeamDirection;
if(perfectBeamTracking) BeamDirection = RC->GetBeamDirection(); if (perfectBeamTracking)
else BeamDirection = TVector3(0,0,1); BeamDirection = RC->GetBeamDirection();
//Beam Energy at vertex else
BeamDirection = TVector3(0, 0, 1);
// Beam Energy at vertex
TA = RC->GetBeamEnergy(); TA = RC->GetBeamEnergy();
if(perfectBeamEnergy) TAcalc = TA; if (perfectBeamEnergy) {
else TAcalc = BeamTarget.Slow(InitialBeamEnergy,abs(VertexZ+75), BeamDirection.Angle(TVector3(0,0,1))); TAcalc = TA;
double beam_mom=sqrt(TAcalc*(TAcalc+2*m_QFS->GetParticleA()->Mass())); }
TVector3 PA=beam_mom*BeamDirection.Unit(); else if (VertexZ > -100) {
TAcalc = BeamTarget.Slow(InitialBeamEnergy, abs(VertexZ + 75), BeamDirection.Angle(TVector3(0, 0, 1)));
}
double beam_mom = sqrt(TAcalc * (TAcalc + 2 * m_QFS->GetParticleA()->Mass()));
TVector3 PA = beam_mom * BeamDirection.Unit();
///////////////////////////////////
///// Angles Protons / Beam //////
///////////////////////////////////
Theta1 = Proton1.Angle(BeamDirection) / deg;
Theta2 = Proton2.Angle(BeamDirection) / deg;
////////////////////////////////////////////// //////////////////////////////////////////////
///// Missing mass using Lorentz Vector ////// ///// Missing mass using Lorentz Vector //////
////////////////////////////////////////////// //////////////////////////////////////////////
LV_A.SetVectM(PA,m_QFS->GetParticleA()->Mass()); LV_A.SetVectM(PA, m_QFS->GetParticleA()->Mass());
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->GetParticleB()->Mass(); Ex = LV_B.M() - m_QFS->GetParticleB()->Mass();
}//endif size==2 } // endif size==2
} }
/* /////////////////////////////////// /* ///////////////////////////////////
...@@ -187,15 +200,15 @@ cout<<"Catana mult="<<Catana->GetEventMultiplicity()<<endl; ...@@ -187,15 +200,15 @@ cout<<"Catana mult="<<Catana->GetEventMultiplicity()<<endl;
//if(i1!=i2){ //if(i1!=i2){
if(true){ if(true){
double TA = BeamTarget.Slow(InitialBeamEnergy,abs(VertexZ-75),RC->GetBeamDirection().Angle(TVector3(0,0,1))); double TA = BeamTarget.Slow(InitialBeamEnergy,abs(VertexZ-75),RC->GetBeamDirection().Angle(TVector3(0,0,1)));
//////////////////////////////////// ////////////////////////////////////
// From Reaction Conditions // From Reaction Conditions
double E1s = RC->GetKineticEnergy(0); double E1s = RC->GetKineticEnergy(0);
double E2s = RC->GetKineticEnergy(1); double E2s = RC->GetKineticEnergy(1);
TVector3 Proton1s=RC->GetParticleMomentum(0).Unit(); TVector3 Proton1s=RC->GetParticleMomentum(0).Unit();
TVector3 Proton2s=RC->GetParticleMomentum(1).Unit(); TVector3 Proton2s=RC->GetParticleMomentum(1).Unit();
// Matching the right energy with the right proton // Matching the right energy with the right proton
if((Proton1s-Proton1).Mag()<(Proton1s-Proton2).Mag()){ if((Proton1s-Proton1).Mag()<(Proton1s-Proton2).Mag()){
E1=E1s; E1=E1s;
E2=E2s; E2=E2s;
...@@ -208,7 +221,7 @@ cout<<"Catana mult="<<Catana->GetEventMultiplicity()<<endl; ...@@ -208,7 +221,7 @@ cout<<"Catana mult="<<Catana->GetEventMultiplicity()<<endl;
Phi1s=Proton1s.Phi(); Phi1s=Proton1s.Phi();
Theta2s=Proton2s.Theta(); Theta2s=Proton2s.Theta();
Phi2s=Proton2s.Phi(); Phi2s=Proton2s.Phi();
} }
else{ else{
E2=E1s; E2=E1s;
...@@ -225,169 +238,161 @@ cout<<"Catana mult="<<Catana->GetEventMultiplicity()<<endl; ...@@ -225,169 +238,161 @@ cout<<"Catana mult="<<Catana->GetEventMultiplicity()<<endl;
} }
// From detectors // 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]);
//E1 = Catana->Energy[i1]; //E1 = Catana->Energy[i1];
//E2 = Catana->Energy[i2]; //E2 = Catana->Energy[i2];
cout<<"------------------------"<<endl; cout<<"------------------------"<<endl;
*/ */
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
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,1,0); 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
double E = protonAl.EvaluateInitialEnergy(Ecatana,0.5*mm,Theta); double E = protonAl.EvaluateInitialEnergy(Ecatana, 0.5 * mm, Theta);
cout<<"Ecatana="<<Ecatana<<endl; cout << "Ecatana=" << Ecatana << endl;
cout<<"Erec0="<<E<<endl; cout << "Erec0=" << E << endl;
// Strasse Chamber // Strasse Chamber
//E = protonAl.EvaluateInitialEnergy(E,3*mm,Theta); // E = protonAl.EvaluateInitialEnergy(E,3*mm,Theta);
// Outer Barrel // Outer Barrel
E = protonSi.EvaluateInitialEnergy(E,300*micrometer,Theta); E = protonSi.EvaluateInitialEnergy(E, 300 * micrometer, Theta);
// Inner Barrel // Inner Barrel
E = protonSi.EvaluateInitialEnergy(E,200*micrometer,Theta); E = protonSi.EvaluateInitialEnergy(E, 200 * micrometer, Theta);
// LH2 target // LH2 target
static TVector3 x1; static TVector3 x1;
x1= x0+dir; x1 = x0 + dir;
TVector3 T1(0,15,0); TVector3 T1(0, 15, 0);
TVector3 T2(0,15,1); TVector3 T2(0, 15, 1);
T1.SetPhi(dir.Phi()); T1.SetPhi(dir.Phi());
T2.SetPhi(dir.Phi()); T2.SetPhi(dir.Phi());
// double d = NPL::MinimumDistancePointLine(T1,T2,x0); // double d = NPL::MinimumDistancePointLine(T1,T2,x0);
double d = 0; double d = 0;
cout<<"d="<<d<<endl; cout << "d=" << d << endl;
cout<<"Theta="<<Theta<<endl; cout << "Theta=" << Theta << endl;
E = protonTarget.EvaluateInitialEnergy(E,d,Theta); E = protonTarget.EvaluateInitialEnergy(E, d, Theta);
cout<<"Erec="<<E<<endl; cout << "Erec=" << E << endl;
return E; return E;
//return 0; // return 0;
} }
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
double Analysis::ApplyCATANAResoProton(double Ein){ double Analysis::ApplyCATANAResoProton(double Ein) {
// Function applying overall CATANA crystal resolution // Function applying overall CATANA crystal resolution
// For protons (1.4% FWHM) // For protons (1.4% FWHM)
double FWHM = 1.4/100 * Ein; double FWHM = 1.4 / 100 * Ein;
double sigma = FWHM/2.35; double sigma = FWHM / 2.35;
double Eout = rand->Gaus(Ein,sigma); double Eout = rand->Gaus(Ein, sigma);
return Eout; return Eout;
} }
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
double Analysis::ApplyCATANAResoGamma(double Ein){ double Analysis::ApplyCATANAResoGamma(double Ein) {
// Function applying overall CATANA crystal resolution // Function applying overall CATANA crystal resolution
// For gammas defined in smsimulator package // For gammas defined in smsimulator package
double a = 0.686569; double a = 0.686569;
double b = 0.564352; double b = 0.564352;
// Ein from MeV to keV // Ein from MeV to keV
Ein = Ein * 1000; Ein = Ein * 1000;
double SigmaE = a * pow(Ein,b); double SigmaE = a * pow(Ein, b);
double Eout = rand->Gaus(Ein,SigmaE); double Eout = rand->Gaus(Ein, SigmaE);
return Eout/1000.; return Eout / 1000.;
} }
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
TVector3 Analysis::InterpolateInPlaneZ(TVector3 V0, TVector3 V1, double Zproj){ TVector3 Analysis::InterpolateInPlaneZ(TVector3 V0, TVector3 V1, double Zproj) {
TVector3 Vproj(-999,-999,-999); TVector3 Vproj(-999, -999, -999);
double t = (Zproj - V1.Z()) / (V1.Z()-V0.Z()); double t = (Zproj - V1.Z()) / (V1.Z() - V0.Z());
double Xproj= V1.X() + (V1.X()-V0.X()) * t; double Xproj = V1.X() + (V1.X() - V0.X()) * t;
double Yproj= V1.Y() + (V1.Y()-V0.Y()) * t; double Yproj = V1.Y() + (V1.Y() - V0.Y()) * t;
Vproj.SetXYZ(Xproj,Yproj,Zproj); Vproj.SetXYZ(Xproj, Yproj, Zproj);
return Vproj; return Vproj;
} }
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
void Analysis::End(){ void Analysis::End() {}
}
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
void Analysis::InitOutputBranch() { void Analysis::InitOutputBranch() {
RootOutput::getInstance()->GetTree()->Branch("Ex",&Ex,"Ex/D"); RootOutput::getInstance()->GetTree()->Branch("Ex", &Ex, "Ex/D");
RootOutput::getInstance()->GetTree()->Branch("E1",&E1,"E1/D"); RootOutput::getInstance()->GetTree()->Branch("E1", &E1, "E1/D");
RootOutput::getInstance()->GetTree()->Branch("E2",&E2,"E2/D"); RootOutput::getInstance()->GetTree()->Branch("E2", &E2, "E2/D");
RootOutput::getInstance()->GetTree()->Branch("Theta12",&Theta12,"Theta12/D"); RootOutput::getInstance()->GetTree()->Branch("Theta12", &Theta12, "Theta12/D");
RootOutput::getInstance()->GetTree()->Branch("ThetaCM",&ThetaCM,"ThetaCM/D"); RootOutput::getInstance()->GetTree()->Branch("ThetaCM", &ThetaCM, "ThetaCM/D");
RootOutput::getInstance()->GetTree()->Branch("VertexX",&VertexX,"VertexX/D"); RootOutput::getInstance()->GetTree()->Branch("VertexX", &VertexX, "VertexX/D");
RootOutput::getInstance()->GetTree()->Branch("VertexY",&VertexY,"VertexY/D"); RootOutput::getInstance()->GetTree()->Branch("VertexY", &VertexY, "VertexY/D");
RootOutput::getInstance()->GetTree()->Branch("VertexZ",&VertexZ,"VertexZ/D"); RootOutput::getInstance()->GetTree()->Branch("VertexZ", &VertexZ, "VertexZ/D");
RootOutput::getInstance()->GetTree()->Branch("deltaX",&deltaX,"deltaX/D"); RootOutput::getInstance()->GetTree()->Branch("deltaX", &deltaX, "deltaX/D");
RootOutput::getInstance()->GetTree()->Branch("deltaY",&deltaY,"deltaY/D"); RootOutput::getInstance()->GetTree()->Branch("deltaY", &deltaY, "deltaY/D");
RootOutput::getInstance()->GetTree()->Branch("deltaZ",&deltaZ,"deltaZ/D"); RootOutput::getInstance()->GetTree()->Branch("deltaZ", &deltaZ, "deltaZ/D");
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("alpha", &alpha, "alpha/D");
RootOutput::getInstance()->GetTree()->Branch("Theta1",&Theta1,"Theta1/D"); RootOutput::getInstance()->GetTree()->Branch("Theta1", &Theta1, "Theta1/D");
RootOutput::getInstance()->GetTree()->Branch("Phi1",&Phi1,"Phi1/D"); RootOutput::getInstance()->GetTree()->Branch("Phi1", &Phi1, "Phi1/D");
RootOutput::getInstance()->GetTree()->Branch("Theta2",&Theta2,"Theta2/D"); RootOutput::getInstance()->GetTree()->Branch("Theta2", &Theta2, "Theta2/D");
RootOutput::getInstance()->GetTree()->Branch("Phi2",&Phi2,"Phi2/D"); RootOutput::getInstance()->GetTree()->Branch("Phi2", &Phi2, "Phi2/D");
RootOutput::getInstance()->GetTree()->Branch("Theta1s",&Theta1s,"Theta1s/D"); RootOutput::getInstance()->GetTree()->Branch("Theta1s", &Theta1s, "Theta1s/D");
RootOutput::getInstance()->GetTree()->Branch("Phi1s",&Phi1s,"Phi1s/D"); RootOutput::getInstance()->GetTree()->Branch("Phi1s", &Phi1s, "Phi1s/D");
RootOutput::getInstance()->GetTree()->Branch("Theta2s",&Theta2s,"Theta2s/D"); RootOutput::getInstance()->GetTree()->Branch("Theta2s", &Theta2s, "Theta2s/D");
RootOutput::getInstance()->GetTree()->Branch("Phi2s",&Phi2s,"Phi2s/D"); RootOutput::getInstance()->GetTree()->Branch("Phi2s", &Phi2s, "Phi2s/D");
RootOutput::getInstance()->GetTree()->Branch("TA",&TA,"TA/D"); RootOutput::getInstance()->GetTree()->Branch("TA", &TA, "TA/D");
RootOutput::getInstance()->GetTree()->Branch("TAcalc",&TAcalc,"TAcalc/D"); RootOutput::getInstance()->GetTree()->Branch("TAcalc", &TAcalc, "TAcalc/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); RootOutput::getInstance()->GetTree()->Branch("ReactionConditions", "TReactionConditions", &RC);
RootOutput::getInstance()->GetTree()->Branch("ReactionConditions","TReactionConditions",&RC);
} }
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
void Analysis::InitInputBranch(){ void Analysis::InitInputBranch() {
RootInput:: getInstance()->GetChain()->SetBranchAddress("InteractionCoordinates",&DC); RootInput::getInstance()->GetChain()->SetBranchAddress("InteractionCoordinates", &DC);
RootInput:: getInstance()->GetChain()->SetBranchAddress("ReactionConditions",&RC); RootInput::getInstance()->GetChain()->SetBranchAddress("ReactionConditions", &RC);
} }
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
void Analysis::ReInitValue(){ void Analysis::ReInitValue() {
Ex = -1000 ; Ex = -1000;
E1= -1000; E1 = -1000;
E2 = -1000; E2 = -1000;
Theta12 = -1000; Theta12 = -1000;
ThetaCM = -1000; ThetaCM = -1000;
VertexX=-1000; VertexX = -1000;
VertexY=-1000; VertexY = -1000;
VertexZ=-1000; VertexZ = -1000;
deltaX=-1000; deltaX = -1000;
deltaY=-1000; deltaY = -1000;
deltaZ=-1000; deltaZ = -1000;
Distance=-1000; Distance = -1000;
sumTheta=-1000; sumTheta = -1000;
deltaPhi=-1000; deltaPhi = -1000;
alpha=-1000; alpha = -1000;
Theta1=-1000; Theta1 = -1000;
Phi1=-1000; Phi1 = -1000;
Theta2=-1000; Theta2 = -1000;
Phi2=-1000; Phi2 = -1000;
Theta1s=-1000; Theta1s = -1000;
Phi1s=-1000; Phi1s = -1000;
Theta2s=-1000; Theta2s = -1000;
Phi2s=-1000; Phi2s = -1000;
TA=-1000; TA = -1000;
TAcalc=-1000; TAcalc = -1000;
} }
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
// Construct Method to be pass to the AnalysisFactory // // Construct Method to be pass to the AnalysisFactory //
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
NPL::VAnalysis* Analysis::Construct(){ NPL::VAnalysis* Analysis::Construct() { return (NPL::VAnalysis*)new Analysis(); }
return (NPL::VAnalysis*) new Analysis();
}
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
// Registering the construct method to the factory // // Registering the construct method to the factory //
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
extern "C"{ extern "C" {
class proxy{ class proxy {
public: public:
proxy(){ proxy() { NPL::AnalysisFactory::getInstance()->SetConstructor(Analysis::Construct); }
NPL::AnalysisFactory::getInstance()->SetConstructor(Analysis::Construct);
}
}; };
proxy p; proxy p;
......
EmPhysicsList Option4 EmPhysicsList Option4
DefaultCutOff 1000 DefaultCutOff 1000000
IonBinaryCascadePhysics 0 IonBinaryCascadePhysics 0
NPIonInelasticPhysics 0 NPIonInelasticPhysics 0
EmExtraPhysics 0 EmExtraPhysics 0
...@@ -8,4 +8,4 @@ StoppingPhysics 0 ...@@ -8,4 +8,4 @@ StoppingPhysics 0
OpticalPhysics 0 OpticalPhysics 0
HadronPhysicsINCLXX 0 HadronPhysicsINCLXX 0
HadronPhysicsQGSP_BIC_HP 0 HadronPhysicsQGSP_BIC_HP 0
Decay 1 Decay 0
This diff is collapsed.
{ {
gStyle->SetPalette(1); gStyle->SetPalette(1);
//TFile *file= new TFile("../../Outputs/Analysis/PhysicsTree.root"); // TFile *file= new TFile("../../Outputs/Analysis/PhysicsTree.root");
TFile *file= new TFile("../../Outputs/Analysis/configJuly2021_ana.root"); // TFile *file= new TFile("../../Outputs/Analysis/configJuly2021_ana_1MeV_perfect_Decay.root");
TTree *tree = (TTree*)file->Get("PhysicsTree"); // TFile* file = new TFile("../../Outputs/Analysis/configJuly2021_ana_1MeV.root");
// TFile* file = new TFile("../../Outputs/Analysis/configJune2022_ana_0MeV_real_100MeVc_50um.root");
// TFile* file = new TFile("../../Outputs/Analysis/configJune2022_ana_0MeV_real_100MeVc_I200um_O300um.root");
// TFile* file = new TFile("../../Outputs/Analysis/configJune2022_ana_0MeV_real_100MeVc_I200um_O300um_newtar.root");
TFile* file = new TFile("../../Outputs/Analysis/test_newtar_ana.root");
// TFile* file = new TFile("../../Outputs/Analysis/test_oldtar_ana.root");
TCanvas *c1 = new TCanvas("c1","c1",1000,1000); TTree* tree = (TTree*)file->Get("PhysicsTree");
c1->Divide(3,3);
c1->cd(1);
tree->Draw("Ex>>h1(200,-10,10)","","");
c1->cd(2);
tree->Draw("fRC_Vertex_Position_X:fRC_Vertex_Position_Z>>h3(400,-150,250,200,-100,100)","","colz");
c1->cd(3);
tree->Draw("VertexX:VertexZ>>h4(400,-150,250,200,-100,100)","","colz");
c1->cd(4);
tree->Draw("fRC_Vertex_Position_X-VertexX>>h5(100,-2,2)","","");
c1->cd(5);
tree->Draw("fRC_Vertex_Position_Y-VertexY>>h6(100,-2,2)","","");
c1->cd(6);
tree->Draw("fRC_Vertex_Position_Z-VertexZ>>h7(100,-2,2)","","");
c1->cd(7);
//tree->Draw("fRC_Kinetic_Energy[1]-Catana.Energy[0]>>h8(100,-10,10)","Catana.GetEventMultiplicity()==2","");
tree->Draw("fRC_Kinetic_Energy[1]-E1>>h8(100,-50,50)","","");
c1->cd(8);
//tree->Draw("fRC_Kinetic_Energy[1]-Catana.Energy[1]>>h9(100,-10,10)","Catana.GetEventMultiplicity()==2","");
tree->Draw("fRC_Kinetic_Energy[0]-E2>>h9(100,-50,50)","","");
c1->cd(9);
tree->Draw("fRC_Kinetic_Energy[0]:E2>>h10(300,0,300,300,0,300)","","colz");
TCanvas* c1 = new TCanvas("c1", "c1", 1000, 1000);
c1->Divide(3, 4);
c1->cd(1);
tree->Draw("Ex>>h1(400,-10,10)", "", "");
h1->SetXTitle("Ex_{missingmass}");
h1->SetYTitle("counts / 50 keV");
c1->cd(2);
tree->Draw("fRC_Vertex_Position_X:fRC_Vertex_Position_Z>>h3(400,-150,250,200,-100,100)", "", "colz");
TBox* box = new TBox(-75, -15, 75, 15);
box->SetFillStyle(0);
box->SetLineColor(kRed);
box->Draw("same");
c1->cd(3);
tree->Draw("VertexX:VertexZ>>h4(400,-150,250,200,-100,100)", "", "colz");
box->Draw("same");
double xmin = -75;
int xmax = 75;
double ymin = -15;
int ymax = 15;
int binx1 = h4->GetXaxis()->FindBin(xmin);
int binx2 = h4->GetXaxis()->FindBin(xmax);
int biny1 = h4->GetYaxis()->FindBin(ymin);
int biny2 = h4->GetYaxis()->FindBin(ymax);
double integralall = h3->Integral(binx1, binx2, biny1, biny2);
double integralmult2 = h4->Integral();
double integralrec = h4->Integral(binx1, binx2, biny1, biny2);
double integralEx = h1->Integral();
cout << "Simulated = " << integralall << endl;
cout << "Mult2 from vertex = " << integralmult2 << " (" << integralmult2 / integralall * 100 << "%)" << endl;
cout << "Rec. within target = " << integralrec << endl;
cout << "Rec. within target(%) = " << integralrec / integralall * 100 << endl;
cout << "Rec. Ex (%) = " << integralEx / integralall * 100 << endl;
c1->cd(4);
tree->Draw("fRC_Vertex_Position_X-VertexX>>h5(100,-2,2)", "", "");
c1->cd(5);
tree->Draw("fRC_Vertex_Position_Y-VertexY>>h6(100,-2,2)", "", "");
c1->cd(6);
tree->Draw("fRC_Vertex_Position_Z-VertexZ>>h7(100,-2,2)", "", "");
c1->cd(7);
// tree->Draw("fRC_Kinetic_Energy[1]-Catana.Energy[0]>>h8(100,-10,10)","Catana.GetEventMultiplicity()==2","");
tree->Draw("fRC_Kinetic_Energy[1]-E1>>h8(100,-50,50)", "", "");
c1->cd(8);
// tree->Draw("fRC_Kinetic_Energy[1]-Catana.Energy[1]>>h9(100,-10,10)","Catana.GetEventMultiplicity()==2","");
tree->Draw("fRC_Kinetic_Energy[0]-E2>>h9(100,-50,50)", "", "");
c1->cd(9);
tree->Draw("fRC_Kinetic_Energy[0]:E2>>h10(300,0,300,300,0,300)", "", "colz");
c1->cd(10);
tree->Draw("Strasse.EventMultiplicity>>h11(10,0,10)", "", "colz");
double integralM2 = h11->Integral(3, 3);
double integralM2bis = h11->Integral(2, 3);
cout << "Mult2 Strasse.EventMult = " << integralM2 << " (" << integralM2 / integralall * 100 << "%)" << endl;
cout << "Mult1and2 Strasse.EventMult = " << integralM2bis << " (" << integralM2bis / integralall * 100 << "%)"
<< endl;
} }
{ {
gStyle->SetPalette(1); gStyle->SetPalette(1);
//TFile *file= new TFile("../../Outputs/Simulation/SimulatedTree.root"); // TFile *file= new TFile("../../Outputs/Simulation/SimulatedTree.root");
TFile *file= new TFile("../../Outputs/Simulation/configJuly2021.root"); // TFile *file= new TFile("../../Outputs/Simulation/configJuly2021_1MeV.root");
TTree *tree = (TTree*)file->Get("SimulatedTree"); // TFile* file = new TFile("../../Outputs/Simulation/configJune2022_0MeV_real_100MeVc_50um.root");
/// TFile* file = new TFile("../../Outputs/Simulation/configJune2022_0MeV_real_100MeVc_I200um_O300um.root");
// TFile* file = new TFile("../../Outputs/Simulation/configJune2022_0MeV_real_100MeVc_Ip50um_Op50um.root");
// TFile* file = new TFile("../../Outputs/Simulation/configJune2022_partial_0MeV_pencil_100MeVc_I200um_O300um.root");
// TFile* file = new TFile("../../Outputs/Simulation/configJune2022_0MeV_pencil.root");
// TFile* file = new TFile("../../Outputs/Simulation/configJune2022_50MeV_pencil.root");
// TFile* file = new TFile("../../Outputs/Simulation/configJune2022_100MeV_pencil.root");
//
// TFile* file = new TFile("../../Outputs/Simulation/test_newtar.root");
TFile* file = new TFile("../../Outputs/Simulation/test_oldtar.root");
TTree* tree = (TTree*)file->Get("SimulatedTree");
TCanvas *c1 = new TCanvas("c1","c1",1000,1000); TCanvas* c1 = new TCanvas("c1", "c1", 1000, 1000);
c1->Divide(4,4); c1->Divide(4, 4);
c1->cd(1); c1->cd(1);
tree->Draw("fRC_Beam_Reaction_Energy:fRC_Vertex_Position_Z>>h1(200,-100,100,650,3200,4500)","","colz"); tree->Draw("fRC_Beam_Reaction_Energy:fRC_Vertex_Position_Z>>h1(200,-100,100,650,3200,4500)", "", "colz");
c1->cd(2); c1->cd(2);
tree->Draw("fRC_Vertex_Position_Y:fRC_Vertex_Position_Z>>h2(200,-100,100,120,-60,60)","","colz"); tree->Draw("fRC_Vertex_Position_Y:fRC_Vertex_Position_Z>>h2(200,-100,100,120,-60,60)", "", "colz");
c1->cd(3); c1->cd(3);
tree->Draw("fRC_Vertex_Position_X:fRC_Vertex_Position_Z>>h3(200,-100,100,120,-60,60)","","colz"); tree->Draw("fRC_Vertex_Position_X:fRC_Vertex_Position_Z>>h3(200,-100,100,120,-60,60)", "", "colz");
c1->cd(4); c1->cd(4);
tree->Draw("fDetected_Position_Y:fDetected_Position_X>>h4","","colz"); tree->Draw("fDetected_Position_Y:fDetected_Position_X>>h4", "", "colz");
c1->cd(5); c1->cd(5);
tree->Draw("fDetected_Position_Y:fDetected_Position_Z>>h5(900,-70,230,700,-70,70)","","colz"); tree->Draw("fDetected_Position_Y:fDetected_Position_Z>>h5(900,-70,230,700,-70,70)", "", "colz");
c1->cd(6); c1->cd(6);
tree->Draw("fInner_TE_StripNbr:fDetected_Position_Z>>h6","","colz"); tree->Draw("fInner_TE_StripNbr:fDetected_Position_Z>>h6", "", "colz");
c1->cd(7); c1->cd(7);
tree->Draw("fInner_LE_StripNbr:fDetected_Position_X>>h7","","colz"); tree->Draw("fInner_LE_StripNbr:fDetected_Position_X>>h7", "", "colz");
c1->cd(8); c1->cd(8);
tree->Draw("Strasse.GetOuterMultTEnergy()+Strasse.GetInnerMultTEnergy()>>h8(6,0,6)","",""); tree->Draw("Strasse.GetOuterMultTEnergy()+Strasse.GetInnerMultTEnergy()>>h8(6,0,6)", "", "");
cout<<"MULT4="<<h8->GetBinContent(5)<<endl; cout << "MULT4=" << h8->GetBinContent(5) << endl;
cout<<"Efficiency2p="<<h8->GetBinContent(5)/50000*100<<endl; cout << "Efficiency2p=" << h8->GetBinContent(5) / h1->Integral() * 100 << endl;
c1->cd(9); c1->cd(9);
tree->Draw("fInner_LE_StripNbr>>h9(1250,0,1250)","",""); tree->Draw("fInner_LE_StripNbr>>h9(1250,0,1250)", "", "");
c1->cd(10); c1->cd(10);
tree->Draw("fInner_TE_StripNbr>>h10(1250,0,1250)","",""); tree->Draw("fInner_TE_StripNbr>>h10(1250,0,1250)", "", "");
c1->cd(11); c1->cd(11);
tree->Draw("fOuter_LE_StripNbr>>h11(1250,0,1250)","",""); tree->Draw("fOuter_LE_StripNbr>>h11(1250,0,1250)", "", "");
c1->cd(12); c1->cd(12);
tree->Draw("fOuter_TE_StripNbr>>h12(1250,0,1250)","",""); tree->Draw("fOuter_TE_StripNbr>>h12(1250,0,1250)", "", "");
c1->cd(13); c1->cd(13);
tree->Draw("fOuter_TE_StripNbr:fDetected_Position_Z>>h13","","colz"); tree->Draw("fOuter_TE_StripNbr:fDetected_Position_Z>>h13", "", "colz");
c1->cd(14); c1->cd(14);
tree->Draw("fOuter_LE_StripNbr:fDetected_Position_X>>h14","","colz"); tree->Draw("fOuter_LE_StripNbr:fDetected_Position_X>>h14", "", "colz");
c1->cd(15); c1->cd(15);
tree->Draw("fRC_Theta[1]>>h15","",""); tree->Draw("fRC_Theta[1]>>h15", "", "");
tree->Draw("fRC_Theta[1]>>h16","Strasse.GetOuterMultTEnergy()+Strasse.GetInnerMultTEnergy()","same"); tree->Draw("fRC_Theta[1]>>h16", "Strasse.GetOuterMultTEnergy()+Strasse.GetInnerMultTEnergy()", "same");
h16->SetLineColor(2); h16->SetLineColor(2);
h16->Scale(1./4);; h16->Scale(1. / 4);
c1->cd(16); ;
tree->Draw("fRC_Theta[0]>>h17","",""); c1->cd(16);
tree->Draw("fRC_Theta[0]>>h18","Strasse.GetOuterMultTEnergy()+Strasse.GetInnerMultTEnergy()","same"); tree->Draw("fRC_Theta[0]>>h17", "", "");
h18->SetLineColor(2); tree->Draw("fRC_Theta[0]>>h18", "Strasse.GetOuterMultTEnergy()+Strasse.GetInnerMultTEnergy()", "same");
h18->Scale(1./4);; h18->SetLineColor(2);
h18->Scale(1. / 4);
;
/* TCanvas* c2 = new TCanvas("c2", "c2", 1000, 1000);
TCanvas *c2 = new TCanvas("c2","c2",1000,1000); c2->Divide(2, 2);
c2->Divide(2,2); c2->cd(1);
c2->cd(1); tree->Draw("Strasse.GetInnerMultTEnergy()>>hMultInnerT(4,0,4)", "", "");
tree->Draw("Strasse.GetInnerMultTEnergy()>>hMultInnerT(4,0,4)","",""); c2->cd(2);
c2->cd(2); tree->Draw("Strasse.GetInnerMultLEnergy()>>hMultInnerL(4,0,4)", "", "");
tree->Draw("Strasse.GetInnerMultLEnergy()>>hMultInnerL(4,0,4)","",""); c2->cd(3);
c2->cd(3); tree->Draw("Strasse.GetOuterMultTEnergy()>>hMultOuterT(4,0,4)", "", "");
tree->Draw("Strasse.GetOuterMultTEnergy()>>hMultOuterT(4,0,4)","",""); c2->cd(4);
c2->cd(4); tree->Draw("Strasse.GetOuterMultLEnergy()>>hMultOuterL(4,0,4)", "", "");
tree->Draw("Strasse.GetOuterMultLEnergy()>>hMultOuterL(4,0,4)","","");
TCanvas* c3 = new TCanvas("c3", "c3", 1000, 1000);
c3->Divide(3, 2);
TCanvas *c3 = new TCanvas("c3","c3",1000,1000); c3->cd(1);
c3->Divide(3,2); tree->Draw("fRC_Kinetic_Energy[0]:fRC_Theta[0]>>hETheta(900,0,90,1000,0,350)", "", "colz");
c3->cd(2);
c3->cd(1); tree->Draw("fRC_Kinetic_Energy[0]:fRC_Theta[0]>>hETheta_MultInnerT(900,0,90,1000,0,350)",
tree->Draw("fRC_Kinetic_Energy[0]:fRC_Theta[0]>>hETheta(900,0,90,1000,0,350)","","colz"); "Strasse.GetInnerMultTEnergy()==1", "colz");
c3->cd(2); c3->cd(3);
tree->Draw("fRC_Kinetic_Energy[0]:fRC_Theta[0]>>hETheta_MultInnerT(900,0,90,1000,0,350)","Strasse.GetInnerMultTEnergy()==1","colz"); tree->Draw("fRC_Kinetic_Energy[0]:fRC_Theta[0]>>hETheta_MultInnerL(900,0,90,1000,0,350)",
c3->cd(3); "Strasse.GetInnerMultLEnergy()==1", "colz");
tree->Draw("fRC_Kinetic_Energy[0]:fRC_Theta[0]>>hETheta_MultInnerL(900,0,90,1000,0,350)","Strasse.GetInnerMultLEnergy()==1","colz"); c3->cd(4);
c3->cd(4); tree->Draw("fRC_Kinetic_Energy[0]:fRC_Theta[0]>>hETheta_MultOuterT(900,0,90,1000,0,350)",
tree->Draw("fRC_Kinetic_Energy[0]:fRC_Theta[0]>>hETheta_MultOuterT(900,0,90,1000,0,350)","Strasse.GetOuterMultTEnergy()==1","colz"); "Strasse.GetOuterMultTEnergy()==1", "colz");
c3->cd(5); c3->cd(5);
tree->Draw("fRC_Kinetic_Energy[0]:fRC_Theta[0]>>hETheta_MultOuterL(900,0,90,1000,0,350)","Strasse.GetOuterMultLEnergy()==1","colz"); tree->Draw("fRC_Kinetic_Energy[0]:fRC_Theta[0]>>hETheta_MultOuterL(900,0,90,1000,0,350)",
c3->cd(6); "Strasse.GetOuterMultLEnergy()==1", "colz");
tree->Draw("fRC_Kinetic_Energy[0]:fRC_Theta[0]>>hETheta_MultAll(900,0,90,1000,0,350)","Strasse.GetInnerMultTEnergy()==1 && Strasse.GetInnerMultLEnergy()==1 && Strasse.GetOuterMultLEnergy()==1 && Strasse.GetOuterMultTEnergy()==1","colz"); c3->cd(6);
tree->Draw("fRC_Kinetic_Energy[0]:fRC_Theta[0]>>hETheta_MultAll(900,0,90,1000,0,350)",
*/ "Strasse.GetInnerMultTEnergy()==1 && Strasse.GetInnerMultLEnergy()==1 && Strasse.GetOuterMultLEnergy()==1 "
"&& Strasse.GetOuterMultTEnergy()==1",
"colz");
} }
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