diff --git a/Inputs/DetectorConfiguration/Sharc.detector b/Inputs/DetectorConfiguration/Sharc.detector
index a16b0589e6de301fb0a3f6c53f253af02510b1fa..09cb6c0e37d62bd7ba02cf34cf25627489c9d21a 100644
--- a/Inputs/DetectorConfiguration/Sharc.detector
+++ b/Inputs/DetectorConfiguration/Sharc.detector
@@ -17,29 +17,29 @@ Sharc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Upstream CD
SharcQQQ
- Z= -62
+ Z= -71
R= 0
Phi= 0
ThicknessDector= 500
SharcQQQ
- Z= -62
+ Z= -71
R= 0
Phi= 90
ThicknessDector= 500
SharcQQQ
- Z= -62
+ Z= -71
R= 0
Phi= 180
ThicknessDector= 500
SharcQQQ
- Z= -62
+ Z= -71
R= 0
Phi= 270
ThicknessDector= 500
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Upstream Box
SharcBOX
- Z= -34.3
+ Z= -28.4
ThicknessDector1= 500
ThicknessDector2= 500
ThicknessDector3= 500
@@ -51,34 +51,34 @@ Sharc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Down Stream Box
SharcBOX
- Z= 34.3
- ThicknessDector1= 1000
- ThicknessDector2= 1000
- ThicknessDector3= 1000
- ThicknessDector4= 1000
- ThicknessPAD1= 0
- ThicknessPAD2= 0
- ThicknessPAD3= 0
- ThicknessPAD4= 0
+ Z= 31.4
+ ThicknessDector1= 100
+ ThicknessDector2= 100
+ ThicknessDector3= 100
+ ThicknessDector4= 100
+ ThicknessPAD1= 1000
+ ThicknessPAD2= 1000
+ ThicknessPAD3= 1000
+ ThicknessPAD4= 1000
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Downstream CD
- %SharcQQQ
- Z= 70
+ SharcQQQ
+ Z= 67
R= 0
Phi= 0
ThicknessDector= 500
- %SharcQQQ
- Z= 70
+ SharcQQQ
+ Z= 67
R= 0
Phi= 90
ThicknessDector= 500
- %SharcQQQ
- Z= 70
+ SharcQQQ
+ Z= 67
R= 0
Phi= 180
ThicknessDector= 500
- %SharcQQQ
- Z= 70
+ SharcQQQ
+ Z= 67
R= 0
Phi= 270
ThicknessDector= 500
diff --git a/Inputs/EventGenerator/10He.reaction b/Inputs/EventGenerator/10He.reaction
index a4fedfb4000e66f2d84374576e285bdf2c6d6ed0..e3dfbf77436852de8fcb6c5b48508020d24bab85 100644
--- a/Inputs/EventGenerator/10He.reaction
+++ b/Inputs/EventGenerator/10He.reaction
@@ -45,20 +45,20 @@ TwoBodyReaction
ParticleDecay 10He
Daughter= 9He n
- ExcitationEnergy= 0 0
+ ExcitationEnergy= 0.1 0
DifferentialCrossSection= 11Li(d,3He)10He.txt particle9He
shoot= 1 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-GammaDecay 10He
+GammaDecay 9He
Cascade
BranchingRatio= 30
Energies= 0.1
DifferentialCrossSection= 11Li(d,3He)10He.txt Gamma9He
Cascade
BranchingRatio= 70
- Energies= 0.2 0.3
+ Energies= 0.05 0.5
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/Inputs/EventGenerator/alpha.source b/Inputs/EventGenerator/alpha.source
index 570339771d49ca47acc4fe91fb807bd2b98d6d6e..7b50a3f9d8bbc4b2a3d33e103b2dee120217e4a4 100644
--- a/Inputs/EventGenerator/alpha.source
+++ b/Inputs/EventGenerator/alpha.source
@@ -4,13 +4,14 @@
% Energy are given in MeV , Position in mm %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Isotropic
- EnergyLow= 5
- EnergyHigh= 5
- HalfOpenAngleMin= 140
- HalfOpenAngleMax= 160
- x0= 0
- y0= 0
- z0= 0
- particle= alpha
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ EnergyLow= 3
+ EnergyHigh= 3
+ HalfOpenAngleMin= 0
+ HalfOpenAngleMax= 180
+ x0= 0
+ y0= 0
+ z0= 0
+ Particle= alpha
+ ExcitationEnergy= 0
+
% Supported particle type: proton, neutron, deuton, triton, He3 , alpha
diff --git a/NPAnalysis/macros/GeometricalEfficiency.C b/NPAnalysis/macros/GeometricalEfficiency.C
index 11bda32f9bbda59f1fd22e9f7d4b49723128ac89..abcf7766102eb426f736776eceda4616f458950d 100644
--- a/NPAnalysis/macros/GeometricalEfficiency.C
+++ b/NPAnalysis/macros/GeometricalEfficiency.C
@@ -65,10 +65,10 @@ void GeometricalEfficiency(const char * fname = "myResult"){
tree->SetBranchStatus("InteractionCoordinates", true);
// Prepare histograms
- TH1F *hDetecTheta = new TH1F("hDetecTheta", "DetecTheta", 180, 0, 180);
+ TH1F *hDetecTheta = new TH1F("hDetecTheta", "DetecTheta", 90, 0, 180);
TH1F *hDetecThetaCM = new TH1F("hDetecThetaCM", "hDetecThetaCM", 90, 0, 180);
- TH1F *hEmittTheta = new TH1F("hEmittTheta", "EmittTheta", 180, 0, 180);
- TH1F *hEmittThetaCM = new TH1F("hEmittThetaCM", "hEmittThetaCM", 180, 0, 180);
+ TH1F *hEmittTheta = new TH1F("hEmittTheta", "EmittTheta", 90, 0, 180);
+ TH1F *hEmittThetaCM = new TH1F("hEmittThetaCM", "hEmittThetaCM", 90, 0, 180);
// Read the TTree
Int_t nentries = tree->GetEntries();
@@ -86,39 +86,17 @@ void GeometricalEfficiency(const char * fname = "myResult"){
}
}
- TCanvas* c1 = new TCanvas("c1", "c1");
- // Compute relative efficiency in %
- TH1F *efficiency = new TH1F("hEfficiency", "Efficiency", 180, 0, 180);
- efficiency->SetTitle("Efficiency;#Theta [deg];#epsilon [%]");
- efficiency->Divide(hDetecTheta,hEmittTheta,100);
- efficiency->SetMaximum(110);
- efficiency->Draw();
-
- TCanvas* c2 = new TCanvas("c2", "c2");
- hEmittTheta->Draw();
- hDetecTheta->SetLineColor(kRed);
- hDetecTheta->Draw("same");
-
- TCanvas* c3 = new TCanvas("c3", "Lab Frame (use with source)");
- TH1F* SolidA = new TH1F(*hDetecTheta);
- SolidA->Sumw2();
- TF1* C = new TF1("C",Form("%i /(4*%f)",nentries,M_PI),0,180);
- SolidA->Divide(C,1);
- SolidA->Draw();
- TF1* f = new TF1("f",Form("2 * %f * sin(x*%f/180.) *2*%f/180.",M_PI,M_PI,M_PI),0,180);
- f->Draw("SAME");
- SolidA->GetXaxis()->SetTitle("#theta_{Lab} (deg)");
- SolidA->GetYaxis()->SetTitle("d#Omega (sr) ");
- c3->Update();
-
- TCanvas* c4 = new TCanvas("c4", "CM Frame (use with reaction)");
+ TCanvas* c4 = new TCanvas("c4", "CM Frame");
TH1F* SolidACM = new TH1F(*hDetecThetaCM);
SolidACM->Sumw2();
TF1* CCM = new TF1("CCM",Form("%i /(4*%f)",nentries,M_PI),0,180);
+ TF1* C = new TF1("C",Form("%i /(4*%f)",nentries,M_PI),0,180);
SolidACM->Divide(C,1);
SolidACM->Draw();
+ TF1* f = new TF1("f",Form("2 * %f * sin(x*%f/180.) *2*%f/180.",M_PI,M_PI,M_PI),0,180);
+ f->Draw("SAME");
f->Draw("SAME");
- SolidACM->GetXaxis()->SetTitle("#theta_{Lab} (deg)");
+ SolidACM->GetXaxis()->SetTitle("#theta_{CM} (deg)");
SolidACM->GetYaxis()->SetTitle("d#Omega (sr) ");
c4->Update();
}
diff --git a/NPLib/Physics/NPEnergyLoss.cxx b/NPLib/Physics/NPEnergyLoss.cxx
index c7c3c32fb7a9adaecb7d537d8a3f88535b64e4f6..659705e3d553e8c075f68e07fe3b1b7d4dab980b 100644
--- a/NPLib/Physics/NPEnergyLoss.cxx
+++ b/NPLib/Physics/NPEnergyLoss.cxx
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009-2013 this file is part of the NPTool Project *
+ * Copyright (C) 2009 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
@@ -119,9 +119,9 @@ EnergyLoss::EnergyLoss(string Path , string Source, int NumberOfSlice=100 , int
if ( unit == "MeV" ) energy = energy*MeV;
if ( unit == "GeV" ) energy = energy*GeV;
fEnergy .push_back ( energy );
- fdEdX_Nuclear .push_back ( nuclear );
- fdEdX_Electronic .push_back ( electronic );
- fdEdX_Total .push_back ( nuclear + electronic );
+ fdEdX_Nuclear .push_back ( nuclear*MeV/mm );
+ fdEdX_Electronic .push_back ( electronic*MeV/mm );
+ fdEdX_Total .push_back ( nuclear*MeV/mm + electronic*MeV/mm );
}
// Close File
@@ -317,10 +317,76 @@ double EnergyLoss::EvaluateMaterialThickness( double InitialEnergy , // Ene
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// Evaluate Total Energy of particle from Energy loss in a giver thickness
+double EnergyLoss::EvaluateEnergyFromDeltaE( double DeltaE , // Energy of the detected particle
+ double TargetThickness , // Target Thickness at 0 degree
+ double Angle , // Particle Angle
+ double EnergyMin , // Starting Energy
+ double EnergyMax , // Maximum Energy allowed
+ double EnergyResolution , // Resolution at which function stop
+ int MaxStep ) // Stop after MaxStep Whatever Precision is reached
+ const
+ {
+
+ if (Angle > halfpi) Angle = pi-Angle;
+ TargetThickness = TargetThickness / ( cos(Angle) );
+
+ double step_size = 10.*MeV;
+ double Energy = EnergyMax;
+ double DE = 0 ;
+ bool check_low = false;
+ bool check_high = false ;
+
+ for(int i = 0 ; i < MaxStep ; i++)
+ {
+ DE = Energy - Slow(Energy,TargetThickness,Angle) ;
+ if(abs(DeltaE-DE)<EnergyResolution) return Energy;
+ else if (DeltaE-DE > 0)
+ {
+ if(Energy - step_size > EnergyMin)
+ {
+ Energy = Energy - step_size ;
+ check_low = true ;
+ }
+ else
+ {
+ step_size = step_size / 10.;
+ Energy = Energy - step_size ;
+ }
+
+ }
+
+ else if (DeltaE-DE < 0)
+ {
+ if(Energy + step_size < EnergyMax)
+ {
+ Energy = Energy + step_size ;
+ check_high = true ;
+ }
+ else
+ {
+ step_size = step_size / 10.;
+ Energy = Energy + step_size ;
+ }
+ }
+
+ if(check_high && check_low)
+ {
+ step_size = step_size/10.;
+ check_high = false;
+ check_low = false;
+ }
+
+ if(step_size < EnergyResolution) return Energy;
+
+ }
+
+// cout << "NPL::NPEnergyLoss::EvaluateEnergyFromDeltaE : Max step was reach before requiered resolution was reach " << endl ;
+ return Energy;
+ }
-
-
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
diff --git a/NPLib/Physics/NPEnergyLoss.h b/NPLib/Physics/NPEnergyLoss.h
index 30cf8e73710f0f28eb359d67bedb6b84f76bfee5..c6e8c07c46535926e93a9dec7aa6357481d89660 100644
--- a/NPLib/Physics/NPEnergyLoss.h
+++ b/NPLib/Physics/NPEnergyLoss.h
@@ -1,7 +1,7 @@
#ifndef __EnergyLoss__
#define __EnergyLoss__
/*****************************************************************************
- * Copyright (C) 2009-2013 this file is part of the NPTool Project *
+ * Copyright (C) 2009 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
@@ -71,7 +71,7 @@ namespace NPL
vector<double> fdEdX_Nuclear ; // Nuclear Stopping Power
vector<double> fdEdX_Electronic ; // Electronic Stopping Power
vector<double> fdEdX_Total ; // Total Stopping Power
- Interpolator* fInter ; // Interpolator Used to evaluate Energy loss at given energy
+ Interpolator* fInter ; // Interpolator Used to evaluate Energy loss at given energy
public : // General Function on dE/dX table
double EvaluateNuclearLoss (double ener) const;
@@ -86,7 +86,7 @@ namespace NPL
double Angle ) // Particle Angle
const;
- // Calculate Energy Loss of a particle inside a material
+ // Calculate Energy Loss of a particle inside a material
double EnergyLossCalulation( double Energy , // Energy of the detected particle
double TargetThickness , // Target Thickness at 0 degree
double Angle ) // Particle Angle
@@ -98,6 +98,16 @@ namespace NPL
double TargetThickness , // Target Thickness at 0 degree
double Angle ) // Particle Angle
const ;
+
+ // Evaluate Total Energy of particle from Energy loss in a giver thickness
+ double EvaluateEnergyFromDeltaE( double DeltaE , // Energy of the detected particle
+ double TargetThickness , // Target Thickness at 0 degree
+ double Angle , // Particle Angle
+ double EnergyMin , // Starting Energy
+ double EnergyMax , // Maximum Energy allowed
+ double EnergyResolution , // Resolution at which function stop
+ int MaxStep = 1000000 ) // Stop after MaxStep Whatever Precision is reached
+ const ;
// Evaluate the thickness the particle has been through using the energy loss and initial energy
// usefull for thickness measurement using particle sources
diff --git a/NPLib/Physics/NPReaction.cxx b/NPLib/Physics/NPReaction.cxx
index c2dc6cd616d65b60bceecaf657c7502ef5bc6109..fb7639e0c596066468eaeb23e47b102a51330977 100644
--- a/NPLib/Physics/NPReaction.cxx
+++ b/NPLib/Physics/NPReaction.cxx
@@ -155,12 +155,21 @@ Reaction::Reaction(string reaction){
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Reaction::~Reaction(){
//------------- Default Destructor ------------
+
+ if(fNuclei1)
+ delete fNuclei1;
+
+ if(fNuclei2)
+ delete fNuclei2;
+
+ if(fNuclei3)
+ delete fNuclei3;
+
+ if(fNuclei4)
+ delete fNuclei4;
- delete fNuclei1;
- delete fNuclei2;
- delete fNuclei3;
- delete fNuclei4;
- delete fCrossSectionHist;
+ if(fCrossSectionHist)
+ delete fCrossSectionHist;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
diff --git a/NPLib/Sharc/TSharcData.h b/NPLib/Sharc/TSharcData.h
index 22f72c2330bb71a1f2996c1b664b797543187436..3644260d094ba7ffa485f9c698aba41f48ee5bfd 100644
--- a/NPLib/Sharc/TSharcData.h
+++ b/NPLib/Sharc/TSharcData.h
@@ -29,9 +29,9 @@
using namespace std ;
// ROOT
-#include "TObject.h"
+#include "TNamed.h"
-class TSharcData : public TObject {
+class TSharcData : public TNamed {
private:
// Sharc
// Energy
@@ -40,18 +40,21 @@ private:
vector<Double_t> fSharc_StripFront_Energy;
vector<Double_t> fSharc_StripFront_TimeCFD;
vector<Double_t> fSharc_StripFront_TimeLED;
+ vector<Double_t> fSharc_StripFront_Time;
vector<UShort_t> fSharc_StripBack_DetectorNbr;
vector<UShort_t> fSharc_StripBack_StripNbr;
vector<Double_t> fSharc_StripBack_Energy;
vector<Double_t> fSharc_StripBack_TimeCFD;
vector<Double_t> fSharc_StripBack_TimeLED;
-
+ vector<Double_t> fSharc_StripBack_Time;
+
vector<UShort_t> fSharc_PAD_DetectorNbr;
vector<Double_t> fSharc_PAD_Energy;
vector<Double_t> fSharc_PAD_TimeCFD;
vector<Double_t> fSharc_PAD_TimeLED;
-
+ vector<Double_t> fSharc_PAD_Time;
+
public:
TSharcData();
virtual ~TSharcData();
@@ -66,18 +69,47 @@ public:
inline void SetFront_Energy(const Double_t &Energy){fSharc_StripFront_Energy.push_back(Energy);}
inline void SetFront_TimeCFD(const Double_t &TimeCFD){fSharc_StripFront_TimeCFD.push_back(TimeCFD);}
inline void SetFront_TimeLED(const Double_t &TimeLED){fSharc_StripFront_TimeLED.push_back(TimeLED);}
-
+ inline void SetFront_Time(const Double_t &Time){fSharc_StripFront_Time.push_back(Time);}
+
+
inline void SetBack_DetectorNbr(const UShort_t &DetNbr){fSharc_StripBack_DetectorNbr.push_back(DetNbr);}
inline void SetBack_StripNbr(const UShort_t &StripNbr){fSharc_StripBack_StripNbr.push_back(StripNbr);}
inline void SetBack_Energy(const Double_t &Energy){fSharc_StripBack_Energy.push_back(Energy);}
inline void SetBack_TimeCFD(const Double_t &TimeCFD){fSharc_StripBack_TimeCFD.push_back(TimeCFD);}
inline void SetBack_TimeLED(const Double_t &TimeLED){fSharc_StripBack_TimeLED.push_back(TimeLED);}
-
+ inline void SetBack_Time(const Double_t &Time){fSharc_StripBack_Time.push_back(Time);}
+
+
inline void SetPAD_DetectorNbr(const UShort_t &DetNbr){fSharc_PAD_DetectorNbr.push_back(DetNbr);}
inline void SetPAD_Energy(const Double_t &Energy){fSharc_PAD_Energy.push_back(Energy);}
inline void SetPAD_TimeCFD(const Double_t &TimeCFD){fSharc_PAD_TimeCFD.push_back(TimeCFD);}
inline void SetPAD_TimeLED(const Double_t &TimeLED){fSharc_PAD_TimeLED.push_back(TimeLED);}
+ inline void SetPAD_Time(const Double_t &Time){fSharc_PAD_Time.push_back(Time);}
+
+ inline void SetFront(const UShort_t &DetNbr,const UShort_t &StripNbr,const Double_t &Energy,const Double_t &TimeCFD,const Double_t &TimeLED,const Double_t &Time = 0) {
+ SetFront_DetectorNbr(DetNbr);
+ SetFront_StripNbr(StripNbr);
+ SetFront_Energy(Energy);
+ SetFront_TimeCFD(TimeCFD);
+ SetFront_TimeLED(TimeLED);
+ SetFront_Time(Time);
+ };
+ inline void SetBack(const UShort_t &DetNbr,const UShort_t &StripNbr,const Double_t &Energy,const Double_t &TimeCFD,const Double_t &TimeLED,const Double_t &Time = 0) {
+ SetBack_DetectorNbr(DetNbr);
+ SetBack_StripNbr(StripNbr);
+ SetBack_Energy(Energy);
+ SetBack_TimeCFD(TimeCFD);
+ SetBack_TimeLED(TimeLED);
+ SetBack_Time(Time);
+ };
+ inline void SetPAD(const UShort_t &DetNbr,const Double_t &Energy,const Double_t &TimeCFD,const Double_t &TimeLED,const Double_t &Time = 0) {
+ SetPAD_DetectorNbr(DetNbr);
+ SetPAD_Energy(Energy);
+ SetPAD_TimeCFD(TimeCFD);
+ SetPAD_TimeLED(TimeLED);
+ SetPAD_Time(Time);
+ };
///////////////////// GETTERS ////////////////////////
inline UShort_t GetFront_DetectorNbr(const unsigned int &i) const {return fSharc_StripFront_DetectorNbr[i];}//!
@@ -85,22 +117,28 @@ public:
inline Double_t GetFront_Energy(const unsigned int &i) const {return fSharc_StripFront_Energy[i];}//!
inline Double_t GetFront_TimeCFD(const unsigned int &i) const {return fSharc_StripFront_TimeCFD[i];}//!
inline Double_t GetFront_TimeLED(const unsigned int &i) const {return fSharc_StripFront_TimeLED[i];}//!
+ inline Double_t GetFront_Time(const unsigned int &i) const {return fSharc_StripFront_Time[i];}//!
+
inline UShort_t GetBack_DetectorNbr(const unsigned int &i) const {return fSharc_StripBack_DetectorNbr[i];}//!
inline UShort_t GetBack_StripNbr(const unsigned int &i) const {return fSharc_StripBack_StripNbr[i];}//!
inline Double_t GetBack_Energy(const unsigned int &i) const {return fSharc_StripBack_Energy[i];}//!
inline Double_t GetBack_TimeCFD(const unsigned int &i) const {return fSharc_StripBack_TimeCFD[i];}//!
inline Double_t GetBack_TimeLED(const unsigned int &i) const {return fSharc_StripBack_TimeLED[i];}//!
+ inline Double_t GetBack_Time(const unsigned int &i) const {return fSharc_StripBack_Time[i];}//!
+
inline UShort_t GetPAD_DetectorNbr(const unsigned int &i) const {return fSharc_PAD_DetectorNbr[i];}//!
inline Double_t GetPAD_Energy(const unsigned int &i) const {return fSharc_PAD_Energy[i];}//!
inline Double_t GetPAD_TimeCFD(const unsigned int &i) const {return fSharc_PAD_TimeCFD[i];}//!
inline Double_t GetPAD_TimeLED(const unsigned int &i) const {return fSharc_PAD_TimeLED[i];}//!
+ inline Double_t GetPAD_Time(const unsigned int &i) const {return fSharc_PAD_Time[i];}//!
+
inline unsigned int GetMultiplicityFront() const {return fSharc_StripFront_DetectorNbr.size();}//!
inline unsigned int GetMultiplicityBack() const {return fSharc_StripBack_DetectorNbr.size();}//!
inline unsigned int GetMultiplicityPAD() const {return fSharc_PAD_DetectorNbr.size();}//!
-
+
ClassDef(TSharcData,1) // SharcData structure
};
diff --git a/NPLib/Sharc/TSharcPhysics.cxx b/NPLib/Sharc/TSharcPhysics.cxx
index e2b3d5b5b8f42a3b9338c81aa1af5d5839ad65c4..edc2be07eae868218183aca497ffb0fe05e15639 100644
--- a/NPLib/Sharc/TSharcPhysics.cxx
+++ b/NPLib/Sharc/TSharcPhysics.cxx
@@ -39,16 +39,15 @@ using namespace Sharc_LOCAL;
ClassImp(TSharcPhysics)
///////////////////////////////////////////////////////////////////////////
-TSharcPhysics::TSharcPhysics()
-{
+TSharcPhysics::TSharcPhysics(){
EventMultiplicity = 0 ;
m_EventData = new TSharcData ;
m_PreTreatedData = new TSharcData ;
m_EventPhysics = this ;
m_NumberOfDetector = 0 ;
m_MaximumStripMultiplicityAllowed = 10;
- //m_StripEnergyMatchingSigma = 0.060 ;
- m_StripEnergyMatchingSigma = 50 ;
+ m_StripEnergyMatchingSigma = 0.060 ;
+ //m_StripEnergyMatchingSigma = 10 ;
m_StripEnergyMatchingNumberOfSigma = 3;
// Threshold
@@ -63,36 +62,36 @@ TSharcPhysics::TSharcPhysics()
}
///////////////////////////////////////////////////////////////////////////
-void TSharcPhysics::BuildSimplePhysicalEvent()
-{
+void TSharcPhysics::BuildSimplePhysicalEvent(){
BuildPhysicalEvent();
}
///////////////////////////////////////////////////////////////////////////
-void TSharcPhysics::BuildPhysicalEvent()
-{
+void TSharcPhysics::BuildPhysicalEvent(){
PreTreat();
bool check_PAD = false ;
- if( CheckEvent() == 1 )
- {
+ if( CheckEvent() == 1 ){
vector< TVector2 > couple = Match_Front_Back() ;
EventMultiplicity = couple.size();
- for(unsigned int i = 0 ; i < couple.size() ; ++i)
- {
+ unsigned int size = couple.size();
+ for(unsigned int i = 0 ; i < size ; ++i){
check_PAD = false ;
int N = m_PreTreatedData->GetFront_DetectorNbr(couple[i].X()) ;
+
int Front = m_PreTreatedData->GetFront_StripNbr(couple[i].X()) ;
int Back = m_PreTreatedData->GetBack_StripNbr(couple[i].Y()) ;
+
double Front_E = m_PreTreatedData->GetFront_Energy( couple[i].X() ) ;
double Back_E = m_PreTreatedData->GetBack_Energy( couple[i].Y() ) ;
+
double Front_T = m_PreTreatedData->GetFront_TimeCFD( couple[i].X() ) ;
double Back_T = m_PreTreatedData->GetBack_TimeCFD ( couple[i].Y() ) ;
-
+
DetectorNumber.push_back(N);
StripFront_E.push_back(Front_E);
StripFront_T.push_back(Front_T) ;
@@ -113,7 +112,8 @@ void TSharcPhysics::BuildPhysicalEvent()
Strip_Back.push_back(Back) ;
// Search for associate PAD
- for(unsigned int j = 0 ; j < m_PreTreatedData-> GetMultiplicityPAD() ; ++j){
+ unsigned int sizePAD = m_PreTreatedData-> GetMultiplicityPAD();
+ for(unsigned int j = 0 ; j < sizePAD ; ++j){
if(m_PreTreatedData->GetPAD_DetectorNbr(j)==N){
PAD_E.push_back( m_PreTreatedData-> GetPAD_Energy(j)) ;
PAD_T.push_back( m_PreTreatedData-> GetPAD_TimeCFD(j) ) ;
@@ -122,23 +122,24 @@ void TSharcPhysics::BuildPhysicalEvent()
}
- if(!check_PAD)
- {
+ if(!check_PAD){
PAD_E.push_back(-1000) ;
PAD_T.push_back(-1000) ;
- }
}
}
+ }
+
+ if(DetectorNumber.size()==1)
return;
}
///////////////////////////////////////////////////////////////////////////
-void TSharcPhysics::PreTreat()
-{
+void TSharcPhysics::PreTreat(){
ClearPreTreatedData();
// Front
- for(unsigned int i = 0 ; i < m_EventData->GetMultiplicityFront() ; ++i){
+ unsigned int sizeFront = m_EventData->GetMultiplicityFront();
+ for(unsigned int i = 0 ; i < sizeFront ; ++i){
if( m_EventData->GetFront_Energy(i)>m_StripFront_E_RAW_Threshold && IsValidChannel("Front", m_EventData->GetFront_DetectorNbr(i), m_EventData->GetFront_StripNbr(i)) ){
double Front_E = fStrip_Front_E(m_EventData , i);
if( Front_E > m_StripFront_E_Threshold ){
@@ -153,7 +154,8 @@ void TSharcPhysics::PreTreat()
// Back
- for(unsigned int i = 0 ; i < m_EventData->GetMultiplicityBack() ; ++i){
+ unsigned int sizeBack = m_EventData->GetMultiplicityBack() ;
+ for(unsigned int i = 0 ; i < sizeBack ; ++i){
if( m_EventData->GetBack_Energy(i)>m_StripBack_E_RAW_Threshold && IsValidChannel("Back", m_EventData->GetBack_DetectorNbr(i), m_EventData->GetBack_StripNbr(i)) ){
double Back_E = fStrip_Back_E(m_EventData , i);
if( Back_E > m_StripBack_E_Threshold ){
@@ -167,7 +169,8 @@ void TSharcPhysics::PreTreat()
// PAD
- for(unsigned int i = 0 ; i < m_EventData->GetMultiplicityPAD() ; ++i){
+ unsigned int sizePAD = m_EventData->GetMultiplicityPAD();
+ for(unsigned int i = 0 ; i < sizePAD ; ++i){
if( m_EventData->GetPAD_Energy(i)>m_PAD_E_RAW_Threshold && IsValidChannel("PAD", m_EventData->GetPAD_DetectorNbr(i),1) ){
double PAD_E = fPAD_E(m_EventData , i);
if( PAD_E > m_PAD_E_Threshold ){
@@ -185,8 +188,7 @@ void TSharcPhysics::PreTreat()
///////////////////////////////////////////////////////////////////////////
-int TSharcPhysics :: CheckEvent()
-{
+int TSharcPhysics :: CheckEvent(){
// Check the size of the different elements
if( m_PreTreatedData->GetMultiplicityBack() == m_PreTreatedData->GetMultiplicityFront() )
return 1 ; // Regular Event
@@ -197,8 +199,7 @@ int TSharcPhysics :: CheckEvent()
}
///////////////////////////////////////////////////////////////////////////
-vector < TVector2 > TSharcPhysics :: Match_Front_Back()
-{
+vector < TVector2 > TSharcPhysics :: Match_Front_Back(){
vector < TVector2 > ArrayOfGoodCouple ;
// Prevent code from treating very high multiplicity Event
@@ -224,8 +225,7 @@ vector < TVector2 > TSharcPhysics :: Match_Front_Back()
////////////////////////////////////////////////////////////////////////////
-bool TSharcPhysics :: IsValidChannel(const string DetectorType, const int telescope , const int channel)
-{
+bool TSharcPhysics :: IsValidChannel(const string DetectorType, const int telescope , const int channel){
if(DetectorType == "Front")
return *(m_FrontChannelStatus[telescope-1].begin()+channel-1);
@@ -240,8 +240,7 @@ bool TSharcPhysics :: IsValidChannel(const string DetectorType, const int telesc
}
///////////////////////////////////////////////////////////////////////////
-void TSharcPhysics::ReadAnalysisConfig()
-{
+void TSharcPhysics::ReadAnalysisConfig(){
bool ReadingStatus = false;
// path to file
@@ -404,8 +403,7 @@ void TSharcPhysics::ReadAnalysisConfig()
///////////////////////////////////////////////////////////////////////////
-void TSharcPhysics::Clear()
-{
+void TSharcPhysics::Clear(){
EventMultiplicity = 0;
// Provide a Classification of Event
@@ -433,8 +431,7 @@ void TSharcPhysics::Clear()
//// Innherited from VDetector Class ////
///////////////////////////////////////////////////////////////////////////
-void TSharcPhysics::ReadConfiguration(string Path)
-{
+void TSharcPhysics::ReadConfiguration(string Path){
ifstream ConfigFile ;
ConfigFile.open(Path.c_str()) ;
string LineBuffer ;
@@ -611,8 +608,7 @@ void TSharcPhysics::ReadConfiguration(string Path)
}
///////////////////////////////////////////////////////////////////////////
-void TSharcPhysics::AddParameterToCalibrationManager()
-{
+void TSharcPhysics::AddParameterToCalibrationManager(){
CalibrationManager* Cal = CalibrationManager::getInstance();
for(int i = 0 ; i < m_NumberOfDetector ; ++i)
@@ -642,8 +638,7 @@ void TSharcPhysics::AddParameterToCalibrationManager()
}
///////////////////////////////////////////////////////////////////////////
-void TSharcPhysics::InitializeRootInputRaw()
-{
+void TSharcPhysics::InitializeRootInputRaw(){
TChain* inputChain = RootInput::getInstance()->GetChain() ;
inputChain->SetBranchStatus( "Sharc" , true ) ;
inputChain->SetBranchStatus( "fSharc_*" , true ) ;
@@ -652,27 +647,25 @@ void TSharcPhysics::InitializeRootInputRaw()
}
///////////////////////////////////////////////////////////////////////////
-void TSharcPhysics::InitializeRootInputPhysics()
-{
+void TSharcPhysics::InitializeRootInputPhysics(){
TChain* inputChain = RootInput::getInstance()->GetChain();
inputChain->SetBranchStatus( "EventMultiplicity" , true );
- inputChain->SetBranchStatus( "EventType " , true );
- inputChain->SetBranchStatus( "DetectorNumber " , true );
- inputChain->SetBranchStatus( "Strip_E " , true );
- inputChain->SetBranchStatus( "Strip_T " , true );
- inputChain->SetBranchStatus( "StripFront_E " , true );
- inputChain->SetBranchStatus( "StripBack_T " , true );
- inputChain->SetBranchStatus( "StripFront_E " , true );
- inputChain->SetBranchStatus( "StripBack_T " , true );
- inputChain->SetBranchStatus( "Strip_Front " , true );
- inputChain->SetBranchStatus( "Strip_Back " , true );
- inputChain->SetBranchStatus( "PAD_E " , true );
- inputChain->SetBranchStatus( "PAD_T " , true );
+ inputChain->SetBranchStatus( "EventType" , true );
+ inputChain->SetBranchStatus( "DetectorNumber" , true );
+ inputChain->SetBranchStatus( "Strip_E" , true );
+ inputChain->SetBranchStatus( "Strip_T" , true );
+ inputChain->SetBranchStatus( "StripFront_E" , true );
+ inputChain->SetBranchStatus( "StripBack_T" , true );
+ inputChain->SetBranchStatus( "StripFront_E" , true );
+ inputChain->SetBranchStatus( "StripBack_T" , true );
+ inputChain->SetBranchStatus( "Strip_Front" , true );
+ inputChain->SetBranchStatus( "Strip_Back" , true );
+ inputChain->SetBranchStatus( "PAD_E" , true );
+ inputChain->SetBranchStatus( "PAD_T" , true );
}
///////////////////////////////////////////////////////////////////////////
-void TSharcPhysics::InitializeRootOutput()
-{
+void TSharcPhysics::InitializeRootOutput(){
TTree* outputTree = RootOutput::getInstance()->GetTree();
outputTree->Branch( "Sharc" , "TSharcPhysics" , &m_EventPhysics );
}
@@ -680,32 +673,35 @@ void TSharcPhysics::InitializeRootOutput()
///// Specific to SharcArray ////
-void TSharcPhysics::AddBoxDetector(double Z)
-{
+void TSharcPhysics::AddBoxDetector(double Z){
double BOX_Wafer_Width = 52.20;
- double BOX_Wafer_Length = 76.20;
+ // double BOX_Wafer_Length = 76.20;
+
+ double BOX_ActiveArea_Length = 72;
+ double BOX_ActiveArea_Width = 42;
+
int BOX_Wafer_Front_NumberOfStrip = 24 ;
int BOX_Wafer_Back_NumberOfStrip = 48 ;
- double StripPitchFront = BOX_Wafer_Length/BOX_Wafer_Front_NumberOfStrip ; //mm
- double StripPitchBack = BOX_Wafer_Width/BOX_Wafer_Back_NumberOfStrip ; //mm
+ double StripPitchFront = BOX_ActiveArea_Length/BOX_Wafer_Front_NumberOfStrip ; //mm
+ double StripPitchBack = BOX_ActiveArea_Width/BOX_Wafer_Back_NumberOfStrip ; //mm
TVector3 U; TVector3 V;TVector3 Strip_1_1;
for(int i = 0 ; i < 4 ; i++){
m_NumberOfDetector++;
if(Z<0){// Up Stream
- if(i==0) {U=TVector3(1,0,0);V=TVector3(0,0,1); Strip_1_1=TVector3(-36.,40.5,-60) ;}
- else if(i==1) {U=TVector3(0,1,0);V=TVector3(0,0,1); Strip_1_1=TVector3(-40.5,-36.,-60.) ;}
- else if(i==2) {U=TVector3(-1,0,0);V=TVector3(0,0,1); Strip_1_1=TVector3(36.,-40.5,-60.) ;}
- else if(i==3) {U=TVector3(0,-1,0);V=TVector3(0,0,1); Strip_1_1=TVector3(40.5,36.,-60.) ;}
+ if(i==0) {U=TVector3(1,0,0);V=TVector3(0,0,1); Strip_1_1=TVector3(-36.,40.5,Z-BOX_Wafer_Width/2.) ;}
+ else if(i==1) {U=TVector3(0,1,0);V=TVector3(0,0,1); Strip_1_1=TVector3(-40.5,-36.,Z-BOX_Wafer_Width/2.) ;}
+ else if(i==2) {U=TVector3(-1,0,0);V=TVector3(0,0,1); Strip_1_1=TVector3(36.,-40.5,Z-BOX_Wafer_Width/2.) ;}
+ else if(i==3) {U=TVector3(0,-1,0);V=TVector3(0,0,1); Strip_1_1=TVector3(40.5,36.,Z-BOX_Wafer_Width/2.) ;}
}
if(Z>0){//Down Stream
- if(i==0) {U=TVector3(-1,0,0);V=TVector3(0,0,-1); Strip_1_1=TVector3(36.,40.5,60.) ;}
- else if(i==1) {U=TVector3(0,-1,0);V=TVector3(0,0,-1); Strip_1_1=TVector3(-40.5,36.,60.) ;}
- else if(i==2) {U=TVector3(1,0,0);V=TVector3(0,0,-1); Strip_1_1=TVector3(-36.,-40.5,60.) ;}
- else if(i==3) {U=TVector3(0,1,0);V=TVector3(0,0,-1); Strip_1_1=TVector3(40.5,-36.,60.) ;}
+ if(i==0) {U=TVector3(-1,0,0);V=TVector3(0,0,-1); Strip_1_1=TVector3(36.,40.5,Z+BOX_Wafer_Width/2.) ;}
+ else if(i==1) {U=TVector3(0,-1,0);V=TVector3(0,0,-1); Strip_1_1=TVector3(-40.5,36.,Z+BOX_Wafer_Width/2.) ;}
+ else if(i==2) {U=TVector3(1,0,0);V=TVector3(0,0,-1); Strip_1_1=TVector3(-36.,-40.5,Z+BOX_Wafer_Width/2.) ;}
+ else if(i==3) {U=TVector3(0,1,0);V=TVector3(0,0,-1); Strip_1_1=TVector3(40.5,-36.,Z+BOX_Wafer_Width/2.) ;}
}
// Buffer object to fill Position Array
@@ -794,8 +790,7 @@ void TSharcPhysics::AddQQQDetector( double R,double Phi,double Z){
return;
}
-TVector3 TSharcPhysics::GetDetectorNormal( const int i) const
-{
+TVector3 TSharcPhysics::GetDetectorNormal( const int i) const{
/* TVector3 U = TVector3 ( GetStripPositionX( DetectorNumber[i] , 24 , 1 ) ,
GetStripPositionY( DetectorNumber[i] , 24 , 1 ) ,
GetStripPositionZ( DetectorNumber[i] , 24 , 1 ) )
@@ -820,11 +815,11 @@ TVector3 TSharcPhysics::GetDetectorNormal( const int i) const
}
-TVector3 TSharcPhysics::GetPositionOfInteraction(const int i) const
-{
- TVector3 Position = TVector3 ( GetStripPositionX( DetectorNumber[i] , Strip_Front[i] , Strip_Back[i] ) ,
- GetStripPositionY( DetectorNumber[i] , Strip_Front[i] , Strip_Back[i] ) ,
- GetStripPositionZ( DetectorNumber[i] , Strip_Front[i] , Strip_Back[i] ) ) ;
+TVector3 TSharcPhysics::GetPositionOfInteraction(const int i) const{
+
+ TVector3 Position = TVector3 ( GetStripPositionX( DetectorNumber[i] , Strip_Front[i] , Strip_Back[i] ) ,
+ GetStripPositionY( DetectorNumber[i] , Strip_Front[i] , Strip_Back[i] ) ,
+ GetStripPositionZ( DetectorNumber[i] , Strip_Front[i] , Strip_Back[i] ) ) ;
return(Position) ;
diff --git a/NPLib/scripts/buildliblist.sh b/NPLib/scripts/buildliblist.sh
index 0f83838fd3f23e6f1ba3f34416e1cf0f71349e8b..eeddcd7281559a33e8c5bf50e72b9554a4897965 100755
--- a/NPLib/scripts/buildliblist.sh
+++ b/NPLib/scripts/buildliblist.sh
@@ -52,7 +52,7 @@ do
# add trailing \
name="$name \\"
# add tab at the beginning
- name=$(echo "\t $name")
- echo "\t $name" >> $outfile
+ name=$(echo "\011 $name")
+ echo "\011 $name" >> $outfile
fi ;
done
diff --git a/NPLib/scripts/makefile.sh b/NPLib/scripts/makefile.sh
index d6183c72bf7a6629b0844a9c9f34af4cbf6ff531..8a2120d1e49aad84319b2cbe3ef14b5285035293 100755
--- a/NPLib/scripts/makefile.sh
+++ b/NPLib/scripts/makefile.sh
@@ -49,7 +49,7 @@ do
# only build defined detector libraries
if echo "$detectorlibs" | grep -q "$lname" ; then
# print informations
- echo "\tEntering $name directory..."
+ echo "\011Entering $name directory..."
# add "-C ./" pattern at the beginning of the name
cmd="-C ./$name"
# execute make command with target specified on command line
diff --git a/NPSimulation/Sharc/SharcScorers.cc b/NPSimulation/Sharc/SharcScorers.cc
index 62612c22430461aaf48e5c480ff1dd9fd2ce4810..29efe7a116aed3a883c606899ee1716cb68c68cb 100644
--- a/NPSimulation/Sharc/SharcScorers.cc
+++ b/NPSimulation/Sharc/SharcScorers.cc
@@ -1,4 +1,4 @@
-/*****************************************************************************
+ /*****************************************************************************
* Copyright (C) 2009-2013 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
@@ -80,8 +80,8 @@ G4bool PS_Silicon_Rectangle::ProcessHits(G4Step* aStep, G4TouchableHistory*){
//Rare case where particle is close to edge of silicon plan
- if (m_StripLengthNumber == m_NumberOfStripLength+1) m_StripLengthNumber = m_StripLengthNumber;
- if (m_StripWidthNumber == m_NumberOfStripWidth+1) m_StripWidthNumber = m_StripWidthNumber;
+ if (m_StripLengthNumber > m_NumberOfStripLength) m_StripLengthNumber = m_StripLengthNumber;
+ if (m_StripWidthNumber > m_NumberOfStripWidth) m_StripWidthNumber = m_StripWidthNumber;
m_Index = aStep->GetTrack()->GetTrackID() + m_DetectorNumber * 1e3 + m_StripLengthNumber * 1e6 + m_StripWidthNumber * 1e9;
@@ -183,12 +183,13 @@ G4bool PS_Silicon_Annular::ProcessHits(G4Step* aStep, G4TouchableHistory*){
m_StripThetaNumber = (int) ((m_Position.phi() - m_PhiStart*0.5) / m_StripPitchRadial ) + 1 ;
EnergyAndTime[2] = m_DetectorNumber;
- EnergyAndTime[3] = m_StripRadialNumber;
- EnergyAndTime[4] = m_StripThetaNumber;
+ EnergyAndTime[3] = m_StripThetaNumber;
+ EnergyAndTime[4] = m_StripRadialNumber;
+
//Rare case where particle is close to edge of silicon plan
- if (m_StripRadialNumber == m_NumberOfStripRadial+1) m_StripRadialNumber = m_NumberOfStripRadial;
- if (m_StripThetaNumber == m_NumberOfStripTheta+1) m_StripThetaNumber = m_NumberOfStripTheta;
+ if (m_StripRadialNumber > m_NumberOfStripRadial) m_StripRadialNumber = m_NumberOfStripRadial;
+ if (m_StripThetaNumber > m_NumberOfStripTheta) m_StripThetaNumber = m_NumberOfStripTheta;
m_Index = aStep->GetTrack()->GetTrackID() + m_DetectorNumber * 1e3 + m_StripRadialNumber * 1e6 + m_NumberOfStripTheta * 1e9;
diff --git a/NPSimulation/Tigress/Tigress.cc b/NPSimulation/Tigress/Tigress.cc
index 860b529df2033ebba2b51ecd6e208e8f6ee51092..3685f489dc10b41f62490c4d36f1b3f6a013cbdf 100644
--- a/NPSimulation/Tigress/Tigress.cc
+++ b/NPSimulation/Tigress/Tigress.cc
@@ -48,6 +48,7 @@
#include "G4RunManager.hh"
#include "G4ios.hh"
#include "G4SubtractionSolid.hh"
+#include "G4IntersectionSolid.hh"
#include "G4UnionSolid.hh"
#include "G4ThreeVector.hh"// NPS
#include "Tigress.hh"
@@ -68,11 +69,28 @@ using namespace CLHEP;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
namespace {
- const G4double CrystalOuterRadius = 30.0*mm; // outer radius for crystal
- const G4double CrystalInnerRadius = 5.0*mm; // inner radius for hole in crystal
- const G4double CrystalLength = 90.0*mm; // crystal length
- const G4double CrystalHoleDepth = 15.0*mm; // depth at which starts the hole
+ // Ge crystal
+ // Cylindrical part
+ const G4double CrystalOuterRadius = 30.0*mm; // outer radius for crystal
+ const G4double CrystalInnerRadius = 5.0*mm; // inner radius for hole in crystal
+ const G4double CrystalLength = 90.0*mm; // crystal length
+ const G4double CrystalHoleDepth = 15.0*mm; // depth at which starts the hole
+ const G4double CrystaHoleRadius = 0*cm;
+ const G4double CrystalInterDistance = 0.6*mm; // Distance between two crystal
+
+ // Squared part
+ const G4double CrystalWidth = 56.5*mm; // Width of one crystal
+ const G4double CrystalsShift = 1.05*mm; // this can't be more than 2.75mm. It is the amount by which one side is cut closer to the center than the other
+
+
+ // Bevel part
+ const G4double CrystalBentLegth = 3.62*cm;
+ // this->germanium_corner_cone_end_length = 3.0*cm; // the ending length of the cones
+
+
+
+ // Exogam Stuff
const G4double CrystalEdgeOffset1 = 26.0*mm; // distance of the edge from the center of the crystal
const G4double CrystalEdgeOffset2 = 28.5*mm; // distance of the edge from the center of the crystal
@@ -99,7 +117,7 @@ Tigress::Tigress(){
GreenVisAtt = new G4VisAttributes(G4Colour(0, 1, 0)) ;
RedVisAtt = new G4VisAttributes(G4Colour(1, 0, 0)) ;
WhiteVisAtt = new G4VisAttributes(G4Colour(1, 1, 1)) ;
- TrGreyVisAtt = new G4VisAttributes(G4Colour(0.5, 0.5, 0.5, 0.5)) ;
+ TrGreyVisAtt = new G4VisAttributes(G4Colour(0.5, 0.5, 0.5, 0.5)) ;
m_LogicClover = 0;
@@ -386,80 +404,29 @@ void Tigress::ReadConfiguration(string Path){
// Return a G4VSolid modeling the Crystal
G4VSolid* Tigress::ConstructCrystal(){
- // Copy Past form Exogam collaboration
-
- char sName[40]; // generic for named objects
-
- // define a coaxial shape that will be modify with SubstractSolid
- G4int nbslice = 4;
- G4double zSliceGe[4] = { 0.0*mm, CrystalHoleDepth, CrystalHoleDepth + 3.0*mm, CrystalLength}; // depth where is the hole
- G4double InnRadGe[4] = { 0.0*mm, 0.0*mm, CrystalInnerRadius, CrystalInnerRadius}; // to define the hole in the crystal
- G4double OutRadGe[4] = { CrystalOuterRadius, CrystalOuterRadius, CrystalOuterRadius, CrystalOuterRadius}; // to define the external surface
-
- G4double Edge[3];
-
- sprintf(sName, "coax");
- G4Polycone *coax
- = new G4Polycone(G4String(sName), 0.*deg, 360.*deg, nbslice, zSliceGe, InnRadGe, OutRadGe);
-
- // substract boxes to remove maters.
- G4RotationMatrix rm; // rm.SetName(G4String("RotationEdge"));
-
- // box definition to remove some matter to the crystal
- sprintf(sName, "LongEdge1");
- Edge[0] = (CrystalOuterRadius-CrystalEdgeOffset1); // x half-width
- Edge[1] = 1.001*CrystalOuterRadius; // y half-width
- Edge[2] = 1.001*CrystalLength/2.0; // z half-width
- G4Box *cutEdge1 = new G4Box(G4String(sName),Edge[0],Edge[1],Edge[2]);
-
- sprintf(sName, "LongEdge2");
- Edge[0] = (CrystalOuterRadius-CrystalEdgeOffset2); // x half-width
- Edge[1] = 1.001*CrystalOuterRadius; // y half-width
- Edge[2] = 1.001*CrystalLength/2.0; // z half-width
- G4Box *cutEdge2 = new G4Box(G4String(sName),Edge[0],Edge[1],Edge[2]);
-
- sprintf(sName, "Bevel");
- Edge[0] = 1.001*CrystalOuterRadius;
- Edge[1] = sin(CrystalEdgeAngle)*(CrystalEdgeDepth);
- Edge[2] = 1.001*CrystalLength/2.0;
-
- G4Box *cutBevel = new G4Box(G4String(sName),Edge[0],Edge[1],Edge[2]);
-
- // now remove previously defined box from coax. The box must be placed correctly before
- // since the box definition goes from negative to positive values.
- sprintf(sName, "coax_cut1_edge");
- G4SubtractionSolid *coax_cut1
- = new G4SubtractionSolid (G4String(sName), coax, cutEdge1, &rm, G4ThreeVector(-CrystalOuterRadius,0.0,CrystalLength/2.0));
-
- sprintf(sName, "coax_cut2_edge");
- G4SubtractionSolid *coax_cut2
- = new G4SubtractionSolid (G4String(sName), coax_cut1, cutEdge2, &rm, G4ThreeVector(CrystalOuterRadius,0.0,CrystalLength/2.0));
-
- sprintf(sName, "coax_cut3_edge");
- rm.rotateZ(90.0*deg);
- G4SubtractionSolid *coax_cut3
- = new G4SubtractionSolid (G4String(sName), coax_cut2, cutEdge2, &rm, G4ThreeVector(0.0,CrystalOuterRadius,CrystalLength/2.0));
-
- sprintf(sName, "coax_cut4_edge");
- G4SubtractionSolid *coax_cut4
- = new G4SubtractionSolid (G4String(sName), coax_cut3, cutEdge1, &rm, G4ThreeVector(0.0,-CrystalOuterRadius,CrystalLength/2.0));
- rm.rotateZ(-90.0*deg);
-
- sprintf(sName, "coax_cut5_edge");
- rm.rotateX(CrystalEdgeAngle);
- G4SubtractionSolid *coax_cut5
- = new G4SubtractionSolid (G4String(sName), coax_cut4, cutBevel, &rm, G4ThreeVector(0.,CrystalEdgeOffset2,0.));
- rm.rotateX(-CrystalEdgeAngle);
-
- sprintf(sName, "coax_cut6_edge");
- rm.rotateZ(90.0*deg);
- rm.rotateX(CrystalEdgeAngle);
- G4SubtractionSolid *coax_cut6
- = new G4SubtractionSolid (G4String(sName), coax_cut5, cutBevel, &rm, G4ThreeVector(CrystalEdgeOffset2,0.,0.));
- rm.rotateX(-CrystalEdgeAngle);
- rm.rotateZ(-90.0*deg);
-
- return coax_cut6;
+
+ G4Tubs* Crystal_Cylinder = new G4Tubs("Crystal_Cylinder", 0, CrystalOuterRadius, CrystalLength*0.5, 0, 2*M_PI);
+ G4Tubs* Crystal_Hole = new G4Tubs("Crystal_Hole", 0, CrystalInnerRadius, (CrystalLength-CrystalHoleDepth)*0.5, 0, 2*M_PI);
+
+ G4SubtractionSolid* Crystal_Stage1 = new G4SubtractionSolid("Crystal_Stage1",Crystal_Cylinder,Crystal_Hole,new G4RotationMatrix,G4ThreeVector(0,0,CrystalHoleDepth));
+
+
+ G4Box* Crystal_Box1 = new G4Box("Crystal_Box1", CrystalWidth*0.5, CrystalWidth*0.5,CrystalLength*0.51);
+
+ G4SubtractionSolid* Crystal_Stage2 = new G4SubtractionSolid("Crystal_Stage2",Crystal_Stage1,Crystal_Box1,new G4RotationMatrix,G4ThreeVector(CrystalWidth,0,0));
+ G4SubtractionSolid* Crystal_Stage3 = new G4SubtractionSolid("Crystal_Stage3",Crystal_Stage2,Crystal_Box1,new G4RotationMatrix,G4ThreeVector(-CrystalWidth,0,0));
+ G4SubtractionSolid* Crystal_Stage4 = new G4SubtractionSolid("Crystal_Stage4",Crystal_Stage3,Crystal_Box1,new G4RotationMatrix,G4ThreeVector(0,CrystalWidth,0));
+ G4SubtractionSolid* Crystal_Stage5 = new G4SubtractionSolid("Crystal_Stage5",Crystal_Stage4,Crystal_Box1,new G4RotationMatrix,G4ThreeVector(0,-CrystalWidth,0));
+
+ /*const G4double CrystalOuterRadius = 30.0*mm; // outer radius for crystal
+ const G4double CrystalInnerRadius = 5.0*mm; // inner radius for hole in crystal
+ const G4double CrystalLength = 90.0*mm; // crystal length
+ const G4double CrystalHoleDepth = 15.0*mm; // depth at which starts the hole
+ const G4double CrystalInterDistance = 0.6*mm; // Distance between two quarter crystal
+ const G4double CrystalWidth = 56.5*mm; // Width of one crystal
+ const G4double CrystalsShift = 1.05*mm; // this can't be more than 2.75mm. It is the amount by which one side is cut closer to the
+*/
+ return Crystal_Stage5;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -490,7 +457,7 @@ G4VSolid* Tigress::ConstructCapsule(){
CrystalEdgeOffset1 + CrystalEdgeOffset2 + CrystalToCapsule,
CrystalEdgeOffset1 + CrystalEdgeOffset2 + CrystalToCapsule};
- G4Polyhedra *caps = new G4Polyhedra(G4String("Capsule"), 0.*deg, 360.*deg, 4, nbslice, zSlice, InnRad, OutRad);
+ G4Polyhedra *caps = new G4Polyhedra(G4String("Capsule"), 45.*deg, 360.*deg, 4, nbslice, zSlice, InnRad, OutRad);
return caps;
}
@@ -499,8 +466,7 @@ G4VSolid* Tigress::ConstructCapsule(){
// Return a G4VSolid modeling the BGO
G4VSolid* Tigress::ConstructBGO(){
- return 0;
-
+ return NULL;
}
@@ -513,7 +479,16 @@ void Tigress::ConstructClover(string){
// Construct the clover itself
G4VSolid* Capsule = ConstructCapsule();
// Place the cristal in the clover
- G4ThreeVector CrystalPosition = G4ThreeVector(+113.6*0.5*0.5*mm,+113.6*0.5*0.5*mm,0);
+
+ /*const G4double CrystalOuterRadius = 30.0*mm; // outer radius for crystal
+ const G4double CrystalInnerRadius = 5.0*mm; // inner radius for hole in crystal
+ const G4double CrystalLength = 90.0*mm; // crystal length
+ const G4double CrystalHoleDepth = 15.0*mm; // depth at which starts the hole
+ const G4double CrystalInterDistance = 0.6*mm; // Distance between two quarter crystal*/
+
+ G4ThreeVector CrystalPosition ;
+ double CrystalOffset = (CrystalWidth+CrystalInterDistance)*0.5;
+
G4VSolid* CrystalB = ConstructCrystal();
m_LogicClover =
@@ -532,25 +507,25 @@ void Tigress::ConstructClover(string){
new G4LogicalVolume(CrystalB,m_MaterialVacuum,"LogicCrystalW", 0, 0, 0);
G4RotationMatrix* CrystalRotation = new G4RotationMatrix(0,0,0);
- CrystalPosition = G4ThreeVector(+113.6*0.5*0.5*mm,+113.6*0.5*0.5*mm,0);
+ CrystalPosition = G4ThreeVector(CrystalOffset,CrystalOffset,0);
new G4PVPlacement(G4Transform3D(*CrystalRotation, CrystalPosition),
logicCrystalB,"LogicCrystalB",m_LogicClover,false,m_CloverId[0]);
logicCrystalB->SetVisAttributes(BlueVisAtt);
CrystalRotation->rotate(-180*deg, G4ThreeVector(0,0,1));
- CrystalPosition = G4ThreeVector(+113.6*0.5*0.5*mm,-113.6*0.5*0.5*mm,0);
+ CrystalPosition = G4ThreeVector(+CrystalOffset,-CrystalOffset,0);
new G4PVPlacement(G4Transform3D(*CrystalRotation, CrystalPosition),
logicCrystalG,"LogicCrystalG",m_LogicClover,false,m_CloverId[0]);
logicCrystalG->SetVisAttributes(GreenVisAtt);
CrystalRotation->rotate(-180*deg, G4ThreeVector(0,0,1));
- CrystalPosition = G4ThreeVector(-113.6*0.5*0.5*mm,-113.6*0.5*0.5*mm,0);
+ CrystalPosition = G4ThreeVector(-CrystalOffset,-CrystalOffset,0);
new G4PVPlacement(G4Transform3D(*CrystalRotation, CrystalPosition),
logicCrystalR,"LogicCrystalR",m_LogicClover,false,m_CloverId[0]);
logicCrystalR->SetVisAttributes(RedVisAtt);
CrystalRotation->rotate(-180*deg, G4ThreeVector(0,0,1));
- CrystalPosition = G4ThreeVector(-113.6*0.5*0.5*mm,+113.6*0.5*0.5*mm,0);
+ CrystalPosition = G4ThreeVector(-CrystalOffset,CrystalOffset,0);
new G4PVPlacement(G4Transform3D(*CrystalRotation, CrystalPosition),
logicCrystalW,"LogicCrystalW",m_LogicClover,false,m_CloverId[0]);
logicCrystalW->SetVisAttributes(WhiteVisAtt);
@@ -593,7 +568,6 @@ void Tigress::ConstructDetector(G4LogicalVolume* world){
DetectorRotation->rotate(m_BetaX[i], U);
DetectorRotation->rotate(m_BetaY[i], V);
DetectorRotation->rotate(m_BetaZ[i], W);
- DetectorRotation->rotate(45*deg, W);
G4ThreeVector DetectorPosition = m_R[i]*W;
new G4PVPlacement(G4Transform3D(*DetectorRotation, DetectorPosition),