diff --git a/Inputs/EventGenerator/isotropic.source b/Inputs/EventGenerator/isotropic.source
index 2d4cc8a0742895e6e27db4eb278195f457d9307c..623c8a24580daf12d48752e5474d2c69c63c8c63 100644
--- a/Inputs/EventGenerator/isotropic.source
+++ b/Inputs/EventGenerator/isotropic.source
@@ -4,13 +4,13 @@
% Energy are given in MeV , Position in mm %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Isotropic
- EnergyLow= 100
- EnergyHigh= 100
- HalfOpenAngleMin= 40
- HalfOpenAngleMax= 40
+ EnergyLow= 0
+ EnergyHigh= 200
+ HalfOpenAngleMin= 30
+ HalfOpenAngleMax= 60
x0= 0
y0= 0
z0= 0
- particle= 14C
+ particle= 11Be
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Supported particle type: proton, neutron, deuton, triton, He3 , alpha
diff --git a/NPLib/Detectors/Strasse/TStrassePhysics.cxx b/NPLib/Detectors/Strasse/TStrassePhysics.cxx
index 8468b3172d313840d723de5f2fd6f17affbf2a53..fa2f77d9459bf0bd9ee1791c434492a1afa7f393 100644
--- a/NPLib/Detectors/Strasse/TStrassePhysics.cxx
+++ b/NPLib/Detectors/Strasse/TStrassePhysics.cxx
@@ -55,7 +55,7 @@ ClassImp(TStrassePhysics)
m_NumberOfInnerDetectors = 0;
m_NumberOfOuterDetectors = 0;
m_MaximumStripMultiplicityAllowed = 10;
- m_StripEnergyMatching = 1.00;
+ m_StripEnergyMatching = 0.100;
////////////////////
// Inner Detector //
@@ -369,32 +369,37 @@ void TStrassePhysics::BuildPhysicalEvent() {
int innerT = m_PreTreatedData->GetInner_TE_StripNbr(inner[i].X());
int innerL = m_PreTreatedData->GetInner_LE_StripNbr(inner[i].Y());
- double TE = m_PreTreatedData->GetInner_TE_Energy(inner[i].X());
+ double innerTE = m_PreTreatedData->GetInner_TE_Energy(inner[i].X());
+ double innerLE = m_PreTreatedData->GetInner_LE_Energy(inner[i].X());
// look for outer
- double outerE = 0;
+ double outerTE = 0;
+ double outerLE = 0;
int outerT=0;
int outerL=0;
for(unsigned int j=0; j<outer.size(); j++){
if(m_PreTreatedData->GetInner_TE_DetectorNbr(outer[j].X())==N){
- outerE = m_PreTreatedData->GetOuter_TE_Energy(outer[j].X());
+ outerTE = m_PreTreatedData->GetOuter_TE_Energy(outer[j].X());
+ outerLE = m_PreTreatedData->GetOuter_LE_Energy(outer[j].X());
outerT = m_PreTreatedData->GetOuter_TE_StripNbr(outer[j].X());
outerL = m_PreTreatedData->GetOuter_LE_StripNbr(outer[j].Y());
}
}
- if(outerE){
+ if(outerTE){
EventMultiplicity++;
DetectorNumber.push_back(N);
InnerStripT.push_back(innerT);
InnerStripL.push_back(innerL);
- DE.push_back(TE);
+ InnerTE.push_back(innerTE);
+ InnerLE.push_back(innerLE);
InnerPosX.push_back(GetInnerPositionOfInteraction(EventMultiplicity-1).x());
InnerPosY.push_back(GetInnerPositionOfInteraction(EventMultiplicity-1).y());
InnerPosZ.push_back(GetInnerPositionOfInteraction(EventMultiplicity-1).z());
OuterStripT.push_back(outerT);
OuterStripL.push_back(outerL);
- E.push_back(outerE);
+ OuterTE.push_back(outerTE);
+ OuterLE.push_back(outerLE);
OuterPosX.push_back(GetOuterPositionOfInteraction(EventMultiplicity-1).x());
OuterPosY.push_back(GetOuterPositionOfInteraction(EventMultiplicity-1).y());
OuterPosZ.push_back(GetOuterPositionOfInteraction(EventMultiplicity-1).z());
@@ -604,6 +609,12 @@ void TStrassePhysics::ReadAnalysisConfig() {
cout << whatToDo << " " << m_E_Threshold << endl;
}
+ else if (whatToDo=="E_FRONTBACK_MATCHING") {
+ AnalysisConfigFile >> DataBuffer;
+ m_StripEnergyMatching = atof(DataBuffer.c_str());
+ cout << whatToDo << " " << m_StripEnergyMatching << endl;
+ }
+
else {
ReadingStatus = false;
}
@@ -618,12 +629,14 @@ void TStrassePhysics::Clear() {
EventMultiplicity = 0;
// DSSD
DetectorNumber.clear();
- E.clear();
InnerStripT.clear();
InnerStripL.clear();
OuterStripT.clear();
OuterStripL.clear();
- DE.clear();
+ InnerTE.clear();
+ OuterTE.clear();
+ InnerLE.clear();
+ OuterLE.clear();
// Position Information
InnerPosX.clear();
@@ -633,7 +646,6 @@ void TStrassePhysics::Clear() {
OuterPosY.clear();
OuterPosZ.clear();
-
}
diff --git a/NPLib/Detectors/Strasse/TStrassePhysics.h b/NPLib/Detectors/Strasse/TStrassePhysics.h
index 49e035e051230fd4a9dad00c1cf3c7f5975cb15b..c220b7bda3a52147cd4264c1da58e971f5cc3b3d 100644
--- a/NPLib/Detectors/Strasse/TStrassePhysics.h
+++ b/NPLib/Detectors/Strasse/TStrassePhysics.h
@@ -68,13 +68,14 @@ class TStrassePhysics : public TObject, public NPL::VDetector {
public:
Int_t EventMultiplicity;
vector<int> DetectorNumber;
- vector<double> E;
- vector<double> DE;
vector<int> InnerStripT;
vector<int> InnerStripL;
+ vector<double> InnerTE;
+ vector<double> InnerLE;
vector<int> OuterStripT;
vector<int> OuterStripL;
-
+ vector<double> OuterTE;
+ vector<double> OuterLE;
vector<double> InnerPosX;
vector<double> InnerPosY;
@@ -82,8 +83,6 @@ class TStrassePhysics : public TObject, public NPL::VDetector {
vector<double> OuterPosX;
vector<double> OuterPosY;
vector<double> OuterPosZ;
-
-
//////////////////////////////////////////////////////////////
// methods inherited from the VDetector ABC class
@@ -164,6 +163,9 @@ class TStrassePhysics : public TObject, public NPL::VDetector {
double GetNumberOfInnerDetector() const {return m_NumberOfInnerDetectors;}
double GetNumberOfOuterDetector() const {return m_NumberOfOuterDetectors;}
int GetEventMultiplicity() const {return EventMultiplicity;}
+ int GetDetNbrSize() const {return DetectorNumber.size();}
+ TVector3 GetInnerPos(const int i) const {return TVector3(InnerPosX[i],InnerPosY[i],InnerPosZ[i]);}
+ TVector3 GetOuterPos(const int i) const {return TVector3(OuterPosX[i],OuterPosY[i],OuterPosZ[i]);}
double GetInnerStripPositionX(const int N, const int X, const int Y){
return m_InnerStripPositionX[N-1][X-1][Y-1];
diff --git a/NPLib/Detectors/Strasse/TStrasseSpectra.cxx b/NPLib/Detectors/Strasse/TStrasseSpectra.cxx
index 64edd19ce58d7f06e918a5132c49cfc8881437fb..e49d11876fa9fd4697664bb5a8cb72f7b173c473 100644
--- a/NPLib/Detectors/Strasse/TStrasseSpectra.cxx
+++ b/NPLib/Detectors/Strasse/TStrasseSpectra.cxx
@@ -140,8 +140,8 @@ void TStrasseSpectra::FillPhysicsSpectra(TStrassePhysics* Physics) {
family= "Strasse/PHY";
// Energy vs time
- unsigned int sizeE = Physics->E.size();
- for(unsigned int i = 0 ; i < sizeE ; i++){
+ unsigned int sizeTE = Physics->InnerTE.size();
+ for(unsigned int i = 0 ; i < sizeTE ; i++){
//name = "Strasse_ENERGY_TIME";
//FillSpectra(family,name,Physics->E[i],Physics->Time[i]);
}
diff --git a/NPLib/Physics/NPParticle.cxx b/NPLib/Physics/NPParticle.cxx
index 733e3a8731e8f2f8fdc59ad578007bb6e98895bf..c34b283b7af0ce0ea58342de7ea3c9c0a383a12e 100644
--- a/NPLib/Physics/NPParticle.cxx
+++ b/NPLib/Physics/NPParticle.cxx
@@ -258,14 +258,16 @@ Particle::Particle(int Z, int A)
{
//----------- Constructor Using nubtab12.asc by Z and A----------
- if(Z==0 && A==4){
- SetUp("4n");
- return;
- }
+ fExcitationEnergy = 0;
+
+ if(Z==0 && A==4){
+ SetUp("4n");
+ return;
+ }
// open the file to read and check if it is open
ifstream inFile;
string Path = getenv("NPTOOL") ;
- string FileName = Path + "/NPLib/Physics/nubtab12.asc";
+ string FileName = Path + "/NPLib/Physics/nubtab16.asc";
inFile.open(FileName.c_str());
// reading the file
@@ -444,7 +446,7 @@ void Particle::Extract(string line){
fLifeTime*=3600*24*365.25*1e12;
fLifeTimeErr*=3600*24*365.25*1e12;
}
-
+
// spin and parity
string s_spinparity = line.substr(79,14);
fSpinParity = s_spinparity.data();
@@ -493,14 +495,14 @@ void Particle::Print() const
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void Particle::GetParticleName() {
- fParticleName.assign(fName);
+ fParticleName.assign(fName);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void Particle::EnergyToBrho(double Q){
if(Q==-1000)
- Q=GetZ();
-
+ Q=GetZ();
+
EnergyToBeta();
BetaToGamma();
//fBrho = sqrt(pow(fKineticEnergy,2) + 2*fKineticEnergy*Mass()) * 1e6 * e_SI / (c_light*1e6) / (Q * e_SI);
@@ -522,7 +524,7 @@ void Particle::TofToEnergy() {
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void Particle::BrhoToEnergy(double Q ){
- if(Q==-1000)
+ if(Q==-1000)
Q=GetZ();
fKineticEnergy = sqrt( pow((fBrho*(c_light*1e6)*Q*e_SI)/(1e6*e_SI),2) + pow(Mass(),2) ) - Mass();
@@ -559,60 +561,60 @@ double Particle::DopplerCorrection(double EnergyLabGamma, double ThetaLabGamma){
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
double Particle::GetEnergyCM(double EnergyLab, double ThetaLab, double PhiLab, double relativisticboost){
- SetKineticEnergy(EnergyLab);
- double EnergyCM;
- double ImpulsionLab;
- double ImpulsionLabX;
- double ImpulsionLabY;
- double ImpulsionLabZ;
- TVector3 VImpulsionLAB;
- TLorentzVector LVEnergyImpulsionLAB;
- TLorentzVector LVEnergyImpulsionCM;
-
- ImpulsionLab = sqrt(EnergyLab*EnergyLab + 2*EnergyLab*Mass());
- ImpulsionLabX = ImpulsionLab*sin(ThetaLab)*cos(PhiLab);
- ImpulsionLabY = ImpulsionLab*sin(ThetaLab)*sin(PhiLab);
- ImpulsionLabZ = ImpulsionLab*cos(ThetaLab);
-
- VImpulsionLAB = TVector3(ImpulsionLabX, ImpulsionLabY, ImpulsionLabZ);
- LVEnergyImpulsionLAB = TLorentzVector(VImpulsionLAB,Mass()+EnergyLab);
-
- LVEnergyImpulsionCM = LVEnergyImpulsionLAB;
- LVEnergyImpulsionCM.Boost(0,0,-relativisticboost);
-
- EnergyCM = LVEnergyImpulsionCM.Energy() - Mass();
-
- return EnergyCM;
+ SetKineticEnergy(EnergyLab);
+ double EnergyCM;
+ double ImpulsionLab;
+ double ImpulsionLabX;
+ double ImpulsionLabY;
+ double ImpulsionLabZ;
+ TVector3 VImpulsionLAB;
+ TLorentzVector LVEnergyImpulsionLAB;
+ TLorentzVector LVEnergyImpulsionCM;
+
+ ImpulsionLab = sqrt(EnergyLab*EnergyLab + 2*EnergyLab*Mass());
+ ImpulsionLabX = ImpulsionLab*sin(ThetaLab)*cos(PhiLab);
+ ImpulsionLabY = ImpulsionLab*sin(ThetaLab)*sin(PhiLab);
+ ImpulsionLabZ = ImpulsionLab*cos(ThetaLab);
+
+ VImpulsionLAB = TVector3(ImpulsionLabX, ImpulsionLabY, ImpulsionLabZ);
+ LVEnergyImpulsionLAB = TLorentzVector(VImpulsionLAB,Mass()+EnergyLab);
+
+ LVEnergyImpulsionCM = LVEnergyImpulsionLAB;
+ LVEnergyImpulsionCM.Boost(0,0,-relativisticboost);
+
+ EnergyCM = LVEnergyImpulsionCM.Energy() - Mass();
+
+ return EnergyCM;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
double Particle::GetThetaCM(double EnergyLab, double ThetaLab, double PhiLab, double relativisticboost){
- SetKineticEnergy(EnergyLab);
- double EnergyCM;
- double ThetaCM;
- double ImpulsionLab;
- double ImpulsionLabX;
- double ImpulsionLabY;
- double ImpulsionLabZ;
- TVector3 VImpulsionLAB;
- TLorentzVector LVEnergyImpulsionLAB;
- TLorentzVector LVEnergyImpulsionCM;
-
- ImpulsionLab = sqrt(EnergyLab*EnergyLab + 2*EnergyLab*Mass());
- ImpulsionLabX = ImpulsionLab*sin(ThetaLab)*cos(PhiLab);
- ImpulsionLabY = ImpulsionLab*sin(ThetaLab)*sin(PhiLab);
- ImpulsionLabZ = ImpulsionLab*cos(ThetaLab);
-
- VImpulsionLAB = TVector3(ImpulsionLabX, ImpulsionLabY, ImpulsionLabZ);
- LVEnergyImpulsionLAB = TLorentzVector(VImpulsionLAB,Mass()+EnergyLab);
-
- LVEnergyImpulsionCM = LVEnergyImpulsionLAB;
- LVEnergyImpulsionCM.Boost(0,0,-relativisticboost);
-
- EnergyCM = LVEnergyImpulsionCM.Energy() - Mass();
- ThetaCM = LVEnergyImpulsionCM.Theta();
-
- return ThetaCM;
+ SetKineticEnergy(EnergyLab);
+ double EnergyCM;
+ double ThetaCM;
+ double ImpulsionLab;
+ double ImpulsionLabX;
+ double ImpulsionLabY;
+ double ImpulsionLabZ;
+ TVector3 VImpulsionLAB;
+ TLorentzVector LVEnergyImpulsionLAB;
+ TLorentzVector LVEnergyImpulsionCM;
+
+ ImpulsionLab = sqrt(EnergyLab*EnergyLab + 2*EnergyLab*Mass());
+ ImpulsionLabX = ImpulsionLab*sin(ThetaLab)*cos(PhiLab);
+ ImpulsionLabY = ImpulsionLab*sin(ThetaLab)*sin(PhiLab);
+ ImpulsionLabZ = ImpulsionLab*cos(ThetaLab);
+
+ VImpulsionLAB = TVector3(ImpulsionLabX, ImpulsionLabY, ImpulsionLabZ);
+ LVEnergyImpulsionLAB = TLorentzVector(VImpulsionLAB,Mass()+EnergyLab);
+
+ LVEnergyImpulsionCM = LVEnergyImpulsionLAB;
+ LVEnergyImpulsionCM.Boost(0,0,-relativisticboost);
+
+ EnergyCM = LVEnergyImpulsionCM.Energy() - Mass();
+ ThetaCM = LVEnergyImpulsionCM.Theta();
+
+ return ThetaCM;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -651,58 +653,58 @@ void Particle::DefineMassByThreshold(const vector<NPL::Particle>& N){
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
double Particle::GetSXn(unsigned int X) const {
- Particle N(GetZ(),GetA()-X);
+ Particle N(GetZ(),GetA()-X);
return N.Mass()+X*neutron_mass_c2-Mass();
- }
+}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
double Particle::GetSXp(unsigned int X) const {
- Particle N(GetZ()-X,GetA()-X);
+ Particle N(GetZ()-X,GetA()-X);
return N.Mass()+X*proton_mass_c2-Mass();
- }
+}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
double Particle::GetSn() const {
return GetSXn(1);
- }
+}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
double Particle::GetSp() const {
return GetSXp(1);
- }
+}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
double Particle::GetS2n() const {
return GetSXn(2);
- }
+}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
double Particle::GetS2p() const {
return GetSXp(2);
- }
+}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
double Particle::GetSt() const {
- Particle N(GetZ()-1,GetA()-3);
- Particle A(1,3);
+ Particle N(GetZ()-1,GetA()-3);
+ Particle A(1,3);
return N.Mass()+A.Mass()-Mass();
- }
+}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
double Particle::GetS3He() const {
- Particle N(GetZ()-2,GetA()-3);
- Particle A(2,3);
+ Particle N(GetZ()-2,GetA()-3);
+ Particle A(2,3);
return N.Mass()+A.Mass()-Mass();
- }
+}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
double Particle::GetSa() const {
- Particle N(GetZ()-2,GetA()-4);
- Particle A(2,4);
+ Particle N(GetZ()-2,GetA()-4);
+ Particle A(2,4);
return N.Mass()+A.Mass()-Mass();
- }
+}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void Particle::PrintThreshold() const {
- cout << GetName() << " thresholds : " << endl;
- cout << " Sn : " << GetSn() << " MeV" << endl;
- cout << " Sp : " << GetSp() << " MeV" << endl;
- cout << " S2n : " << GetS2n() << " MeV" << endl;
- cout << " S2p : " << GetS2p() << " MeV" << endl;
- cout << " St : " << GetSt() << " MeV" << endl;
- cout << " S3He : " << GetS3He() << " MeV" << endl;
- cout << " Sa : " << GetSa() << " MeV" << endl;
- }
+ cout << GetName() << " thresholds : " << endl;
+ cout << " Sn : " << GetSn() << " MeV" << endl;
+ cout << " Sp : " << GetSp() << " MeV" << endl;
+ cout << " S2n : " << GetS2n() << " MeV" << endl;
+ cout << " S2p : " << GetS2p() << " MeV" << endl;
+ cout << " St : " << GetSt() << " MeV" << endl;
+ cout << " S3He : " << GetS3He() << " MeV" << endl;
+ cout << " Sa : " << GetSa() << " MeV" << endl;
+}
diff --git a/NPSimulation/Detectors/PISTA/PISTA.cc b/NPSimulation/Detectors/PISTA/PISTA.cc
index 89cbc8e7c8db033c79ea3f033ab8850b44af760c..5a22806e550132c67f119a3de49f1b1b4b58ae63 100644
--- a/NPSimulation/Detectors/PISTA/PISTA.cc
+++ b/NPSimulation/Detectors/PISTA/PISTA.cc
@@ -60,8 +60,8 @@ namespace PISTA_NS{
// Energy and time Resolution
const double EnergyThreshold = 0.1*MeV;
const double ResoTime = 0.2*ns ;
- const double ResoEnergy = 0.015*MeV ;
- const double DE_ResoEnergy = 0.015*MeV ;
+ const double DE_ResoEnergy = 0.040*MeV ;
+ const double E_ResoEnergy = 0.018*MeV ;
// Trapezoid dimension
const double TrapezoidBaseLarge = 74.1*mm;
@@ -282,7 +282,7 @@ void PISTA::ReadSensitive(const G4Event* ){
if(Energy>EnergyThreshold){
double Time = RandGauss::shoot(FirstStageScorer->GetTimeLength(i), ResoTime);
int DetNbr = FirstStageScorer->GetDetectorLength(i);
- int StripFront = FirstStageScorer->GetStripLength(i);
+ int StripFront = FirstStageScorer->GetStripWidth(i);
m_Event->SetPISTA_DE(DetNbr, StripFront, Energy, Energy, Time, Time);
}
}
@@ -294,7 +294,7 @@ void PISTA::ReadSensitive(const G4Event* ){
unsigned int sizeE = SecondStageScorer->GetLengthMult();
for(unsigned int i = 0 ; i < sizeE ; i++){
- double Energy = RandGauss::shoot(SecondStageScorer->GetEnergyLength(i), ResoEnergy);
+ double Energy = RandGauss::shoot(SecondStageScorer->GetEnergyLength(i), E_ResoEnergy);
if(Energy>EnergyThreshold){
double Time = RandGauss::shoot(SecondStageScorer->GetTimeLength(i), ResoTime);
int DetNbr = SecondStageScorer->GetDetectorLength(i);
diff --git a/NPSimulation/Detectors/Strasse/Strasse.cc b/NPSimulation/Detectors/Strasse/Strasse.cc
index 230a0e44e20d9ab8989a5372cc79150ee873d9f6..bd4061cd8863caa01a041d21fba9ee9e57707e59 100644
--- a/NPSimulation/Detectors/Strasse/Strasse.cc
+++ b/NPSimulation/Detectors/Strasse/Strasse.cc
@@ -64,7 +64,7 @@ using namespace CLHEP;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
namespace Strasse_NS {
// Energy and time Resolution
- const double EnergyThreshold = 10 * keV;
+ const double EnergyThreshold = 1 * keV;
const double ResoEnergy = 0.015 * MeV;
////////////////////
@@ -139,12 +139,15 @@ using namespace Strasse_NS;
Strasse::Strasse() {
InitializeMaterial();
m_Event = new TStrasseData();
+ m_ReactionRegion = NULL;
m_InnerScorer1 = 0;
m_OuterScorer1 = 0;
m_InnerScorer2 = 0;
m_OuterScorer2 = 0;
m_InnerDetector = 0;
m_OuterDetector = 0;
+ m_Target = 0;
+ m_TargetCell = 0;
m_Chamber = 0;
m_Blades = 0;
m_Stars = 0;
@@ -168,6 +171,9 @@ Strasse::Strasse() {
GuardRingVisAtt = new G4VisAttributes(G4Colour(0.85, 0.85, 0.85, 0.5));
// Light Blue
BladeVisAtt = new G4VisAttributes(G4Colour(1, 0.65, 0.0, 0.7));
+ // Target: Transparent grey for the cell + Transparent blue for the LH2
+ TargetVisAtt = new G4VisAttributes(G4Colour(0.15, 0.85, 0.85, 0.1));
+ TargetCellVisAtt = new G4VisAttributes(G4Colour(0.8, 0.8, 0.8, 0.5));
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -194,6 +200,16 @@ void Strasse::AddOuterDetector(double R, double Z, double Phi, double Shift, G4T
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void Strasse::AddChamber(double Z) { m_Chamber_Z.push_back(Z); }
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void Strasse::AddTarget(double R, double L, string MaterialName, string CellMaterialName, double CellThickness,
+ G4ThreeVector Pos) {
+ m_Target_R.push_back(R);
+ m_Target_L.push_back(L);
+ m_Target_MaterialName.push_back(MaterialName);
+ m_Target_CellMaterialName.push_back(CellMaterialName);
+ m_Target_CellThickness.push_back(CellThickness);
+ m_Target_Pos.push_back(Pos);
+}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4LogicalVolume* Strasse::BuildInnerDetector() {
if (!m_InnerDetector) {
@@ -549,6 +565,51 @@ G4LogicalVolume* Strasse::BuildOuterDetector() {
return m_OuterDetector;
}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+G4LogicalVolume* Strasse::BuildTarget(int i) {
+ // NOT YET FUNCTIONAL: multiple target blocks
+ if (i > 0) {
+ cout << "ERROR: Multiple STRASSE target blocks defined in detector file" << endl;
+ }
+
+ G4Material* TargetMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary(m_Target_MaterialName[i]);
+ G4Tubs* solidTarget = new G4Tubs("Target", // its name
+ 0., m_Target_R[i], m_Target_L[i] / 2., 0, 360.); // size
+ m_Target = new G4LogicalVolume(solidTarget, // its solid
+ TargetMaterial, // its material
+ "Target"); // its name
+ m_Target->SetVisAttributes(TargetVisAtt);
+
+ return m_Target;
+}
+
+G4LogicalVolume* Strasse::BuildTargetCell(int i) {
+ // NOT YET FUNCTIONAL: multiple target blocks
+ if (i > 0) {
+ cout << "ERROR: Multiple STRASSE target blocks defined in detector file" << endl;
+ }
+
+ // if (!m_TargetCell && !m_Target) {
+ // Start with the full target cell volume
+ G4Material* CellMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary(m_Target_CellMaterialName[i]);
+ G4Tubs* CellFull =
+ new G4Tubs("CellFull", // its name
+ 0., m_Target_R[i] + m_Target_CellThickness[i], m_Target_L[i] / 2 + m_Target_CellThickness[i], 0, 360.);
+ // For subtraction of internal part, to only keep the cell
+ G4Tubs* CellInterior = new G4Tubs("CellInterior", // its name
+ 0., m_Target_R[i], m_Target_L[i] / 2., 0, 360.);
+ // Cell creation
+ G4SubtractionSolid* solidCell = new G4SubtractionSolid("solidCell", CellFull, CellInterior,
+ new G4RotationMatrix(0, 0, 0), G4ThreeVector(0, 0, 0));
+
+ m_TargetCell = new G4LogicalVolume(solidCell, // its solid
+ CellMaterial, // its material
+ "TargetCell"); // name
+ m_TargetCell->SetVisAttributes(TargetCellVisAtt);
+
+ return m_TargetCell;
+}
+
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4LogicalVolume* Strasse::BuildChamber() {
if (!m_Chamber) {
@@ -789,6 +850,32 @@ void Strasse::ReadConfiguration(NPL::InputParser parser) {
exit(1);
}
+ // Strasse Custom Target Block
+ vector<NPL::InputBlock*> blocks_target = parser.GetAllBlocksWithTokenAndValue("Strasse", "Target");
+ if (NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << "//// " << blocks_target.size() << " Custom strasse target found " << endl;
+
+ vector<string> targettoken = {"Radius", "Length", "TargetMaterial", "CellMaterial", "CellThickness", "Pos"};
+
+ for (unsigned int i = 0; i < blocks_target.size(); i++) {
+ if (blocks_target[i]->HasTokenList(targettoken)) {
+ if (NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << endl << "//// Strasse custom target" << i + 1 << endl;
+
+ double R = blocks_target[i]->GetDouble("Radius", "mm");
+ double L = blocks_target[i]->GetDouble("Length", "mm");
+ string TargetMaterialName = blocks_target[i]->GetString("TargetMaterial");
+ string CellMaterialName = blocks_target[i]->GetString("CellMaterial");
+ double CellThickness = blocks_target[i]->GetDouble("CellThickness", "mm");
+ G4ThreeVector Pos = NPS::ConvertVector(blocks_target[i]->GetTVector3("Pos", "mm"));
+ AddTarget(R, L, TargetMaterialName, CellMaterialName, CellThickness, Pos);
+ }
+ else {
+ cout << "ERROR: check your input file formatting on " << i + 1 << " custom strasse target block " << endl;
+ exit(1);
+ }
+ }
+
// Inner Barrel
vector<NPL::InputBlock*> blocks_inner = parser.GetAllBlocksWithTokenAndValue("Strasse", "Inner");
if (NPOptionManager::getInstance()->GetVerboseLevel())
@@ -916,7 +1003,6 @@ void Strasse::ConstructDetector(G4LogicalVolume* world) {
}
// Chamber
-
for (unsigned short i = 0; i < m_Chamber_Z.size(); i++) {
G4ThreeVector Det_pos = G4ThreeVector(0, 0, -m_Chamber_Z[i]);
G4RotationMatrix* Rot = new G4RotationMatrix();
@@ -924,6 +1010,42 @@ void Strasse::ConstructDetector(G4LogicalVolume* world) {
new G4PVPlacement(G4Transform3D(*Rot, Det_pos), BuildChamber(), "Strasse", world, false, i + 1);
}
+ // Active Target volume (LH2 only)
+ G4LogicalVolume* logicTarget[m_Target_R.size()];
+ G4LogicalVolume* logicTargetCell[m_Target_R.size()];
+ for (unsigned short i = 0; i < m_Target_R.size(); i++) {
+ G4ThreeVector Tar_pos = m_Target_Pos[i];
+ cout << "TargetPos:" << m_Target_Pos[i].z() << endl;
+ G4RotationMatrix* Rot = new G4RotationMatrix();
+ logicTarget[i] = BuildTarget(i);
+ new G4PVPlacement(Rot, Tar_pos, logicTarget[i], "Strasse_Target", world, false, i + 1);
+ logicTargetCell[i] = BuildTargetCell(i);
+ new G4PVPlacement(Rot, Tar_pos, logicTargetCell[i], "Strasse_TargetCell", world, false, i + 1);
+
+ if (!m_ReactionRegion) {
+ m_ReactionRegion = new G4Region("NPSimulationProcess");
+ }
+ // Reaction only possible LH2 volume, not in the cell window
+ // so use "m_Target" defined in BuildTarget() and not logicTarget[i]/m_TargetCell
+ m_ReactionRegion->AddRootLogicalVolume(m_Target);
+ m_ReactionRegion->SetUserLimits(new G4UserLimits(1. * mm)); //???
+
+ G4FastSimulationManager* mng = m_ReactionRegion->GetFastSimulationManager();
+ unsigned int size = m_ReactionModel.size();
+ for (unsigned int o = 0; o < size; o++) {
+ mng->RemoveFastSimulationModel(m_ReactionModel[o]);
+ }
+ m_ReactionModel.clear();
+
+ G4VFastSimulationModel* fsm;
+ fsm = new NPS::BeamReaction("BeamReaction", m_ReactionRegion);
+ ((NPS::BeamReaction*)fsm)->SetStepSize(1. * mm);
+ m_ReactionModel.push_back(fsm);
+
+ fsm = new NPS::Decay("Decay", m_ReactionRegion);
+ m_ReactionModel.push_back(fsm);
+ }
+
// G4ThreeVector Det_pos = G4ThreeVector(0,0,+11.5) ;
G4ThreeVector Det_pos = G4ThreeVector(0, 0, 0);
G4RotationMatrix* Rot = new G4RotationMatrix();
@@ -933,7 +1055,6 @@ void Strasse::ConstructDetector(G4LogicalVolume* world) {
if (found_chamber) {
new G4PVPlacement(Rot, Det_pos, BuildChamberFromCAD(ChamberPath), "Strasse_Chamber", world, false, 0);
}
-
if (found_blades) {
new G4PVPlacement(Rot, Det_pos, BuildBlades(BladesPath), "Strasse_Blades", world, false, 0);
}
diff --git a/NPSimulation/Detectors/Strasse/Strasse.hh b/NPSimulation/Detectors/Strasse/Strasse.hh
index 0250f105e330526b09a2647d71867cc38f665c5c..76ab02f78f2b6babb3249580240dedb295c40fa2 100644
--- a/NPSimulation/Detectors/Strasse/Strasse.hh
+++ b/NPSimulation/Detectors/Strasse/Strasse.hh
@@ -27,149 +27,167 @@
using namespace std;
// G4 headers
-#include "G4ThreeVector.hh"
-#include "G4RotationMatrix.hh"
+#include "G4FastSimulationManager.hh"
#include "G4LogicalVolume.hh"
#include "G4MultiFunctionalDetector.hh"
+#include "G4RotationMatrix.hh"
+#include "G4ThreeVector.hh"
+#include "G4UserLimits.hh"
+#include "G4VFastSimulationModel.hh"
// NPTool header
+#include "BeamReaction.hh"
+#include "Decay.hh"
+#include "NPInputParser.h"
#include "NPSVDetector.hh"
#include "TStrasseData.h"
-#include "NPInputParser.h"
-class Strasse : public NPS::VDetector{
+class Strasse : public NPS::VDetector {
////////////////////////////////////////////////////
/////// Default Constructor and Destructor /////////
////////////////////////////////////////////////////
- public:
- Strasse() ;
- virtual ~Strasse() ;
-
- ////////////////////////////////////////////////////
- /////// Specific Function of this Class ///////////
- ////////////////////////////////////////////////////
- public:
- // Cylindrical coordinate
- void AddInnerDetector(double R,double Z,double Phi, double Shift, G4ThreeVector Ref);
- void AddOuterDetector(double R,double Z,double Phi, double Shift, G4ThreeVector Ref);
- void AddChamber(double Z);
-
- G4LogicalVolume* BuildInnerDetector();
- G4LogicalVolume* BuildOuterDetector();
- G4LogicalVolume* BuildElectronic();
- G4LogicalVolume* BuildChamber();
- G4LogicalVolume* BuildChamberFromCAD(string path);
- G4LogicalVolume* BuildStars(string path);
- G4LogicalVolume* BuildBlades(string path);
- G4LogicalVolume* BuildBase(string path);
-
- private:
- G4LogicalVolume* m_InnerDetector;
- G4LogicalVolume* m_OuterDetector;
- G4LogicalVolume* m_Electronic;
- G4LogicalVolume* m_Stars;
- G4LogicalVolume* m_Chamber;
- G4LogicalVolume* m_Blades;
- G4LogicalVolume* m_Base;
-
- string ChamberPath;
- string BasePath;
- string StarsPath;
- string BladesPath;
- bool found_chamber;
- bool found_blades;
- bool found_stars;
- bool found_base;
-
- private:
- // Initialize material used in detector definition
- void InitializeMaterial();
-
-
- // List of material
- G4Material* m_MaterialSilicon ;
- G4Material* m_MaterialAl ;
- G4Material* m_MaterialVacuum ;
- G4Material* m_MaterialPCB ;
- G4Material* m_MaterialCu ;
-
- // calculated dimension
- double m_Active_InnerWafer_Width;
- double m_Active_InnerWafer_Length;
- double m_Active_OuterWafer_Width;
- double m_Active_OuterWafer_Length;
-
-
- ////////////////////////////////////////////////////
- ////// Inherite from NPS::VDetector class /////////
- ////////////////////////////////////////////////////
- public:
- // Read stream at Configfile to pick-up parameters of detector (Position,...)
- // Called in DetecorConstruction::ReadDetextorConfiguration Method
- void ReadConfiguration(NPL::InputParser) ;
-
- // Construct detector and inialise sensitive part.
- // Called After DetecorConstruction::AddDetector Method
- void ConstructDetector(G4LogicalVolume* world) ;
-
- // Add Detector branch to the EventTree.
- // Called After DetecorConstruction::AddDetector Method
- void InitializeRootOutput() ;
-
- // Read sensitive part and fill the Root tree.
- // Called at in the EventAction::EndOfEventAvtion
- void ReadSensitive(const G4Event* event) ;
-
- public: // Scorer
- // Initialize all Scorer used by the MUST2Array
- void InitializeScorers() ;
-
- // Associated Scorer
- G4MultiFunctionalDetector* m_InnerScorer1 ;
- G4MultiFunctionalDetector* m_OuterScorer1 ;
- G4MultiFunctionalDetector* m_InnerScorer2 ;
- G4MultiFunctionalDetector* m_OuterScorer2 ;
-
- ////////////////////////////////////////////////////
- ///////////Event class to store Data////////////////
- ////////////////////////////////////////////////////
- private:
- TStrasseData* m_Event ;
-
- ////////////////////////////////////////////////////
- ///////////////Private intern Data//////////////////
- ////////////////////////////////////////////////////
- private: // Geometry
- // Detector Coordinate
- vector<double> m_Inner_R;
- vector<double> m_Inner_Z;
- vector<double> m_Inner_Phi;
- vector<double> m_Inner_Shift;
- vector<G4ThreeVector> m_Inner_Ref;
-
- vector<double> m_Outer_R;
- vector<double> m_Outer_Z;
- vector<double> m_Outer_Phi;
- vector<double> m_Outer_Shift;
- vector<G4ThreeVector> m_Outer_Ref;
-
- vector<double> m_Chamber_Z;
-
-
- // Needed for dynamic loading of the library
- public:
- static NPS::VDetector* Construct();
-
-
- private: // Visualisation
- G4VisAttributes* SiliconVisAtt ;
- G4VisAttributes* PCBVisAtt;
- G4VisAttributes* PADVisAtt ;
- G4VisAttributes* StarsVisAtt ;
- G4VisAttributes* ChamberVisAtt ;
- G4VisAttributes* GuardRingVisAtt ;
- G4VisAttributes* BladeVisAtt ;
+ public:
+ Strasse();
+ virtual ~Strasse();
+ ////////////////////////////////////////////////////
+ /////// Specific Function of this Class ///////////
+ ////////////////////////////////////////////////////
+ public:
+ // Cylindrical coordinate
+ void AddInnerDetector(double R, double Z, double Phi, double Shift, G4ThreeVector Ref);
+ void AddOuterDetector(double R, double Z, double Phi, double Shift, G4ThreeVector Ref);
+ void AddTarget(double R, double L, string MaterialName, string CellMaterialName, double CellThickness,
+ G4ThreeVector Pos);
+ void AddChamber(double Z);
+
+ G4LogicalVolume* BuildInnerDetector();
+ G4LogicalVolume* BuildOuterDetector();
+ G4LogicalVolume* BuildTarget(int i);
+ G4LogicalVolume* BuildTargetCell(int i);
+ G4LogicalVolume* BuildElectronic();
+ G4LogicalVolume* BuildChamber();
+ G4LogicalVolume* BuildChamberFromCAD(string path);
+ G4LogicalVolume* BuildStars(string path);
+ G4LogicalVolume* BuildBlades(string path);
+ G4LogicalVolume* BuildBase(string path);
+
+ private:
+ G4LogicalVolume* m_InnerDetector;
+ G4LogicalVolume* m_OuterDetector;
+ G4LogicalVolume* m_Target;
+ G4LogicalVolume* m_TargetCell;
+ G4LogicalVolume* m_Electronic;
+ G4LogicalVolume* m_Stars;
+ G4LogicalVolume* m_Chamber;
+ G4LogicalVolume* m_Blades;
+ G4LogicalVolume* m_Base;
+
+ string ChamberPath;
+ string BasePath;
+ string StarsPath;
+ string BladesPath;
+ bool found_chamber;
+ bool found_blades;
+ bool found_stars;
+ bool found_base;
+
+ private:
+ // Initialize material used in detector definition
+ void InitializeMaterial();
+
+ // List of material
+ G4Material* m_MaterialSilicon;
+ G4Material* m_MaterialAl;
+ G4Material* m_MaterialVacuum;
+ G4Material* m_MaterialPCB;
+ G4Material* m_MaterialCu;
+
+ // calculated dimension
+ double m_Active_InnerWafer_Width;
+ double m_Active_InnerWafer_Length;
+ double m_Active_OuterWafer_Width;
+ double m_Active_OuterWafer_Length;
+
+ ////////////////////////////////////////////////////
+ ////// Inherite from NPS::VDetector class /////////
+ ////////////////////////////////////////////////////
+ public:
+ // Read stream at Configfile to pick-up parameters of detector (Position,...)
+ // Called in DetecorConstruction::ReadDetextorConfiguration Method
+ void ReadConfiguration(NPL::InputParser);
+ // Construct detector and inialise sensitive part.
+ // Called After DetecorConstruction::AddDetector Method
+ void ConstructDetector(G4LogicalVolume* world);
+
+ // Add Detector branch to the EventTree.
+ // Called After DetecorConstruction::AddDetector Method
+ void InitializeRootOutput();
+
+ // Read sensitive part and fill the Root tree.
+ // Called at in the EventAction::EndOfEventAvtion
+ void ReadSensitive(const G4Event* event);
+
+ public: // Scorer
+ // Initialize all Scorer used by the MUST2Array
+ void InitializeScorers();
+
+ // Associated Scorer
+ G4MultiFunctionalDetector* m_InnerScorer1;
+ G4MultiFunctionalDetector* m_OuterScorer1;
+ G4MultiFunctionalDetector* m_InnerScorer2;
+ G4MultiFunctionalDetector* m_OuterScorer2;
+
+ ////////////////////////////////////////////////////
+ ///////////Event class to store Data////////////////
+ ////////////////////////////////////////////////////
+ private:
+ TStrasseData* m_Event;
+
+ ////////////////////////////////////////////////////
+ ///////////////Private intern Data//////////////////
+ ////////////////////////////////////////////////////
+ private: // Geometry
+ // Detector Coordinate
+ vector<double> m_Inner_R;
+ vector<double> m_Inner_Z;
+ vector<double> m_Inner_Phi;
+ vector<double> m_Inner_Shift;
+ vector<G4ThreeVector> m_Inner_Ref;
+
+ vector<double> m_Outer_R;
+ vector<double> m_Outer_Z;
+ vector<double> m_Outer_Phi;
+ vector<double> m_Outer_Shift;
+ vector<G4ThreeVector> m_Outer_Ref;
+
+ vector<double> m_Target_R;
+ vector<double> m_Target_L;
+ vector<string> m_Target_MaterialName;
+ vector<string> m_Target_CellMaterialName;
+ vector<double> m_Target_CellThickness;
+ vector<G4ThreeVector> m_Target_Pos;
+
+ vector<double> m_Chamber_Z;
+
+ // Region were reaction can occure:
+ G4Region* m_ReactionRegion;
+ vector<G4VFastSimulationModel*> m_ReactionModel;
+
+ // Needed for dynamic loading of the library
+ public:
+ static NPS::VDetector* Construct();
+
+ private: // Visualisation
+ G4VisAttributes* SiliconVisAtt;
+ G4VisAttributes* PCBVisAtt;
+ G4VisAttributes* PADVisAtt;
+ G4VisAttributes* StarsVisAtt;
+ G4VisAttributes* ChamberVisAtt;
+ G4VisAttributes* GuardRingVisAtt;
+ G4VisAttributes* BladeVisAtt;
+ G4VisAttributes* TargetVisAtt;
+ G4VisAttributes* TargetCellVisAtt;
};
#endif
diff --git a/NPSimulation/Process/BeamReaction.cc b/NPSimulation/Process/BeamReaction.cc
index 73ae9ab1de7432b17f208c9262322fa38b4b8c04..d1fd91c7c81adb7f00a804ef6356ce82dc078dab 100644
--- a/NPSimulation/Process/BeamReaction.cc
+++ b/NPSimulation/Process/BeamReaction.cc
@@ -500,14 +500,14 @@ void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack, G4FastStep& fastStep)
// Momentum in beam and world frame for light particle 1
G4ThreeVector momentum_kine1_beam(sin(Theta1) * cos(Phi1), sin(Theta1) * sin(Phi1), cos(Theta1));
G4ThreeVector momentum_kine1_world = momentum_kine1_beam;
- momentum_kine1_world.rotate(Beam_theta, uy); // rotation of Beam_theta around Y axis
- momentum_kine1_world.rotate(Beam_phi, uz); // rotation of Beam_phi around Z axis
+ // momentum_kine1_world.rotate(Beam_theta, uy); // rotation of Beam_theta around Y axis
+ // momentum_kine1_world.rotate(Beam_phi, uz); // rotation of Beam_phi around Z axis
// Momentum in beam and world frame for light particle 2
G4ThreeVector momentum_kine2_beam(sin(Theta2) * cos(Phi2), sin(Theta2) * sin(Phi2), cos(Theta2));
G4ThreeVector momentum_kine2_world = momentum_kine2_beam;
- momentum_kine2_world.rotate(Beam_theta, uy); // rotation of Beam_theta on Y axis
- momentum_kine2_world.rotate(Beam_phi, uz); // rotation of Beam_phi on Z axis
+ // momentum_kine2_world.rotate(Beam_theta, uy); // rotation of Beam_theta on Y axis
+ // momentum_kine2_world.rotate(Beam_phi, uz); // rotation of Beam_phi on Z axis
// Momentum in beam and world frame for heavy residual
//
diff --git a/Projects/Strasse/Analysis.cxx b/Projects/Strasse/Analysis.cxx
index 61182b5f30ee430848c6c40732b61c3d7db759f8..36991d64d018cd427d97d62ad5b4a48b68318e90 100644
--- a/Projects/Strasse/Analysis.cxx
+++ b/Projects/Strasse/Analysis.cxx
@@ -19,69 +19,69 @@
* *
* *
*****************************************************************************/
-#include<iostream>
+#include <iostream>
using namespace std;
-#include"Analysis.h"
-#include"NPAnalysisFactory.h"
-#include"NPDetectorManager.h"
-#include"NPOptionManager.h"
-#include"NPFunction.h"
-#include"NPTrackingUtility.h"
-#include"NPPhysicalConstants.h"
-#include"TRandom3.h"
+#include "Analysis.h"
+#include "NPAnalysisFactory.h"
+#include "NPDetectorManager.h"
+#include "NPFunction.h"
+#include "NPOptionManager.h"
+#include "NPPhysicalConstants.h"
+#include "NPTrackingUtility.h"
+#include "TRandom3.h"
////////////////////////////////////////////////////////////////////////////////
-Analysis::Analysis(){
-}
+Analysis::Analysis() {}
////////////////////////////////////////////////////////////////////////////////
-Analysis::~Analysis(){
-}
+Analysis::~Analysis() {}
////////////////////////////////////////////////////////////////////////////////
-void Analysis::Init(){
- IC= new TInitialConditions;
- DC= new TInteractionCoordinates;
- RC= new TReactionConditions;
-
+void Analysis::Init() {
+ IC = new TInitialConditions;
+ DC = new TInteractionCoordinates;
+ RC = new TReactionConditions;
+
InitOutputBranch();
InitInputBranch();
-
- Strasse = (TStrassePhysics*) m_DetectorManager -> GetDetector("Strasse");
-// Catana = (TCatanaPhysics*) m_DetectorManager -> GetDetector("Catana");
+
+ Strasse = (TStrassePhysics*)m_DetectorManager->GetDetector("Strasse");
+ // Catana = (TCatanaPhysics*) m_DetectorManager -> GetDetector("Catana");
// reaction properties
m_QFS = new NPL::QFS();
m_QFS->ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
// reaction properties
myBeam = new NPL::Beam();
myBeam->ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
- InitialBeamEnergy = myBeam->GetEnergy()/* * myBeam->GetA()*/;
+ InitialBeamEnergy = myBeam->GetEnergy() /* * myBeam->GetA()*/;
// target thickness
- TargetThickness = m_DetectorManager->GetTargetThickness();
- string TargetMaterial = m_DetectorManager->GetTargetMaterial();
+ // TargetThickness = m_DetectorManager->GetTargetThickness();
+ // string TargetMaterial = m_DetectorManager->GetTargetMaterial();
+ TargetThickness = 150 * mm;
+ string TargetMaterial = "LH2";
// EnergyLoss Tables
string BeamName = NPL::ChangeNameToG4Standard(myBeam->GetName());
- BeamTarget = NPL::EnergyLoss(BeamName+"_"+TargetMaterial+".G4table","G4Table",100);
- protonTarget = NPL::EnergyLoss("proton_"+TargetMaterial+".G4table","G4Table",100);
- protonAl = NPL::EnergyLoss("proton_Al.G4table","G4Table",100);
- protonSi = NPL::EnergyLoss("proton_Si.G4table","G4Table",100);
- LV_T.SetVectM(TVector3(0,0,0),NPUNITS::proton_mass_c2);
+ BeamTarget = NPL::EnergyLoss(BeamName + "_" + TargetMaterial + ".G4table", "G4Table", 100);
+ protonTarget = NPL::EnergyLoss("proton_" + TargetMaterial + ".G4table", "G4Table", 100);
+ protonAl = NPL::EnergyLoss("proton_Al.G4table", "G4Table", 100);
+ protonSi = NPL::EnergyLoss("proton_Si.G4table", "G4Table", 100);
+ 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
ReInitValue();
- bool perfectBeamEnergy=false;
- bool perfectBeamTracking=false;
- bool perfectProtonTracking=false;
- bool perfectProtonEnergy=false;
- bool CATANAEnergyReco=false;
+ bool perfectBeamEnergy = false;
+ bool perfectBeamTracking = false;
+ bool perfectProtonTracking = false;
+ bool perfectProtonEnergy = false;
+ bool CATANAEnergyReco = false;
unsigned int size = Strasse->GetEventMultiplicity();
- if(size==2){ // 2 proton detected
+ if (size == 2) { // 2 proton detected
////////////////////////////////////////////////
////////////Protons (E,P) Reconstruction ///////
@@ -89,86 +89,99 @@ void Analysis::TreatEvent(){
// Proton 1
TVector3 InnerPos1 = Strasse->GetInnerPositionOfInteraction(0);
TVector3 OuterPos1 = Strasse->GetOuterPositionOfInteraction(0);
- TVector3 Proton1 = OuterPos1-InnerPos1;
- Proton1=Proton1.Unit();
+ 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();
+ TVector3 Proton2 = OuterPos2 - InnerPos2;
+ Proton2 = Proton2.Unit();
- deltaPhi = abs(Proton1.Phi()/deg-Proton2.Phi()/deg);
- sumTheta = Proton1.Theta()/deg+Proton2.Theta()/deg;
- Theta12 = Proton1.Angle(Proton2)/deg;
+ deltaPhi = abs(Proton1.Phi() / deg - Proton2.Phi() / deg);
+ sumTheta = Proton1.Theta() / deg + Proton2.Theta() / deg;
+ Theta12 = Proton1.Angle(Proton2) / deg;
// computing minimum distance of the two lines
TVector3 Vertex;
TVector3 delta;
- Distance = NPL::MinimumDistanceTwoLines(InnerPos1,OuterPos1,InnerPos2,OuterPos2,Vertex,delta);
- VertexX=Vertex.X();
- VertexY=Vertex.Y();
- VertexZ=Vertex.Z();
- deltaX=delta.X();
- deltaY=delta.Y();
- deltaZ=delta.Z();
-
- if(perfectProtonTracking){
- // From Reaction Conditions
- VertexX=RC->GetVertexPositionX();
- VertexY=RC->GetVertexPositionY();
- VertexZ=RC->GetVertexPositionZ();
- Proton2=RC->GetParticleMomentum(0).Unit();
- Proton1=RC->GetParticleMomentum(1).Unit();
- deltaX=0;
- deltaY=0;
- deltaZ=0;
+ Distance = NPL::MinimumDistanceTwoLines(InnerPos1, OuterPos1, InnerPos2, OuterPos2, Vertex, delta);
+ VertexX = Vertex.X();
+ VertexY = Vertex.Y();
+ VertexZ = Vertex.Z();
+ deltaX = delta.X();
+ deltaY = delta.Y();
+ deltaZ = delta.Z();
+
+ if (perfectProtonTracking) {
+ // From Reaction Conditions
+ VertexX = RC->GetVertexPositionX();
+ VertexY = RC->GetVertexPositionY();
+ VertexZ = RC->GetVertexPositionZ();
+ Proton2 = RC->GetParticleDirection(0).Unit();
+ Proton1 = RC->GetParticleDirection(1).Unit();
+ deltaX = 0;
+ deltaY = 0;
+ deltaZ = 0;
}
- if(perfectProtonEnergy){
- E1 = RC->GetKineticEnergy(1);
- E2 = RC->GetKineticEnergy(0);
+ if (perfectProtonEnergy) {
+ E1 = RC->GetKineticEnergy(1);
+ E2 = RC->GetKineticEnergy(0);
}
else {
- //General CATANA resolution for proton (1.4 % FWHM)
- // simply applied to E, crude approximation (no CATANA reco.)
- E1 = ApplyCATANAResoProton(RC->GetKineticEnergy(1));
- E2 = ApplyCATANAResoProton(RC->GetKineticEnergy(0));
+ // General CATANA resolution for proton (1.4 % FWHM)
+ // simply applied to E, crude approximation (no CATANA reco.)
+ E1 = ApplyCATANAResoProton(RC->GetKineticEnergy(1));
+ E2 = ApplyCATANAResoProton(RC->GetKineticEnergy(0));
}
- if(CATANAEnergyReco){
- //if true, bypass previous proton energy calc. and use full CATANA Reco.
- ///////////
- // TO DO //
- ///////////
+ if (CATANAEnergyReco) {
+ // if true, bypass previous proton energy calc. and use full CATANA Reco.
+ ///////////
+ // TO DO //
+ ///////////
}
- double P1= sqrt(E1*(E1+2*NPUNITS::proton_mass_c2));
- double P2= sqrt(E2*(E2+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));
///////////////////////////////////////
//////Beam (E,P) Reconstruction ///////
///////////////////////////////////////
TVector3 BeamDirection;
- if(perfectBeamTracking) BeamDirection = RC->GetBeamDirection();
- else BeamDirection = TVector3(0,0,1);
- //Beam Energy at vertex
+ if (perfectBeamTracking)
+ BeamDirection = RC->GetBeamDirection();
+ else
+ BeamDirection = TVector3(0, 0, 1);
+
+ // Beam Energy at vertex
TA = RC->GetBeamEnergy();
- if(perfectBeamEnergy) TAcalc = TA;
- else 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();
+ if (perfectBeamEnergy) {
+ TAcalc = TA;
+ }
+ 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 //////
//////////////////////////////////////////////
- LV_A.SetVectM(PA,m_QFS->GetParticleA()->Mass());
- LV_p1.SetVectM(Proton1.Unit()*P1,NPUNITS::proton_mass_c2);
- LV_p2.SetVectM(Proton2.Unit()*P2,NPUNITS::proton_mass_c2);
+ LV_A.SetVectM(PA, m_QFS->GetParticleA()->Mass());
+ 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);
+ // LV_B = RC->GetParticleMomentum(2);
Ex = LV_B.M() - m_QFS->GetParticleB()->Mass();
- }//endif size==2
+ } // endif size==2
}
/* ///////////////////////////////////
@@ -187,15 +200,15 @@ cout<<"Catana mult="<<Catana->GetEventMultiplicity()<<endl;
//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);
+ 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;
@@ -208,7 +221,7 @@ cout<<"Catana mult="<<Catana->GetEventMultiplicity()<<endl;
Phi1s=Proton1s.Phi();
Theta2s=Proton2s.Theta();
Phi2s=Proton2s.Phi();
-
+
}
else{
E2=E1s;
@@ -225,169 +238,161 @@ cout<<"Catana mult="<<Catana->GetEventMultiplicity()<<endl;
}
// From detectors
- E1 = ReconstructProtonEnergy(Vertex,Proton1,Catana->Energy[i1]);
+ E1 = ReconstructProtonEnergy(Vertex,Proton1,Catana->Energy[i1]);
E2 = ReconstructProtonEnergy(Vertex,Proton2,Catana->Energy[i2]);
- //E1 = Catana->Energy[i1];
+ //E1 = Catana->Energy[i1];
//E2 = Catana->Energy[i2];
cout<<"------------------------"<<endl;
*/
////////////////////////////////////////////////////////////////////////////////
-double Analysis::ReconstructProtonEnergy(const TVector3& x0, const TVector3& dir,const double& Ecatana){
-
- TVector3 Normal = TVector3(0,1,0);
- Normal.SetPhi(dir.Phi());
- double Theta = dir.Angle(Normal);
- // Catana Al housing
- double E = protonAl.EvaluateInitialEnergy(Ecatana,0.5*mm,Theta);
- cout<<"Ecatana="<<Ecatana<<endl;
- cout<<"Erec0="<<E<<endl;
- // Strasse Chamber
- //E = protonAl.EvaluateInitialEnergy(E,3*mm,Theta);
- // Outer Barrel
- E = protonSi.EvaluateInitialEnergy(E,300*micrometer,Theta);
- // Inner Barrel
- E = protonSi.EvaluateInitialEnergy(E,200*micrometer,Theta);
- // LH2 target
- static TVector3 x1;
- x1= x0+dir;
- TVector3 T1(0,15,0);
- TVector3 T2(0,15,1);
- T1.SetPhi(dir.Phi());
- T2.SetPhi(dir.Phi());
- // double d = NPL::MinimumDistancePointLine(T1,T2,x0);
- double d = 0;
-
- cout<<"d="<<d<<endl;
- cout<<"Theta="<<Theta<<endl;
- E = protonTarget.EvaluateInitialEnergy(E,d,Theta);
- cout<<"Erec="<<E<<endl;
- return E;
-
- //return 0;
+double Analysis::ReconstructProtonEnergy(const TVector3& x0, const TVector3& dir, const double& Ecatana) {
+
+ TVector3 Normal = TVector3(0, 1, 0);
+ Normal.SetPhi(dir.Phi());
+ double Theta = dir.Angle(Normal);
+ // Catana Al housing
+ double E = protonAl.EvaluateInitialEnergy(Ecatana, 0.5 * mm, Theta);
+ cout << "Ecatana=" << Ecatana << endl;
+ cout << "Erec0=" << E << endl;
+ // Strasse Chamber
+ // E = protonAl.EvaluateInitialEnergy(E,3*mm,Theta);
+ // Outer Barrel
+ E = protonSi.EvaluateInitialEnergy(E, 300 * micrometer, Theta);
+ // Inner Barrel
+ E = protonSi.EvaluateInitialEnergy(E, 200 * micrometer, Theta);
+ // LH2 target
+ static TVector3 x1;
+ x1 = x0 + dir;
+ TVector3 T1(0, 15, 0);
+ TVector3 T2(0, 15, 1);
+ T1.SetPhi(dir.Phi());
+ T2.SetPhi(dir.Phi());
+ // double d = NPL::MinimumDistancePointLine(T1,T2,x0);
+ double d = 0;
+
+ cout << "d=" << d << endl;
+ cout << "Theta=" << Theta << endl;
+ E = protonTarget.EvaluateInitialEnergy(E, d, Theta);
+ cout << "Erec=" << E << endl;
+ return E;
+
+ // return 0;
}
-
////////////////////////////////////////////////////////////////////////////////
-double Analysis::ApplyCATANAResoProton(double Ein){
- // Function applying overall CATANA crystal resolution
- // For protons (1.4% FWHM)
- double FWHM = 1.4/100 * Ein;
- double sigma = FWHM/2.35;
- double Eout = rand->Gaus(Ein,sigma);
- return Eout;
+double Analysis::ApplyCATANAResoProton(double Ein) {
+ // Function applying overall CATANA crystal resolution
+ // For protons (1.4% FWHM)
+ double FWHM = 1.4 / 100 * Ein;
+ double sigma = FWHM / 2.35;
+ double Eout = rand->Gaus(Ein, sigma);
+ return Eout;
}
////////////////////////////////////////////////////////////////////////////////
-double Analysis::ApplyCATANAResoGamma(double Ein){
- // Function applying overall CATANA crystal resolution
- // For gammas defined in smsimulator package
- double a = 0.686569;
- double b = 0.564352;
- // Ein from MeV to keV
- Ein = Ein * 1000;
- double SigmaE = a * pow(Ein,b);
- double Eout = rand->Gaus(Ein,SigmaE);
- return Eout/1000.;
+double Analysis::ApplyCATANAResoGamma(double Ein) {
+ // Function applying overall CATANA crystal resolution
+ // For gammas defined in smsimulator package
+ double a = 0.686569;
+ double b = 0.564352;
+ // Ein from MeV to keV
+ Ein = Ein * 1000;
+ double SigmaE = a * pow(Ein, b);
+ double Eout = rand->Gaus(Ein, SigmaE);
+ return Eout / 1000.;
}
////////////////////////////////////////////////////////////////////////////////
-TVector3 Analysis::InterpolateInPlaneZ(TVector3 V0, TVector3 V1, double Zproj){
- TVector3 Vproj(-999,-999,-999);
- double t = (Zproj - V1.Z()) / (V1.Z()-V0.Z());
- double Xproj= V1.X() + (V1.X()-V0.X()) * t;
- double Yproj= V1.Y() + (V1.Y()-V0.Y()) * t;
- Vproj.SetXYZ(Xproj,Yproj,Zproj);
- return Vproj;
+TVector3 Analysis::InterpolateInPlaneZ(TVector3 V0, TVector3 V1, double Zproj) {
+ TVector3 Vproj(-999, -999, -999);
+ double t = (Zproj - V1.Z()) / (V1.Z() - V0.Z());
+ double Xproj = V1.X() + (V1.X() - V0.X()) * t;
+ double Yproj = V1.Y() + (V1.Y() - V0.Y()) * t;
+ Vproj.SetXYZ(Xproj, Yproj, Zproj);
+ return Vproj;
}
////////////////////////////////////////////////////////////////////////////////
-void Analysis::End(){
-}
+void Analysis::End() {}
////////////////////////////////////////////////////////////////////////////////
void Analysis::InitOutputBranch() {
- RootOutput::getInstance()->GetTree()->Branch("Ex",&Ex,"Ex/D");
- RootOutput::getInstance()->GetTree()->Branch("E1",&E1,"E1/D");
- RootOutput::getInstance()->GetTree()->Branch("E2",&E2,"E2/D");
- RootOutput::getInstance()->GetTree()->Branch("Theta12",&Theta12,"Theta12/D");
- RootOutput::getInstance()->GetTree()->Branch("ThetaCM",&ThetaCM,"ThetaCM/D");
- RootOutput::getInstance()->GetTree()->Branch("VertexX",&VertexX,"VertexX/D");
- RootOutput::getInstance()->GetTree()->Branch("VertexY",&VertexY,"VertexY/D");
- RootOutput::getInstance()->GetTree()->Branch("VertexZ",&VertexZ,"VertexZ/D");
- RootOutput::getInstance()->GetTree()->Branch("deltaX",&deltaX,"deltaX/D");
- RootOutput::getInstance()->GetTree()->Branch("deltaY",&deltaY,"deltaY/D");
- RootOutput::getInstance()->GetTree()->Branch("deltaZ",&deltaZ,"deltaZ/D");
- 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("TA",&TA,"TA/D");
- RootOutput::getInstance()->GetTree()->Branch("TAcalc",&TAcalc,"TAcalc/D");
-
-
- RootOutput::getInstance()->GetTree()->Branch("Distance",&Distance,"Distance/D");
- RootOutput::getInstance()->GetTree()->Branch("InteractionCoordinates","TInteractionCoordinates",&DC);
- RootOutput::getInstance()->GetTree()->Branch("ReactionConditions","TReactionConditions",&RC);
+ RootOutput::getInstance()->GetTree()->Branch("Ex", &Ex, "Ex/D");
+ RootOutput::getInstance()->GetTree()->Branch("E1", &E1, "E1/D");
+ RootOutput::getInstance()->GetTree()->Branch("E2", &E2, "E2/D");
+ RootOutput::getInstance()->GetTree()->Branch("Theta12", &Theta12, "Theta12/D");
+ RootOutput::getInstance()->GetTree()->Branch("ThetaCM", &ThetaCM, "ThetaCM/D");
+ RootOutput::getInstance()->GetTree()->Branch("VertexX", &VertexX, "VertexX/D");
+ RootOutput::getInstance()->GetTree()->Branch("VertexY", &VertexY, "VertexY/D");
+ RootOutput::getInstance()->GetTree()->Branch("VertexZ", &VertexZ, "VertexZ/D");
+ RootOutput::getInstance()->GetTree()->Branch("deltaX", &deltaX, "deltaX/D");
+ RootOutput::getInstance()->GetTree()->Branch("deltaY", &deltaY, "deltaY/D");
+ RootOutput::getInstance()->GetTree()->Branch("deltaZ", &deltaZ, "deltaZ/D");
+ 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("TA", &TA, "TA/D");
+ RootOutput::getInstance()->GetTree()->Branch("TAcalc", &TAcalc, "TAcalc/D");
+
+ RootOutput::getInstance()->GetTree()->Branch("Distance", &Distance, "Distance/D");
+ RootOutput::getInstance()->GetTree()->Branch("InteractionCoordinates", "TInteractionCoordinates", &DC);
+ RootOutput::getInstance()->GetTree()->Branch("ReactionConditions", "TReactionConditions", &RC);
}
////////////////////////////////////////////////////////////////////////////////
-void Analysis::InitInputBranch(){
- RootInput:: getInstance()->GetChain()->SetBranchAddress("InteractionCoordinates",&DC);
- RootInput:: getInstance()->GetChain()->SetBranchAddress("ReactionConditions",&RC);
+void Analysis::InitInputBranch() {
+ RootInput::getInstance()->GetChain()->SetBranchAddress("InteractionCoordinates", &DC);
+ RootInput::getInstance()->GetChain()->SetBranchAddress("ReactionConditions", &RC);
}
////////////////////////////////////////////////////////////////////////////////
-void Analysis::ReInitValue(){
- Ex = -1000 ;
- E1= -1000;
+void Analysis::ReInitValue() {
+ Ex = -1000;
+ E1 = -1000;
E2 = -1000;
Theta12 = -1000;
ThetaCM = -1000;
- VertexX=-1000;
- VertexY=-1000;
- VertexZ=-1000;
- deltaX=-1000;
- deltaY=-1000;
- deltaZ=-1000;
- Distance=-1000;
- sumTheta=-1000;
- deltaPhi=-1000;
- alpha=-1000;
- Theta1=-1000;
- Phi1=-1000;
- Theta2=-1000;
- Phi2=-1000;
- Theta1s=-1000;
- Phi1s=-1000;
- Theta2s=-1000;
- Phi2s=-1000;
- TA=-1000;
- TAcalc=-1000;
+ VertexX = -1000;
+ VertexY = -1000;
+ VertexZ = -1000;
+ deltaX = -1000;
+ deltaY = -1000;
+ deltaZ = -1000;
+ Distance = -1000;
+ sumTheta = -1000;
+ deltaPhi = -1000;
+ alpha = -1000;
+ Theta1 = -1000;
+ Phi1 = -1000;
+ Theta2 = -1000;
+ Phi2 = -1000;
+ Theta1s = -1000;
+ Phi1s = -1000;
+ Theta2s = -1000;
+ Phi2s = -1000;
+ TA = -1000;
+ TAcalc = -1000;
}
-
////////////////////////////////////////////////////////////////////////////////
// Construct Method to be pass to the AnalysisFactory //
////////////////////////////////////////////////////////////////////////////////
-NPL::VAnalysis* Analysis::Construct(){
- return (NPL::VAnalysis*) new Analysis();
-}
+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);
- }
+extern "C" {
+class proxy {
+ public:
+ proxy() { NPL::AnalysisFactory::getInstance()->SetConstructor(Analysis::Construct); }
};
proxy p;
diff --git a/Projects/Strasse/PhysicsListOption.txt b/Projects/Strasse/PhysicsListOption.txt
index 495166f434c7f2377d07f89ba844aa37411514a3..56772eef01e5db197b2641b663872b7d1506c70f 100644
--- a/Projects/Strasse/PhysicsListOption.txt
+++ b/Projects/Strasse/PhysicsListOption.txt
@@ -1,5 +1,5 @@
EmPhysicsList Option4
-DefaultCutOff 1000
+DefaultCutOff 1000000
IonBinaryCascadePhysics 0
NPIonInelasticPhysics 0
EmExtraPhysics 0
@@ -8,4 +8,4 @@ StoppingPhysics 0
OpticalPhysics 0
HadronPhysicsINCLXX 0
HadronPhysicsQGSP_BIC_HP 0
-Decay 1
+Decay 0
diff --git a/Projects/Strasse/geometry/strasse_June2022.detector b/Projects/Strasse/geometry/strasse_June2022.detector
new file mode 100644
index 0000000000000000000000000000000000000000..2da7c74f1bd3be5d0f1d1c75728c57499cb29ce9
--- /dev/null
+++ b/Projects/Strasse/geometry/strasse_June2022.detector
@@ -0,0 +1,101 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Target
+ THICKNESS= 150 mm
+ ANGLE= 0 deg
+ RADIUS= 15 mm
+ MATERIAL= LH2
+ X= 0 mm
+ Y= 0 mm
+ Z= 0 mm
+ NbSlices= 150
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1
+Strasse Info
+ Inner_Wafer_Length= 124 mm
+ Inner_Wafer_Width= 32 mm
+ Inner_Wafer_Thickness= 200 micrometer
+ Inner_Wafer_AlThickness= 0.4 micrometer
+ Inner_Wafer_PADExternal= 0 mm
+ Inner_Wafer_PADInternal= 0 mm
+ Inner_Wafer_GuardRing= 1.0 mm
+ Inner_Wafer_TransverseStrips= 610
+ Inner_Wafer_LongitudinalStrips= 150
+ Inner_PCB_PortWidth= 2.0 mm
+ Inner_PCB_StarboardWidth= 2.0 mm
+ Inner_PCB_BevelAngle= 90 deg
+ Inner_PCB_UpstreamWidth= 12 mm
+ Inner_PCB_DownstreamWidth= 39 mm
+ Inner_PCB_MidWidth= 0.001 mm
+ Inner_PCB_Thickness= 2.4 mm
+ Inner_PCB_Ledge= 1.0 mm
+ Inner_PCB_Step= 2.1 mm
+ Outer_Wafer_Length= 123 mm
+ Outer_Wafer_Width= 64.6 mm
+ Outer_Wafer_Thickness= 300 micrometer
+ Outer_Wafer_AlThickness= 0.4 micrometer
+ Outer_Wafer_PADExternal= 0 mm
+ Outer_Wafer_PADInternal= 0 mm
+ Outer_Wafer_GuardRing= 1.0 mm
+ Outer_PCB_PortWidth= 2.0 mm
+ Outer_PCB_StarboardWidth= 2.0 mm
+ Outer_PCB_BevelAngle= 90 deg
+ Outer_PCB_UpstreamWidth= 40 mm
+ Outer_PCB_DownstreamWidth= 13 mm
+ Outer_PCB_MidWidth= 0.001 mm
+ Outer_PCB_Thickness= 2.4 mm
+ Outer_PCB_Ledge= 1.0 mm
+ Outer_PCB_Step= 2.0 mm
+ Outer_Wafer_TransverseStrips= 605
+ Outer_Wafer_LongitudinalStrips= 313
+ % unused if using CAD file (.stl) chamber
+ Chamber_Thickness= 3 mm
+ Chamber_Cylinder_Length= 360 mm
+ Chamber_Radius= 180 mm
+ Chamber_ExitTube_Radius= 79.5 mm
+ Chamber_ExitTube_Length= 100 mm
+ Chamber_Flange_Inner_Radius= 50 mm
+ Chamber_Sphere_Radius= 220 mm
+ Chamber_Sphere_Shift= 60 mm
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Alias InnerPhi
+ Action= Copy
+ %Value= 0
+ Value= -5.5 54.5 114.5 174.5 234.5 294.5
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Strasse Inner
+ Radius= 28.7 mm
+ Z= 76.5 mm
+ Phi= @InnerPhi deg
+ Shift= 3.5 mm
+ Ref= 0 0 0 mm
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Alias OuterPhi
+ Action= Copy
+ %Value= 0
+ Value= 0 60 120 180 240 300
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Strasse Outer
+ Radius= 60.7 mm
+ Z= 76.5 mm
+ Phi= @OuterPhi deg
+ Shift= 0 mm
+ Ref= 0 0 0 mm
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%Strasse InactiveMaterial
+% Chamber= ./geometry/STRASSE_Chamber.stl
+% Stars= ./geometry/STRASSE_StarSupports.stl
+% Base= ./geometry/STRASSE_Base.stl
+% Blades= ./geometry/STRASSE_Blades.stl
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1
+%Catana CSV
+% Path= geometry/Catana.csv
+% Pos= 0 0 100 mm
+% Rshift= 100 micrometer
+
+
diff --git a/Projects/Strasse/geometry/strasse_last.detector b/Projects/Strasse/geometry/strasse_last.detector
new file mode 100644
index 0000000000000000000000000000000000000000..2ad1541154f47500e6c15296d513b23731092bcb
--- /dev/null
+++ b/Projects/Strasse/geometry/strasse_last.detector
@@ -0,0 +1,96 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1
+Strasse Info
+ Inner_Wafer_Length= 124 mm
+ Inner_Wafer_Width= 32 mm
+ Inner_Wafer_Thickness= 200 micrometer
+ Inner_Wafer_AlThickness= 0.4 micrometer
+ Inner_Wafer_PADExternal= 0 mm
+ Inner_Wafer_PADInternal= 0 mm
+ Inner_Wafer_GuardRing= 1.0 mm
+ Inner_Wafer_TransverseStrips= 610
+ Inner_Wafer_LongitudinalStrips= 150
+ Inner_PCB_PortWidth= 2.0 mm
+ Inner_PCB_StarboardWidth= 2.0 mm
+ Inner_PCB_BevelAngle= 90 deg
+ Inner_PCB_UpstreamWidth= 12 mm
+ Inner_PCB_DownstreamWidth= 39 mm
+ Inner_PCB_MidWidth= 0.001 mm
+ Inner_PCB_Thickness= 2.4 mm
+ Inner_PCB_Ledge= 1.0 mm
+ Inner_PCB_Step= 2.1 mm
+ Outer_Wafer_Length= 123 mm
+ Outer_Wafer_Width= 64.6 mm
+ Outer_Wafer_Thickness= 300 micrometer
+ Outer_Wafer_AlThickness= 0.4 micrometer
+ Outer_Wafer_PADExternal= 0 mm
+ Outer_Wafer_PADInternal= 0 mm
+ Outer_Wafer_GuardRing= 1.0 mm
+ Outer_PCB_PortWidth= 2.0 mm
+ Outer_PCB_StarboardWidth= 2.0 mm
+ Outer_PCB_BevelAngle= 90 deg
+ Outer_PCB_UpstreamWidth= 40 mm
+ Outer_PCB_DownstreamWidth= 13 mm
+ Outer_PCB_MidWidth= 0.001 mm
+ Outer_PCB_Thickness= 2.4 mm
+ Outer_PCB_Ledge= 1.0 mm
+ Outer_PCB_Step= 2.0 mm
+ Outer_Wafer_TransverseStrips= 605
+ Outer_Wafer_LongitudinalStrips= 313
+ % unused if using CAD file (.stl) chamber
+ Chamber_Thickness= 3 mm
+ Chamber_Cylinder_Length= 360 mm
+ Chamber_Radius= 180 mm
+ Chamber_ExitTube_Radius= 79.5 mm
+ Chamber_ExitTube_Length= 100 mm
+ Chamber_Flange_Inner_Radius= 50 mm
+ Chamber_Sphere_Radius= 220 mm
+ Chamber_Sphere_Shift= 60 mm
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Strasse Target
+ Radius= 15 mm
+ Length= 150 mm
+ TargetMaterial= LH2
+ CellMaterial= Mylar
+ CellThickness= 0.15 mm
+ Pos= 0 0 0 mm
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Alias InnerPhi
+ Action= Copy
+ %Value= 0
+ Value= -5.5 54.5 114.5 174.5 234.5 294.5
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Strasse Inner
+ Radius= 28.7 mm
+ Z= 76.5 mm
+ Phi= @InnerPhi deg
+ Shift= 3.5 mm
+ Ref= 0 0 0 mm
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Alias OuterPhi
+ Action= Copy
+ %Value= 0
+ Value= 0 60 120 180 240 300
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Strasse Outer
+ Radius= 60.7 mm
+ Z= 76.5 mm
+ Phi= @OuterPhi deg
+ Shift= 0 mm
+ Ref= 0 0 0 mm
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%Strasse InactiveMaterial
+% Chamber= ./geometry/STRASSE_Chamber.stl
+% Stars= ./geometry/STRASSE_StarSupports.stl
+% Base= ./geometry/STRASSE_Base.stl
+% Blades= ./geometry/STRASSE_Blades.stl
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%1
+%Catana CSV
+% Path= geometry/Catana.csv
+% Pos= 0 0 100 mm
+% Rshift= 100 micrometer
+
+
diff --git a/Projects/Strasse/macro/CalculateDetectorOffset.cxx b/Projects/Strasse/macro/CalculateDetectorOffset.cxx
index bef189e72f89ebf42768c79884c3ab4e9ca01e19..68d056a933c7570b555541b737ce5191af82a36f 100644
--- a/Projects/Strasse/macro/CalculateDetectorOffset.cxx
+++ b/Projects/Strasse/macro/CalculateDetectorOffset.cxx
@@ -1,6 +1,6 @@
-#include <iostream>
-#include <ctime>
#include <cstdlib>
+#include <ctime>
+#include <iostream>
using namespace std;
// ROOT headers
@@ -11,245 +11,220 @@ using namespace std;
using namespace NPL;
-
-void CalculateDetectorOffset(const char * fname = "./geometry/strasse_July2021"){
+// void CalculateDetectorOffset(const char * fname = "./geometry/strasse_July2021"){
+void CalculateDetectorOffset(const char* fname = "./geometry/strasse_June2022") {
// Open output ROOT file from NPTool simulation run
string path = "";
string inFileName = fname;
inFileName += ".detector";
+ InputParser* inParser = new InputParser(inFileName, true);
+ vector<NPL::InputBlock*> blocks_info = inParser->GetAllBlocksWithTokenAndValue("Strasse", "Info");
- InputParser *inParser = new InputParser(inFileName,true);
- vector<NPL::InputBlock*> blocks_info = inParser->GetAllBlocksWithTokenAndValue("Strasse","Info");
-
- if(blocks_info.size()>1){
- cout << "ERROR: can only accepte one info block, " << blocks_info.size() << " info block founds." << endl;
- exit(1);
+ if (blocks_info.size() > 1) {
+ cout << "ERROR: can only accepte one info block, " << blocks_info.size() << " info block founds." << endl;
+ exit(1);
}
-
- vector<string> info = {
- "Inner_Wafer_Length",
- "Inner_Wafer_Width",
- "Inner_Wafer_Thickness",
- "Inner_Wafer_AlThickness",
- "Inner_Wafer_PADExternal",
- "Inner_Wafer_PADInternal",
- "Inner_Wafer_GuardRing",
- "Inner_PCB_PortWidth",
- "Inner_PCB_StarboardWidth",
- "Inner_PCB_BevelAngle",
- "Inner_PCB_UpstreamWidth",
- "Inner_PCB_DownstreamWidth",
- "Inner_PCB_MidWidth",
- "Inner_PCB_Thickness",
- "Inner_Wafer_TransverseStrips",
- "Inner_Wafer_LongitudinalStrips",
- "Outer_Wafer_Length",
- "Outer_Wafer_Width",
- "Outer_Wafer_Thickness",
- "Outer_Wafer_AlThickness",
- "Outer_Wafer_PADExternal",
- "Outer_Wafer_PADInternal",
- "Outer_Wafer_GuardRing",
- "Outer_PCB_PortWidth",
- "Outer_PCB_StarboardWidth",
- "Outer_PCB_BevelAngle",
- "Outer_PCB_UpstreamWidth",
- "Outer_PCB_DownstreamWidth",
- "Outer_PCB_MidWidth",
- "Outer_PCB_Thickness",
- "Outer_Wafer_TransverseStrips",
- "Outer_Wafer_LongitudinalStrips",
- "Chamber_Thickness",
- "Chamber_Cylinder_Length",
- "Chamber_Radius",
- "Chamber_ExitTube_Radius",
- "Chamber_ExitTube_Length",
- "Chamber_Flange_Inner_Radius",
- "Chamber_Sphere_Radius",
- "Chamber_Sphere_Shift"
- };
+ vector<string> info = {"Inner_Wafer_Length",
+ "Inner_Wafer_Width",
+ "Inner_Wafer_Thickness",
+ "Inner_Wafer_AlThickness",
+ "Inner_Wafer_PADExternal",
+ "Inner_Wafer_PADInternal",
+ "Inner_Wafer_GuardRing",
+ "Inner_PCB_PortWidth",
+ "Inner_PCB_StarboardWidth",
+ "Inner_PCB_BevelAngle",
+ "Inner_PCB_UpstreamWidth",
+ "Inner_PCB_DownstreamWidth",
+ "Inner_PCB_MidWidth",
+ "Inner_PCB_Thickness",
+ "Inner_Wafer_TransverseStrips",
+ "Inner_Wafer_LongitudinalStrips",
+ "Outer_Wafer_Length",
+ "Outer_Wafer_Width",
+ "Outer_Wafer_Thickness",
+ "Outer_Wafer_AlThickness",
+ "Outer_Wafer_PADExternal",
+ "Outer_Wafer_PADInternal",
+ "Outer_Wafer_GuardRing",
+ "Outer_PCB_PortWidth",
+ "Outer_PCB_StarboardWidth",
+ "Outer_PCB_BevelAngle",
+ "Outer_PCB_UpstreamWidth",
+ "Outer_PCB_DownstreamWidth",
+ "Outer_PCB_MidWidth",
+ "Outer_PCB_Thickness",
+ "Outer_Wafer_TransverseStrips",
+ "Outer_Wafer_LongitudinalStrips",
+ "Chamber_Thickness",
+ "Chamber_Cylinder_Length",
+ "Chamber_Radius",
+ "Chamber_ExitTube_Radius",
+ "Chamber_ExitTube_Length",
+ "Chamber_Flange_Inner_Radius",
+ "Chamber_Sphere_Radius",
+ "Chamber_Sphere_Shift"};
////////////////////
// Inner Detector //
////////////////////
// Wafer parameter
- double Inner_Wafer_Length=-999;
- double Inner_Wafer_Width=-999;
- double Inner_Wafer_Thickness=-999;
- double Inner_Wafer_AlThickness=-999;
- double Inner_Wafer_PADExternal=-999;
- double Inner_Wafer_PADInternal=-999;
- double Inner_Wafer_GuardRing=-999;
+ double Inner_Wafer_Length = -999;
+ double Inner_Wafer_Width = -999;
+ double Inner_Wafer_Thickness = -999;
+ double Inner_Wafer_AlThickness = -999;
+ double Inner_Wafer_PADExternal = -999;
+ double Inner_Wafer_PADInternal = -999;
+ double Inner_Wafer_GuardRing = -999;
// PCB parameter
- double Inner_PCB_PortWidth=-999;
- double Inner_PCB_StarboardWidth=-999;
- double Inner_PCB_BevelAngle=-999;
- double Inner_PCB_UpstreamWidth=-999;
- double Inner_PCB_DownstreamWidth=-999;
- double Inner_PCB_MidWidth=-999;
- double Inner_PCB_Thickness=-999;
- double Inner_Wafer_TransverseStrips=-999;
- double Inner_Wafer_LongitudinalStrips=-999;
+ double Inner_PCB_PortWidth = -999;
+ double Inner_PCB_StarboardWidth = -999;
+ double Inner_PCB_BevelAngle = -999;
+ double Inner_PCB_UpstreamWidth = -999;
+ double Inner_PCB_DownstreamWidth = -999;
+ double Inner_PCB_MidWidth = -999;
+ double Inner_PCB_Thickness = -999;
+ double Inner_Wafer_TransverseStrips = -999;
+ double Inner_Wafer_LongitudinalStrips = -999;
////////////////////
// Outer Detector //
////////////////////
// Wafer parameter
- double Outer_Wafer_Length=-999;
- double Outer_Wafer_Width=-999;
- double Outer_Wafer_Thickness=-999;
- double Outer_Wafer_AlThickness=-999;
- double Outer_Wafer_PADExternal=-999;
- double Outer_Wafer_PADInternal=-999;
- double Outer_Wafer_GuardRing=-999;
+ double Outer_Wafer_Length = -999;
+ double Outer_Wafer_Width = -999;
+ double Outer_Wafer_Thickness = -999;
+ double Outer_Wafer_AlThickness = -999;
+ double Outer_Wafer_PADExternal = -999;
+ double Outer_Wafer_PADInternal = -999;
+ double Outer_Wafer_GuardRing = -999;
// PCB parameter
- double Outer_PCB_PortWidth=-999;
- double Outer_PCB_StarboardWidth=-999;
- double Outer_PCB_BevelAngle=-999;
- double Outer_PCB_UpstreamWidth=-999;
- double Outer_PCB_DownstreamWidth=-999;
- double Outer_PCB_MidWidth=-999;
- double Outer_PCB_Thickness=-999;
- double Outer_Wafer_TransverseStrips=-999;
- double Outer_Wafer_LongitudinalStrips=-999;
+ double Outer_PCB_PortWidth = -999;
+ double Outer_PCB_StarboardWidth = -999;
+ double Outer_PCB_BevelAngle = -999;
+ double Outer_PCB_UpstreamWidth = -999;
+ double Outer_PCB_DownstreamWidth = -999;
+ double Outer_PCB_MidWidth = -999;
+ double Outer_PCB_Thickness = -999;
+ double Outer_Wafer_TransverseStrips = -999;
+ double Outer_Wafer_LongitudinalStrips = -999;
// Vacuum Chamber //
- double Chamber_Thickness=-999;
- double Chamber_Cylinder_Length=-999;
- double Chamber_Radius=-999;
- double Chamber_ExitTube_Radius=-999;
- double Chamber_ExitTube_Length=-999;
- double Chamber_Flange_Inner_Radius=-999;
- double Chamber_Sphere_Radius=-999;
- double Chamber_Sphere_Shift=-999;
-
- if(blocks_info[0]->HasTokenList(info)){
- cout << endl << "//// Strasse info block" << endl;
- Inner_Wafer_Length = blocks_info[0]->GetDouble("Inner_Wafer_Length","mm");
- Inner_Wafer_Width = blocks_info[0]->GetDouble("Inner_Wafer_Width","mm");
- Inner_Wafer_Thickness = blocks_info[0]->GetDouble("Inner_Wafer_Thickness","micrometer");
- Inner_Wafer_AlThickness = blocks_info[0]->GetDouble("Inner_Wafer_AlThickness","micrometer");
- Inner_Wafer_PADExternal = blocks_info[0]->GetDouble("Inner_Wafer_PADExternal","mm");
- Inner_Wafer_PADInternal = blocks_info[0]->GetDouble("Inner_Wafer_PADInternal","mm");
- Inner_Wafer_GuardRing = blocks_info[0]->GetDouble("Inner_Wafer_GuardRing","mm");
- Inner_Wafer_TransverseStrips = blocks_info[0]->GetInt("Inner_Wafer_TransverseStrips");
- Inner_Wafer_LongitudinalStrips = blocks_info[0]->GetInt("Inner_Wafer_LongitudinalStrips");
- Inner_PCB_PortWidth = blocks_info[0]->GetDouble("Inner_PCB_PortWidth","mm");
- Inner_PCB_StarboardWidth = blocks_info[0]->GetDouble("Inner_PCB_StarboardWidth","mm");
- Inner_PCB_BevelAngle = blocks_info[0]->GetDouble("Inner_PCB_BevelAngle","mm");
- Inner_PCB_UpstreamWidth = blocks_info[0]->GetDouble("Inner_PCB_UpstreamWidth","mm");
- Inner_PCB_DownstreamWidth = blocks_info[0]->GetDouble("Inner_PCB_DownstreamWidth","mm");
- Inner_PCB_MidWidth = blocks_info[0]->GetDouble("Inner_PCB_MidWidth","mm");
- Inner_PCB_Thickness = blocks_info[0]->GetDouble("Inner_PCB_Thickness","mm");
- Outer_Wafer_Length = blocks_info[0]->GetDouble("Outer_Wafer_Length","mm");
- Outer_Wafer_Width = blocks_info[0]->GetDouble("Outer_Wafer_Width","mm");
- Outer_Wafer_Thickness = blocks_info[0]->GetDouble("Outer_Wafer_Thickness","mm");
- Outer_Wafer_AlThickness = blocks_info[0]->GetDouble("Outer_Wafer_AlThickness","micrometer");
- Outer_Wafer_PADExternal = blocks_info[0]->GetDouble("Outer_Wafer_PADExternal","mm");
- Outer_Wafer_PADInternal = blocks_info[0]->GetDouble("Outer_Wafer_PADInternal","mm");
- Outer_Wafer_GuardRing = blocks_info[0]->GetDouble("Outer_Wafer_GuardRing","mm");
- Outer_Wafer_TransverseStrips = blocks_info[0]->GetInt("Outer_Wafer_TransverseStrips");
- Outer_Wafer_LongitudinalStrips = blocks_info[0]->GetInt("Outer_Wafer_LongitudinalStrips");
- Outer_PCB_PortWidth = blocks_info[0]->GetDouble("Outer_PCB_PortWidth","mm");
- Outer_PCB_StarboardWidth = blocks_info[0]->GetDouble("Outer_PCB_StarboardWidth","mm");
- Outer_PCB_BevelAngle = blocks_info[0]->GetDouble("Outer_PCB_BevelAngle","deg");
- Outer_PCB_UpstreamWidth = blocks_info[0]->GetDouble("Outer_PCB_UpstreamWidth","mm");
- Outer_PCB_DownstreamWidth = blocks_info[0]->GetDouble("Outer_PCB_DownstreamWidth","mm");
- Outer_PCB_MidWidth = blocks_info[0]->GetDouble("Outer_PCB_MidWidth","mm");
- Outer_PCB_Thickness = blocks_info[0]->GetDouble("Outer_PCB_Thickness","mm");
- Chamber_Thickness= blocks_info[0]->GetDouble("Chamber_Thickness","mm");
- Chamber_Cylinder_Length= blocks_info[0]->GetDouble("Chamber_Cylinder_Length","mm");
- Chamber_Radius= blocks_info[0]->GetDouble("Chamber_Radius","mm");
- Chamber_ExitTube_Radius=blocks_info[0]->GetDouble("Chamber_ExitTube_Radius","mm");
- Chamber_ExitTube_Length=blocks_info[0]->GetDouble("Chamber_ExitTube_Length","mm");
- Chamber_Flange_Inner_Radius=blocks_info[0]->GetDouble("Chamber_Flange_Inner_Radius","mm");
- Chamber_Sphere_Radius=blocks_info[0]->GetDouble("Chamber_Sphere_Radius","mm");
- Chamber_Sphere_Shift=blocks_info[0]->GetDouble("Chamber_Sphere_Shift","mm");
+ double Chamber_Thickness = -999;
+ double Chamber_Cylinder_Length = -999;
+ double Chamber_Radius = -999;
+ double Chamber_ExitTube_Radius = -999;
+ double Chamber_ExitTube_Length = -999;
+ double Chamber_Flange_Inner_Radius = -999;
+ double Chamber_Sphere_Radius = -999;
+ double Chamber_Sphere_Shift = -999;
+
+ if (blocks_info[0]->HasTokenList(info)) {
+ cout << endl << "//// Strasse info block" << endl;
+ Inner_Wafer_Length = blocks_info[0]->GetDouble("Inner_Wafer_Length", "mm");
+ Inner_Wafer_Width = blocks_info[0]->GetDouble("Inner_Wafer_Width", "mm");
+ Inner_Wafer_Thickness = blocks_info[0]->GetDouble("Inner_Wafer_Thickness", "micrometer");
+ Inner_Wafer_AlThickness = blocks_info[0]->GetDouble("Inner_Wafer_AlThickness", "micrometer");
+ Inner_Wafer_PADExternal = blocks_info[0]->GetDouble("Inner_Wafer_PADExternal", "mm");
+ Inner_Wafer_PADInternal = blocks_info[0]->GetDouble("Inner_Wafer_PADInternal", "mm");
+ Inner_Wafer_GuardRing = blocks_info[0]->GetDouble("Inner_Wafer_GuardRing", "mm");
+ Inner_Wafer_TransverseStrips = blocks_info[0]->GetInt("Inner_Wafer_TransverseStrips");
+ Inner_Wafer_LongitudinalStrips = blocks_info[0]->GetInt("Inner_Wafer_LongitudinalStrips");
+ Inner_PCB_PortWidth = blocks_info[0]->GetDouble("Inner_PCB_PortWidth", "mm");
+ Inner_PCB_StarboardWidth = blocks_info[0]->GetDouble("Inner_PCB_StarboardWidth", "mm");
+ Inner_PCB_BevelAngle = blocks_info[0]->GetDouble("Inner_PCB_BevelAngle", "mm");
+ Inner_PCB_UpstreamWidth = blocks_info[0]->GetDouble("Inner_PCB_UpstreamWidth", "mm");
+ Inner_PCB_DownstreamWidth = blocks_info[0]->GetDouble("Inner_PCB_DownstreamWidth", "mm");
+ Inner_PCB_MidWidth = blocks_info[0]->GetDouble("Inner_PCB_MidWidth", "mm");
+ Inner_PCB_Thickness = blocks_info[0]->GetDouble("Inner_PCB_Thickness", "mm");
+ Outer_Wafer_Length = blocks_info[0]->GetDouble("Outer_Wafer_Length", "mm");
+ Outer_Wafer_Width = blocks_info[0]->GetDouble("Outer_Wafer_Width", "mm");
+ Outer_Wafer_Thickness = blocks_info[0]->GetDouble("Outer_Wafer_Thickness", "mm");
+ Outer_Wafer_AlThickness = blocks_info[0]->GetDouble("Outer_Wafer_AlThickness", "micrometer");
+ Outer_Wafer_PADExternal = blocks_info[0]->GetDouble("Outer_Wafer_PADExternal", "mm");
+ Outer_Wafer_PADInternal = blocks_info[0]->GetDouble("Outer_Wafer_PADInternal", "mm");
+ Outer_Wafer_GuardRing = blocks_info[0]->GetDouble("Outer_Wafer_GuardRing", "mm");
+ Outer_Wafer_TransverseStrips = blocks_info[0]->GetInt("Outer_Wafer_TransverseStrips");
+ Outer_Wafer_LongitudinalStrips = blocks_info[0]->GetInt("Outer_Wafer_LongitudinalStrips");
+ Outer_PCB_PortWidth = blocks_info[0]->GetDouble("Outer_PCB_PortWidth", "mm");
+ Outer_PCB_StarboardWidth = blocks_info[0]->GetDouble("Outer_PCB_StarboardWidth", "mm");
+ Outer_PCB_BevelAngle = blocks_info[0]->GetDouble("Outer_PCB_BevelAngle", "deg");
+ Outer_PCB_UpstreamWidth = blocks_info[0]->GetDouble("Outer_PCB_UpstreamWidth", "mm");
+ Outer_PCB_DownstreamWidth = blocks_info[0]->GetDouble("Outer_PCB_DownstreamWidth", "mm");
+ Outer_PCB_MidWidth = blocks_info[0]->GetDouble("Outer_PCB_MidWidth", "mm");
+ Outer_PCB_Thickness = blocks_info[0]->GetDouble("Outer_PCB_Thickness", "mm");
+ Chamber_Thickness = blocks_info[0]->GetDouble("Chamber_Thickness", "mm");
+ Chamber_Cylinder_Length = blocks_info[0]->GetDouble("Chamber_Cylinder_Length", "mm");
+ Chamber_Radius = blocks_info[0]->GetDouble("Chamber_Radius", "mm");
+ Chamber_ExitTube_Radius = blocks_info[0]->GetDouble("Chamber_ExitTube_Radius", "mm");
+ Chamber_ExitTube_Length = blocks_info[0]->GetDouble("Chamber_ExitTube_Length", "mm");
+ Chamber_Flange_Inner_Radius = blocks_info[0]->GetDouble("Chamber_Flange_Inner_Radius", "mm");
+ Chamber_Sphere_Radius = blocks_info[0]->GetDouble("Chamber_Sphere_Radius", "mm");
+ Chamber_Sphere_Shift = blocks_info[0]->GetDouble("Chamber_Sphere_Shift", "mm");
}
-
-
-
vector<NPL::InputBlock*> starget = inParser->GetAllBlocksWithToken("Target");
- double TargetThickness = -999;
- double TargetX = -999;
- double TargetY = -999;
- double TargetZ = -999;
+ double TargetThickness = -999;
+ double TargetX = -999;
+ double TargetY = -999;
+ double TargetZ = -999;
- if(starget.size()==1){
+ if (starget.size() == 1) {
cout << "//// TARGET ////" << endl;
cout << "//// Solid Target found " << endl;
- vector<string> token = {"Thickness","Radius","Material","Angle","X","Y","Z"};
- if(starget[0]->HasTokenList(token)){
- TargetThickness= starget[0]->GetDouble("Thickness","mm");
- TargetX=starget[0]->GetDouble("X","mm");
- TargetY=starget[0]->GetDouble("Y","mm");
- TargetZ=starget[0]->GetDouble("Z","mm");
+ vector<string> token = {"Thickness", "Radius", "Material", "Angle", "X", "Y", "Z"};
+ if (starget[0]->HasTokenList(token)) {
+ TargetThickness = starget[0]->GetDouble("Thickness", "mm");
+ TargetX = starget[0]->GetDouble("X", "mm");
+ TargetY = starget[0]->GetDouble("Y", "mm");
+ TargetZ = starget[0]->GetDouble("Z", "mm");
}
- else{
+ else {
cout << "ERROR: Target token list incomplete, check your input file" << endl;
exit(1);
}
}
-
-
-//////////////////////////////////////////////////
-
-
-double d_TargetFront_InnerActive = 15; //mm
-double d_TargetFront_OuterActive = 43; //mm
-
-double InnerTotalLength =
- Inner_PCB_UpstreamWidth
- + Inner_Wafer_Length * 2
- + Inner_PCB_MidWidth
- + Inner_PCB_DownstreamWidth;
-
-double d_TargetCenter_InnerCenter =
- - TargetThickness/2.
- + d_TargetFront_InnerActive
- - Inner_Wafer_GuardRing
- - Inner_PCB_UpstreamWidth
- + InnerTotalLength/2.;
-
-double OuterTotalLength =
- Outer_PCB_UpstreamWidth
- + Outer_Wafer_Length * 2
- + Outer_PCB_MidWidth
- + Outer_PCB_DownstreamWidth;
-
-double d_TargetCenter_OuterCenter =
- - TargetThickness/2.
- + d_TargetFront_OuterActive
- + OuterTotalLength/2.
- - Outer_Wafer_GuardRing
- - Outer_PCB_UpstreamWidth;
-
-
-cout<<endl;
-cout<< "---------- INPUT DISTANCES (mm) -------------"<<endl;
-cout<<endl;
-cout<<"Target Thickness : "<<TargetThickness<<endl;
-cout<<"Beginning Target - Beginning Inner Active : "<<d_TargetFront_InnerActive<<endl;
-cout<<"Beginning Target - Beginning Outer Active : "<<d_TargetFront_OuterActive<<endl;
-cout<<endl;
-cout<< "--------- CALCULATED DISTANCES (mm) -----------"<<endl;
-
-cout << "InnerTotalLength = "<<InnerTotalLength<<endl;
-cout << "d_TargetCenter_InnerCenter = "<<d_TargetCenter_InnerCenter<<endl;
-cout<<endl;
-
-cout << "OuterTotalLength = "<<OuterTotalLength<<endl;
-cout << "d_TargetCenter_OuterCenter = "<<d_TargetCenter_OuterCenter<<endl;
-cout<<endl;
-cout<< "---------------------------------- -----------"<<endl;
-cout<<"Remark: this calculation assumes that the center of target is at (0,0,0)"<<endl;
-}
+ //////////////////////////////////////////////////
+
+ double d_TargetFront_InnerActive = 15; // mm
+ double d_TargetFront_OuterActive = 43; // mm
+
+ double InnerTotalLength =
+ Inner_PCB_UpstreamWidth + Inner_Wafer_Length * 2 + Inner_PCB_MidWidth + Inner_PCB_DownstreamWidth;
+
+ double d_TargetCenter_InnerCenter = -TargetThickness / 2. + d_TargetFront_InnerActive - Inner_Wafer_GuardRing -
+ Inner_PCB_UpstreamWidth + InnerTotalLength / 2.;
+
+ double OuterTotalLength =
+ Outer_PCB_UpstreamWidth + Outer_Wafer_Length * 2 + Outer_PCB_MidWidth + Outer_PCB_DownstreamWidth;
+
+ double d_TargetCenter_OuterCenter = -TargetThickness / 2. + d_TargetFront_OuterActive + OuterTotalLength / 2. -
+ Outer_Wafer_GuardRing - Outer_PCB_UpstreamWidth;
+
+ cout << endl;
+ cout << "---------- INPUT DISTANCES (mm) -------------" << endl;
+ cout << endl;
+ cout << "Target Thickness : " << TargetThickness << endl;
+ cout << "Beginning Target - Beginning Inner Active : " << d_TargetFront_InnerActive << endl;
+ cout << "Beginning Target - Beginning Outer Active : " << d_TargetFront_OuterActive << endl;
+ cout << endl;
+ cout << "--------- CALCULATED DISTANCES (mm) -----------" << endl;
+
+ cout << "InnerTotalLength = " << InnerTotalLength << endl;
+ cout << "d_TargetCenter_InnerCenter = " << d_TargetCenter_InnerCenter << endl;
+ cout << endl;
+
+ cout << "OuterTotalLength = " << OuterTotalLength << endl;
+ cout << "d_TargetCenter_OuterCenter = " << d_TargetCenter_OuterCenter << endl;
+ cout << endl;
+ cout << "---------------------------------- -----------" << endl;
+ cout << "Remark: this calculation assumes that the center of target is at (0,0,0)" << endl;
+}
diff --git a/Projects/Strasse/macro/CheckAna.cxx b/Projects/Strasse/macro/CheckAna.cxx
index 6e32e08d62582008b20605942d9b6c755090f479..596dfcf434222c19725d165492be1f6636d76e10 100644
--- a/Projects/Strasse/macro/CheckAna.cxx
+++ b/Projects/Strasse/macro/CheckAna.cxx
@@ -1,30 +1,67 @@
{
-gStyle->SetPalette(1);
-//TFile *file= new TFile("../../Outputs/Analysis/PhysicsTree.root");
-TFile *file= new TFile("../../Outputs/Analysis/configJuly2021_ana.root");
-TTree *tree = (TTree*)file->Get("PhysicsTree");
+ gStyle->SetPalette(1);
+ // TFile *file= new TFile("../../Outputs/Analysis/PhysicsTree.root");
+ // TFile *file= new TFile("../../Outputs/Analysis/configJuly2021_ana_1MeV_perfect_Decay.root");
+ // 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);
-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");
+ TTree* tree = (TTree*)file->Get("PhysicsTree");
+ 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;
}
diff --git a/Projects/Strasse/macro/CheckSim.cxx b/Projects/Strasse/macro/CheckSim.cxx
index c7e2a762f830f0799da52eae67f30aa19867125d..519fc722fe8efaa64a4a6b14440303f9675a3aa1 100644
--- a/Projects/Strasse/macro/CheckSim.cxx
+++ b/Projects/Strasse/macro/CheckSim.cxx
@@ -1,80 +1,95 @@
{
-gStyle->SetPalette(1);
-//TFile *file= new TFile("../../Outputs/Simulation/SimulatedTree.root");
-TFile *file= new TFile("../../Outputs/Simulation/configJuly2021.root");
-TTree *tree = (TTree*)file->Get("SimulatedTree");
+ gStyle->SetPalette(1);
+ // TFile *file= new TFile("../../Outputs/Simulation/SimulatedTree.root");
+ // TFile *file= new TFile("../../Outputs/Simulation/configJuly2021_1MeV.root");
+ // 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);
-c1->Divide(4,4);
-c1->cd(1);
-tree->Draw("fRC_Beam_Reaction_Energy:fRC_Vertex_Position_Z>>h1(200,-100,100,650,3200,4500)","","colz");
-c1->cd(2);
-tree->Draw("fRC_Vertex_Position_Y:fRC_Vertex_Position_Z>>h2(200,-100,100,120,-60,60)","","colz");
-c1->cd(3);
-tree->Draw("fRC_Vertex_Position_X:fRC_Vertex_Position_Z>>h3(200,-100,100,120,-60,60)","","colz");
-c1->cd(4);
-tree->Draw("fDetected_Position_Y:fDetected_Position_X>>h4","","colz");
-c1->cd(5);
-tree->Draw("fDetected_Position_Y:fDetected_Position_Z>>h5(900,-70,230,700,-70,70)","","colz");
-c1->cd(6);
-tree->Draw("fInner_TE_StripNbr:fDetected_Position_Z>>h6","","colz");
-c1->cd(7);
-tree->Draw("fInner_LE_StripNbr:fDetected_Position_X>>h7","","colz");
-c1->cd(8);
-tree->Draw("Strasse.GetOuterMultTEnergy()+Strasse.GetInnerMultTEnergy()>>h8(6,0,6)","","");
-cout<<"MULT4="<<h8->GetBinContent(5)<<endl;
-cout<<"Efficiency2p="<<h8->GetBinContent(5)/50000*100<<endl;
-c1->cd(9);
-tree->Draw("fInner_LE_StripNbr>>h9(1250,0,1250)","","");
-c1->cd(10);
-tree->Draw("fInner_TE_StripNbr>>h10(1250,0,1250)","","");
-c1->cd(11);
-tree->Draw("fOuter_LE_StripNbr>>h11(1250,0,1250)","","");
-c1->cd(12);
-tree->Draw("fOuter_TE_StripNbr>>h12(1250,0,1250)","","");
-c1->cd(13);
-tree->Draw("fOuter_TE_StripNbr:fDetected_Position_Z>>h13","","colz");
-c1->cd(14);
-tree->Draw("fOuter_LE_StripNbr:fDetected_Position_X>>h14","","colz");
-c1->cd(15);
-tree->Draw("fRC_Theta[1]>>h15","","");
-tree->Draw("fRC_Theta[1]>>h16","Strasse.GetOuterMultTEnergy()+Strasse.GetInnerMultTEnergy()","same");
-h16->SetLineColor(2);
-h16->Scale(1./4);;
-c1->cd(16);
-tree->Draw("fRC_Theta[0]>>h17","","");
-tree->Draw("fRC_Theta[0]>>h18","Strasse.GetOuterMultTEnergy()+Strasse.GetInnerMultTEnergy()","same");
-h18->SetLineColor(2);
-h18->Scale(1./4);;
+ TCanvas* c1 = new TCanvas("c1", "c1", 1000, 1000);
+ c1->Divide(4, 4);
+ c1->cd(1);
+ tree->Draw("fRC_Beam_Reaction_Energy:fRC_Vertex_Position_Z>>h1(200,-100,100,650,3200,4500)", "", "colz");
+ c1->cd(2);
+ tree->Draw("fRC_Vertex_Position_Y:fRC_Vertex_Position_Z>>h2(200,-100,100,120,-60,60)", "", "colz");
+ c1->cd(3);
+ tree->Draw("fRC_Vertex_Position_X:fRC_Vertex_Position_Z>>h3(200,-100,100,120,-60,60)", "", "colz");
+ c1->cd(4);
+ tree->Draw("fDetected_Position_Y:fDetected_Position_X>>h4", "", "colz");
+ c1->cd(5);
+ tree->Draw("fDetected_Position_Y:fDetected_Position_Z>>h5(900,-70,230,700,-70,70)", "", "colz");
+ c1->cd(6);
+ tree->Draw("fInner_TE_StripNbr:fDetected_Position_Z>>h6", "", "colz");
+ c1->cd(7);
+ tree->Draw("fInner_LE_StripNbr:fDetected_Position_X>>h7", "", "colz");
+ c1->cd(8);
+ tree->Draw("Strasse.GetOuterMultTEnergy()+Strasse.GetInnerMultTEnergy()>>h8(6,0,6)", "", "");
+ cout << "MULT4=" << h8->GetBinContent(5) << endl;
+ cout << "Efficiency2p=" << h8->GetBinContent(5) / h1->Integral() * 100 << endl;
+ c1->cd(9);
+ tree->Draw("fInner_LE_StripNbr>>h9(1250,0,1250)", "", "");
+ c1->cd(10);
+ tree->Draw("fInner_TE_StripNbr>>h10(1250,0,1250)", "", "");
+ c1->cd(11);
+ tree->Draw("fOuter_LE_StripNbr>>h11(1250,0,1250)", "", "");
+ c1->cd(12);
+ tree->Draw("fOuter_TE_StripNbr>>h12(1250,0,1250)", "", "");
+ c1->cd(13);
+ tree->Draw("fOuter_TE_StripNbr:fDetected_Position_Z>>h13", "", "colz");
+ c1->cd(14);
+ tree->Draw("fOuter_LE_StripNbr:fDetected_Position_X>>h14", "", "colz");
+ c1->cd(15);
+ tree->Draw("fRC_Theta[1]>>h15", "", "");
+ tree->Draw("fRC_Theta[1]>>h16", "Strasse.GetOuterMultTEnergy()+Strasse.GetInnerMultTEnergy()", "same");
+ h16->SetLineColor(2);
+ h16->Scale(1. / 4);
+ ;
+ c1->cd(16);
+ tree->Draw("fRC_Theta[0]>>h17", "", "");
+ tree->Draw("fRC_Theta[0]>>h18", "Strasse.GetOuterMultTEnergy()+Strasse.GetInnerMultTEnergy()", "same");
+ h18->SetLineColor(2);
+ h18->Scale(1. / 4);
+ ;
-/*
-TCanvas *c2 = new TCanvas("c2","c2",1000,1000);
-c2->Divide(2,2);
-c2->cd(1);
-tree->Draw("Strasse.GetInnerMultTEnergy()>>hMultInnerT(4,0,4)","","");
-c2->cd(2);
-tree->Draw("Strasse.GetInnerMultLEnergy()>>hMultInnerL(4,0,4)","","");
-c2->cd(3);
-tree->Draw("Strasse.GetOuterMultTEnergy()>>hMultOuterT(4,0,4)","","");
-c2->cd(4);
-tree->Draw("Strasse.GetOuterMultLEnergy()>>hMultOuterL(4,0,4)","","");
+ TCanvas* c2 = new TCanvas("c2", "c2", 1000, 1000);
+ c2->Divide(2, 2);
+ c2->cd(1);
+ tree->Draw("Strasse.GetInnerMultTEnergy()>>hMultInnerT(4,0,4)", "", "");
+ c2->cd(2);
+ tree->Draw("Strasse.GetInnerMultLEnergy()>>hMultInnerL(4,0,4)", "", "");
+ c2->cd(3);
+ tree->Draw("Strasse.GetOuterMultTEnergy()>>hMultOuterT(4,0,4)", "", "");
+ c2->cd(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->Divide(3,2);
-
-c3->cd(1);
-tree->Draw("fRC_Kinetic_Energy[0]:fRC_Theta[0]>>hETheta(900,0,90,1000,0,350)","","colz");
-c3->cd(2);
-tree->Draw("fRC_Kinetic_Energy[0]:fRC_Theta[0]>>hETheta_MultInnerT(900,0,90,1000,0,350)","Strasse.GetInnerMultTEnergy()==1","colz");
-c3->cd(3);
-tree->Draw("fRC_Kinetic_Energy[0]:fRC_Theta[0]>>hETheta_MultInnerL(900,0,90,1000,0,350)","Strasse.GetInnerMultLEnergy()==1","colz");
-c3->cd(4);
-tree->Draw("fRC_Kinetic_Energy[0]:fRC_Theta[0]>>hETheta_MultOuterT(900,0,90,1000,0,350)","Strasse.GetOuterMultTEnergy()==1","colz");
-c3->cd(5);
-tree->Draw("fRC_Kinetic_Energy[0]:fRC_Theta[0]>>hETheta_MultOuterL(900,0,90,1000,0,350)","Strasse.GetOuterMultLEnergy()==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");
-
-*/
+ c3->cd(1);
+ tree->Draw("fRC_Kinetic_Energy[0]:fRC_Theta[0]>>hETheta(900,0,90,1000,0,350)", "", "colz");
+ c3->cd(2);
+ tree->Draw("fRC_Kinetic_Energy[0]:fRC_Theta[0]>>hETheta_MultInnerT(900,0,90,1000,0,350)",
+ "Strasse.GetInnerMultTEnergy()==1", "colz");
+ c3->cd(3);
+ tree->Draw("fRC_Kinetic_Energy[0]:fRC_Theta[0]>>hETheta_MultInnerL(900,0,90,1000,0,350)",
+ "Strasse.GetInnerMultLEnergy()==1", "colz");
+ c3->cd(4);
+ tree->Draw("fRC_Kinetic_Energy[0]:fRC_Theta[0]>>hETheta_MultOuterT(900,0,90,1000,0,350)",
+ "Strasse.GetOuterMultTEnergy()==1", "colz");
+ c3->cd(5);
+ tree->Draw("fRC_Kinetic_Energy[0]:fRC_Theta[0]>>hETheta_MultOuterL(900,0,90,1000,0,350)",
+ "Strasse.GetOuterMultLEnergy()==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");
}