From 3bdceb623af8d9e3cf1d3bea0984eecbac812495 Mon Sep 17 00:00:00 2001
From: moukaddam <mhd.moukaddam@gmail.com>
Date: Tue, 20 Dec 2016 14:15:42 +0000
Subject: [PATCH] Fixing the geometry of the Barrel and Hyball
---
NPLib/Detectors/Tiara/TTiaraBarrelPhysics.cxx | 8 +-
NPLib/Detectors/Tiara/TTiaraHyballPhysics.cxx | 94 ++++++++-----------
Projects/T40/Analysis.cxx | 25 +++--
3 files changed, 60 insertions(+), 67 deletions(-)
diff --git a/NPLib/Detectors/Tiara/TTiaraBarrelPhysics.cxx b/NPLib/Detectors/Tiara/TTiaraBarrelPhysics.cxx
index 3cc316faa..8f24cf727 100644
--- a/NPLib/Detectors/Tiara/TTiaraBarrelPhysics.cxx
+++ b/NPLib/Detectors/Tiara/TTiaraBarrelPhysics.cxx
@@ -431,8 +431,8 @@ TVector3 TTiaraBarrelPhysics::GetPositionOfInteraction(const int i) const{
double Y = INNERBARREL_PCB_Width*(0.5+sin(45*deg));
double Z = Strip_Pos[i]*(0.5*INNERBARREL_ActiveWafer_Length);
//original version of the line above = double Z = (Strip_Pos[i]-0.5)*INNERBARREL_ActiveWafer_Length;
- TVector3 POS(X,Y,-Z);
- POS.RotateZ((DetectorNumber[i]-1)*45*deg);
+ TVector3 POS(X,Y,-Z); // since RowPos = (U-D)/(U+D) => Downstream hit (i.e. Z>0) has RowPos<0, thus the sign
+ POS.RotateZ((5-DetectorNumber[i])*45*deg);// looking downstream Detector 1 is at 3'oclock
return( POS ) ;
}
///////////////////////////////////////////////////////////////////////////////
@@ -440,8 +440,8 @@ TVector3 TTiaraBarrelPhysics::GetRandomisedPositionOfInteraction(const int i) co
TVector3 RandomPOS = GetPositionOfInteraction(i);
TVector3 v1(-12.0, 27.76*(0.5+sin(45*deg)), 0.0); // the numbers used in this line and the one below are related to those in lines 594-597
TVector3 v2(12.0, 27.76*(0.5+sin(45*deg)), 0.0);
- v1.RotateZ((DetectorNumber[i]-1)*45*deg);
- v2.RotateZ((DetectorNumber[i]-1)*45*deg);
+ v1.RotateZ((5-DetectorNumber[i])*45*deg);
+ v2.RotateZ((5-DetectorNumber[i])*45*deg);
TVector3 u = (v2-v1).Unit();
double RandomNumber = Random->Rndm();
TVector3 DeltaHolder((RandomNumber*6.0)-3.0,(RandomNumber*6.0)-3.0,0.0);
diff --git a/NPLib/Detectors/Tiara/TTiaraHyballPhysics.cxx b/NPLib/Detectors/Tiara/TTiaraHyballPhysics.cxx
index 0543e5a8a..9b401c2e1 100644
--- a/NPLib/Detectors/Tiara/TTiaraHyballPhysics.cxx
+++ b/NPLib/Detectors/Tiara/TTiaraHyballPhysics.cxx
@@ -429,7 +429,6 @@ void TTiaraHyballPhysics::ReadConfiguration(NPL::InputParser parser){
double R = blocks[i]->GetDouble("R","mm");
double Phi = blocks[i]->GetDouble("Phi","deg");
AddWedgeDetector(R,Phi,Z);
-
}
else{
@@ -531,52 +530,47 @@ void TTiaraHyballPhysics::InitializeRootOutput(){
///// Specific to TiaraHyballArray ////
void TTiaraHyballPhysics::AddWedgeDetector( double R,double Phi,double Z){
- double Wedge_R_Min = 32.6+R;
- double Wedge_R_Max = 135.1+R;
- double Wedge_Phi_Min = -27.2*deg/rad;
- double Wedge_Phi_Max = 27.2*deg/rad;
- Phi= Phi*deg/rad;
-
- int Wedge_Ring_NumberOfStrip = 16 ;
- int Wedge_Sector_NumberOfStrip = 8 ;
-
- double StripPitchSector = (Wedge_Phi_Max-Wedge_Phi_Min)/Wedge_Sector_NumberOfStrip ;
- double StripPitchRing = (Wedge_R_Max-Wedge_R_Min)/Wedge_Ring_NumberOfStrip ;
-
- TVector3 Strip_1_1;
-
- m_NumberOfDetector++;
- Strip_1_1=TVector3(0,0,Z);
+ m_NumberOfDetector++;
+ double r_min = R+32.6;
+ double r_max = R+135.1;
+ double phi_offset = 5.6*deg;
+ double phi_min = -27.2*deg; // momo: check the offset
+ double phi_max = 27.2*deg;
+ int NumberOfRings = 16 ;
+ int NumberOfSectors = 8 ;
+ double ring_pitch = (r_max-r_min)/NumberOfRings ;
+ double sec_pitch = (phi_max-phi_min)/NumberOfSectors ;
// Buffer object to fill Position Array
- vector<double> lineX ; vector<double> lineY ; vector<double> lineZ ;
-
- vector< vector< double > > OneWedgeStripPositionX ;
- vector< vector< double > > OneWedgeStripPositionY ;
- vector< vector< double > > OneWedgeStripPositionZ ;
-
- TVector3 StripCenter = Strip_1_1;
- for(int f = 0 ; f < Wedge_Ring_NumberOfStrip ; f++){
- lineX.clear() ;
- lineY.clear() ;
- lineZ.clear() ;
-
- for(int b = 0 ; b < Wedge_Sector_NumberOfStrip ; b++){
- StripCenter = Strip_1_1;
- StripCenter.SetY(Wedge_R_Max-f*StripPitchRing-0.5*StripPitchRing);
+ TVector3 StripCenter(0,0,0);
+ vector<double> lineX ;
+ vector<double> lineY ;
+ vector<double> lineZ ;
+ vector< vector< double > > WedgeStripPositionX ;
+ vector< vector< double > > WedgeStripPositionY ;
+ vector< vector< double > > WedgeStripPositionZ ;
+
+ for(int iRing = 0 ; iRing < NumberOfRings ; iRing++){
+ lineX.clear() ;
+ lineY.clear() ;
+ lineZ.clear() ;
+ for(int iSec = 0 ; iSec < NumberOfSectors ; iSec++){ // momo: check the order
+ StripCenter.SetX(0);
+ StripCenter.SetY(r_max - (iRing+0.5)*ring_pitch);
StripCenter.SetZ(Z);
- StripCenter.RotateZ(Phi+Wedge_Phi_Min+b*StripPitchSector);
+ StripCenter.RotateZ(Phi + (iSec-4+0.5)*sec_pitch + phi_offset); // momo : check signs and offset
lineX.push_back( StripCenter.X() );
lineY.push_back( StripCenter.Y() );
lineZ.push_back( StripCenter.Z() );
}
- OneWedgeStripPositionX.push_back(lineX);
- OneWedgeStripPositionY.push_back(lineY);
- OneWedgeStripPositionZ.push_back(lineZ);
+ WedgeStripPositionX.push_back(lineX);
+ WedgeStripPositionY.push_back(lineY);
+ WedgeStripPositionZ.push_back(lineZ);
}
- m_StripPositionX.push_back( OneWedgeStripPositionX ) ;
- m_StripPositionY.push_back( OneWedgeStripPositionY ) ;
- m_StripPositionZ.push_back( OneWedgeStripPositionZ ) ;
+
+ m_StripPositionX.push_back( WedgeStripPositionX ) ;
+ m_StripPositionY.push_back( WedgeStripPositionY ) ;
+ m_StripPositionZ.push_back( WedgeStripPositionZ ) ;
return;
}
@@ -603,20 +597,14 @@ TVector3 TTiaraHyballPhysics::GetPositionOfInteraction(const int i) const{
}
///////////////////////////////////////////////////////////////////////////////////////////////////////////
TVector3 TTiaraHyballPhysics::GetRandomisedPositionOfInteraction(const int i) const{
- TVector3 RandomPosition = GetPositionOfInteraction(i);
- double Zholder = RandomPosition.Z();
- RandomPosition.SetZ(0.0);
- double R = RandomPosition.Mag();
- double Theta = RandomPosition.Theta(); // defines the inclination coordinate in a spherical coordinate system
- double Phi = RandomPosition.Phi(); // defines the azimuthal coordinate in a spherical coordinate system
- double RandomNumber1 = Rand->Rndm();
- double DeltaR = ((RandomNumber1 * 6.4)-3.2);
- R = R + DeltaR; // randomises R within a given detector ring
- double RandomNumber2 = Rand->Rndm();
- double DeltaAngle = ((RandomNumber2 * 0.118682389)-0.0593411946);
- Phi = Phi + DeltaAngle; // randomises Phi within a given detector sector
- RandomPosition.SetXYZ(R*(sin(Theta))*(cos(Phi)),R*(sin(Phi))*(sin(Theta)),Zholder);
- return(RandomPosition) ;
+ TVector3 OriginalPosition = GetPositionOfInteraction(i);
+ double rho = OriginalPosition.Perp();
+ double phi = OriginalPosition.Phi();
+ double z = OriginalPosition.Z();
+ // randomises within a given detector ring and sector
+ double rho_rand = (rho-3.2) + 6.4*sqrt(Rand->Uniform(0,1));// sqrt is necessary for realistic randomise!
+ double phi_rand = phi + Rand->Uniform(-3.4*deg, +3.4*deg);
+ return( TVector3(rho_rand*cos(phi_rand),rho_rand*sin(phi_rand),z) ) ;
}
///////////////////////////////////////////////////////////////////////////////////////////////////////////
diff --git a/Projects/T40/Analysis.cxx b/Projects/T40/Analysis.cxx
index 2c68ec226..f723791e4 100644
--- a/Projects/T40/Analysis.cxx
+++ b/Projects/T40/Analysis.cxx
@@ -99,13 +99,13 @@ void Analysis::Init(){
string beam=NPL::ChangeNameToG4Standard(myReaction->GetNucleus1().GetName());
LightTarget = NPL::EnergyLoss(light+"_"+TargetMaterial+".G4table","G4Table",10 );
LightAl = NPL::EnergyLoss(light+"_Al.G4table","G4Table",10);
- LightSi = NPL::EnergyLoss(light+"_Si.G4table","G4Table",10);
+ //LightSi = NPL::EnergyLoss(light+"_Si.G4table","G4Table",10);
+ LightSi = NPL::EnergyLoss("He4_Si.SRIM","SRIM",10);
BeamTarget = NPL::EnergyLoss(beam+"_"+TargetMaterial+".G4table","G4Table",10);
FinalBeamEnergy = BeamTarget.Slow(OriginalBeamEnergy, TargetThickness*0.5, 0);
myReaction->SetBeamEnergy(FinalBeamEnergy);
cout << "Final Beam energy (middle of target): " << FinalBeamEnergy << endl;
-
Initial = new TInitialConditions();
Rand = new TRandom3();
ThetaNormalTarget = 0 ;
@@ -183,11 +183,11 @@ void Analysis::TreatEvent(){
// Part 2 : Impact Energy
Energy = ELab = 0;
Si_E_TH = TH->Strip_E[countTiaraHyball];
- Energy = Si_E_TH;
+ Energy = Si_E_TH; // calibration for hyball is in MeV
// Correct for energy loss using the thickness of the target and the dead layer
- ELab = LightSi.EvaluateInitialEnergy( Energy*keV ,0.0*micrometer , ThetaTHSurface);
- ELab = LightTarget.EvaluateInitialEnergy( ELab ,TargetThickness/2., ThetaNormalTarget);
+ ELab = LightSi.EvaluateInitialEnergy( Energy ,0.0*micrometer , ThetaTHSurface);
+ //ELab = LightTarget.EvaluateInitialEnergy( ELab ,TargetThickness/2., ThetaNormalTarget);
/////////////////////////////
// Part 3 : Excitation Energy Calculation
Ex = myReaction -> ReconstructRelativistic( ELab , ThetaLab );
@@ -208,7 +208,6 @@ void Analysis::TreatEvent(){
/////////////////////////// LOOP on TiaraBarrel /////////////////////////////
for(unsigned int countTiaraBarrel = 0 ; countTiaraBarrel < TB->Strip_E.size() ; countTiaraBarrel++){
-
/////////////////////////////
// Part 1 : Impact Angle
ThetaTBSurface = 0;
@@ -217,7 +216,10 @@ void Analysis::TreatEvent(){
TVector3 BeamImpact(XTarget,YTarget,0);
TVector3 HitDirection = TB -> GetRandomisedPositionOfInteraction(countTiaraBarrel) - BeamImpact ;
ThetaLab = HitDirection.Angle( BeamDirection );
- ThetaTBSurface = HitDirection.Angle(TVector3(0,0,-1) ) - TMath::PiOver2();
+ TVector3 NormalOnBarrel(+1,0,0); // Normal on detector 5
+ int det = TB->DetectorNumber[countTiaraBarrel];
+ NormalOnBarrel.RotateZ((5-det)*45*deg);
+ ThetaTBSurface = HitDirection.Angle(NormalOnBarrel);
ThetaNormalTarget = HitDirection.Angle( TVector3(0,0,1) ) ;
}
else{
@@ -230,7 +232,7 @@ void Analysis::TreatEvent(){
// Part 2 : Impact Energy
Energy = ELab = 0;
Si_E_InnerTB = TB->Strip_E[countTiaraBarrel];
- Energy = Si_E_InnerTB;
+ Energy = Si_E_InnerTB*keV;// calibration for barrel is in keV
//treat the back detector (in progress)
if(false && TB->Outer_Strip_E.size()>0)
if(TB->Outer_Strip_E[countTiaraBarrel]>0){
@@ -238,9 +240,12 @@ void Analysis::TreatEvent(){
Energy = Si_E_InnerTB + Si_E_OuterTB;
}
+ if (ThetaTBSurface>89*deg && ThetaTBSurface<91*deg ) { //this case produces problems with the Eloss calc
+ //ThetaTBSurface += 2*deg; //temporary solution
+ }
// Evaluate energy using the thickness, Target and Si dead layer Correction
- ELab = LightSi.EvaluateInitialEnergy( Energy*keV ,0.0*micrometer , ThetaTBSurface);
- ELab = LightTarget.EvaluateInitialEnergy( ELab ,TargetThickness/2., ThetaNormalTarget);
+ ELab = LightSi.EvaluateInitialEnergy( Energy ,0.0*micrometer , ThetaTBSurface);
+ //ELab = LightTarget.EvaluateInitialEnergy( ELab ,TargetThickness/2., ThetaNormalTarget);
/////////////////////////////
// Part 3 : Excitation Energy Calculation
--
GitLab