From 50741c3a314af611c76e011a9adb14088bb3fbf2 Mon Sep 17 00:00:00 2001
From: Clenain <lenain@lpccaen.in2p3.fr>
Date: Wed, 18 Mar 2020 18:00:50 +0100
Subject: [PATCH] updating Minos and Dali geometries and materials, adding new
materials in NPMaterials
---
NPLib/Detectors/Dali/TDaliData.cxx | 3 +
NPLib/Detectors/Dali/TDaliData.h | 13 +-
NPLib/Detectors/Dali/TDaliPhysics.cxx | 60 +-
NPLib/Detectors/Dali/TDaliPhysics.h | 11 +
NPLib/Detectors/Minos/TMinosPhysics.cxx | 1285 +++++++++++------------
NPLib/Detectors/Minos/TMinosPhysics.h | 41 +-
NPLib/TrackReconstruction/Tracking.cxx | 779 +++++++-------
NPLib/TrackReconstruction/Tracking.h | 2 +-
NPSimulation/Core/MaterialManager.cc | 371 ++++---
NPSimulation/Detectors/Dali/Dali.cc | 396 +++----
NPSimulation/Detectors/Dali/Dali.hh | 30 +-
NPSimulation/Detectors/Minos/Minos.cc | 753 +++++--------
NPSimulation/Detectors/Minos/Minos.hh | 40 +-
Projects/Dali/Dali.detector | 81 +-
Projects/Dali/DaliMinos.detector | 0
15 files changed, 1751 insertions(+), 2114 deletions(-)
mode change 100644 => 100755 Projects/Dali/DaliMinos.detector
diff --git a/NPLib/Detectors/Dali/TDaliData.cxx b/NPLib/Detectors/Dali/TDaliData.cxx
index 542de349e..70bc43ecb 100644
--- a/NPLib/Detectors/Dali/TDaliData.cxx
+++ b/NPLib/Detectors/Dali/TDaliData.cxx
@@ -52,10 +52,13 @@ void TDaliData::Clear() {
fDali_T_DetectorNbr.clear();
fDali_TDC.clear();
fDali_Time.clear();
+
+ /* fDali_ParticleID.clear(); */
}
+
//////////////////////////////////////////////////////////////////////
void TDaliData::Dump() const {
// This method is very useful for debuging and worth the dev.
diff --git a/NPLib/Detectors/Dali/TDaliData.h b/NPLib/Detectors/Dali/TDaliData.h
index 43c26549f..6c58a2fd8 100644
--- a/NPLib/Detectors/Dali/TDaliData.h
+++ b/NPLib/Detectors/Dali/TDaliData.h
@@ -44,13 +44,12 @@ class TDaliData : public TObject {
vector<Double_t> fDali_TDC;
vector<Double_t> fDali_Time;
-
+ /* vector<Int_t> fDali_ParticleID; */
//////////////////////////////////////////////////////////////
// Constructor and destructor
public:
TDaliData();
~TDaliData();
-
//////////////////////////////////////////////////////////////
// Inherited from TObject and overriden to avoid warnings
@@ -59,7 +58,6 @@ class TDaliData : public TObject {
void Clear(const Option_t*) {};
void Dump() const;
-
//////////////////////////////////////////////////////////////
// Getters and Setters
// Prefer inline declaration to avoid unnecessary called of
@@ -93,6 +91,11 @@ class TDaliData : public TObject {
fDali_T_DetectorNbr.push_back(DetNbr);
fDali_Time.push_back(Time);
};//!
+
+ // ID for simulation
+ /* inline void SetParticleID(const Int_t& ID) { */
+ /* fDali_ParticleID.push_back(ID); */
+ /* };//! */
// (A&T DC)
inline void SetADCAndTDC(const UShort_t& DetNbr,const Double_t& Energy,const Double_t& Time){
@@ -133,7 +136,9 @@ class TDaliData : public TObject {
inline Double_t Get_Time(const unsigned int &i) const
{return fDali_Time[i];}//!
-
+
+ /* inline Int_t GetParticleID(const unsigned int &i) const */
+ /* {return fDali_ParticleID[i];}//! */
//////////////////////////////////////////////////////////////
// Required for ROOT dictionnary
ClassDef(TDaliData,1) // DaliData structure
diff --git a/NPLib/Detectors/Dali/TDaliPhysics.cxx b/NPLib/Detectors/Dali/TDaliPhysics.cxx
index 11df06300..72f870f57 100644
--- a/NPLib/Detectors/Dali/TDaliPhysics.cxx
+++ b/NPLib/Detectors/Dali/TDaliPhysics.cxx
@@ -75,8 +75,6 @@ void TDaliPhysics::BuildSimplePhysicalEvent() {
BuildPhysicalEvent();
}
-
-
///////////////////////////////////////////////////////////////////////////
void TDaliPhysics::BuildPhysicalEvent() {
// apply thresholds and calibration
@@ -94,6 +92,8 @@ void TDaliPhysics::BuildPhysicalEvent() {
}
}
}
+
+
}
///////////////////////////////////////////////////////////////////////////
@@ -110,6 +110,8 @@ void TDaliPhysics::PreTreat() {
// Energy
unsigned int mysize = m_EventData->GetMultEnergy();
for (UShort_t i = 0; i < mysize ; ++i) {
+
+ /* ParticleID.push_back(m_EventData->GetParticleID(i)); */
if (m_EventData->Get_Energy(i) > m_E_RAW_Threshold) {
Double_t Energy = Cal->ApplyCalibration("Dali/ENERGY"+NPL::itoa(m_EventData->GetE_DetectorNbr(i)),m_EventData->Get_Energy(i));
if (Energy > m_E_Threshold) {
@@ -124,10 +126,11 @@ void TDaliPhysics::PreTreat() {
Double_t Time= Cal->ApplyCalibration("Dali/TIME"+NPL::itoa(m_EventData->GetT_DetectorNbr(i)),m_EventData->Get_Time(i));
m_PreTreatedData->SetTime(m_EventData->GetT_DetectorNbr(i), Time);
}
-}
+}
+
///////////////////////////////////////////////////////////////////////////
void TDaliPhysics::ReadAnalysisConfig() {
bool ReadingStatus = false;
@@ -190,13 +193,12 @@ void TDaliPhysics::ReadAnalysisConfig() {
}
}
-
-
///////////////////////////////////////////////////////////////////////////
void TDaliPhysics::Clear() {
DetectorNumber.clear();
Energy.clear();
Time.clear();
+ /* ParticleID.clear(); */
}
@@ -207,32 +209,32 @@ void TDaliPhysics::ReadConfiguration(NPL::InputParser parser) {
if(NPOptionManager::getInstance()->GetVerboseLevel())
cout << "//// " << blocks.size() << " detectors found " << endl;
- vector<string> cart = {"POS","Shape"};
- vector<string> sphe = {"R","Theta","Phi","Shape"};
+ /* vector<string> cart = {"POS","Shape"}; */
+ /* vector<string> sphe = {"R","Theta","Phi","Shape"}; */
- for(unsigned int i = 0 ; i < blocks.size() ; i++){
- if(blocks[i]->HasTokenList(cart)){
- if(NPOptionManager::getInstance()->GetVerboseLevel())
- cout << endl << "//// Dali " << i+1 << endl;
+ /* for(unsigned int i = 0 ; i < blocks.size() ; i++){ */
+ /* if(blocks[i]->HasTokenList(cart)){ */
+ /* if(NPOptionManager::getInstance()->GetVerboseLevel()) */
+ /* cout << endl << "//// Dali " << i+1 << endl; */
- TVector3 Pos = blocks[i]->GetTVector3("POS","mm");
- string Shape = blocks[i]->GetString("Shape");
- AddDetector(Pos,Shape);
- }
- else if(blocks[i]->HasTokenList(sphe)){
- if(NPOptionManager::getInstance()->GetVerboseLevel())
- cout << endl << "//// Dali " << i+1 << endl;
- double R = blocks[i]->GetDouble("R","mm");
- double Theta = blocks[i]->GetDouble("Theta","deg");
- double Phi = blocks[i]->GetDouble("Phi","deg");
- string Shape = blocks[i]->GetString("Shape");
- AddDetector(R,Theta,Phi,Shape);
- }
- else{
- cout << "ERROR: check your input file formatting " << endl;
- exit(1);
- }
- }
+ /* TVector3 Pos = blocks[i]->GetTVector3("POS","mm"); */
+ /* string Shape = blocks[i]->GetString("Shape"); */
+ /* AddDetector(Pos,Shape); */
+ /* } */
+ /* else if(blocks[i]->HasTokenList(sphe)){ */
+ /* if(NPOptionManager::getInstance()->GetVerboseLevel()) */
+ /* cout << endl << "//// Dali " << i+1 << endl; */
+ /* double R = blocks[i]->GetDouble("R","mm"); */
+ /* double Theta = blocks[i]->GetDouble("Theta","deg"); */
+ /* double Phi = blocks[i]->GetDouble("Phi","deg"); */
+ /* string Shape = blocks[i]->GetString("Shape"); */
+ /* AddDetector(R,Theta,Phi,Shape); */
+ /* } */
+ /* else{ */
+ /* cout << "ERROR: check your input file formatting " << endl; */
+ /* exit(1); */
+ /* } */
+ /* } */
}
///////////////////////////////////////////////////////////////////////////
diff --git a/NPLib/Detectors/Dali/TDaliPhysics.h b/NPLib/Detectors/Dali/TDaliPhysics.h
index 60d4a5cad..39a60b4e3 100644
--- a/NPLib/Detectors/Dali/TDaliPhysics.h
+++ b/NPLib/Detectors/Dali/TDaliPhysics.h
@@ -67,6 +67,7 @@ class TDaliPhysics : public TObject, public NPL::VDetector {
vector<double> Energy;
vector<double> TDC;
vector<double> Time;
+ /* vector<int> ParticleID; */
/// A usefull method to bundle all operation to add a detector
void AddDetector(TVector3 POS, string shape);
@@ -155,6 +156,16 @@ class TDaliPhysics : public TObject, public NPL::VDetector {
TDaliData* GetRawData() const {return m_EventData;}
TDaliData* GetPreTreatedData() const {return m_PreTreatedData;}
+ // Use to access energies et det number
+
+ int GetDetNumber(const int &i) const {return DetectorNumber[i];};
+ int GetMultEnergy() const {return Energy.size();};
+ double GetEnergy(const int &i) const {return Energy[i];};
+ /* int GetParticleID(const int &i) const {return ParticleID[i];}; */
+
+
+
+
// parameters used in the analysis
private:
// thresholds
diff --git a/NPLib/Detectors/Minos/TMinosPhysics.cxx b/NPLib/Detectors/Minos/TMinosPhysics.cxx
index 6435da703..ceda3da52 100644
--- a/NPLib/Detectors/Minos/TMinosPhysics.cxx
+++ b/NPLib/Detectors/Minos/TMinosPhysics.cxx
@@ -53,22 +53,22 @@ TMinosPhysics* current_phy = 0;
///////////////////////////////////////////////////////////////////////////
TMinosPhysics::TMinosPhysics()
- : m_EventData(new TMinosData),
- m_PreTreatedData(new TMinosData),
- m_EventPhysics(this),
- Tracking_functions(new Tracking),
- m_E_RAW_Threshold(0), // adc channels
- m_E_Threshold(0), // MeV
- hfit(new TH1F("hfit","hfit",512,0,512)),
- grxztmp(new TGraph()),
- gryztmp(new TGraph()),
- npoint_temp(0),
- allevt_2pfiltered(0),
- m_NumberOfDetectors(0) {
+ : m_EventData(new TMinosData),
+ m_PreTreatedData(new TMinosData),
+ m_EventPhysics(this),
+ Tracking_functions(new Tracking),
+ m_E_RAW_Threshold(0), // adc channels
+ m_E_Threshold(0), // MeV
+ hfit(new TH1F("hfit","hfit",512,0,512)),
+ grxztmp(new TGraph()),
+ gryztmp(new TGraph()),
+ npoint_temp(0),
+ allevt_2pfiltered(0),
+ m_NumberOfDetectors(0) {
current_phy = this;
data_result.SetClass("TMinosResult");
fitdata.SetClass("TMinosClust");
- }
+ }
int NclusterFit;
@@ -86,7 +86,7 @@ void TMinosPhysics::AddDetector(double R, double Theta, TVector3 POS, double Phi
// Compute the TVector3 corresponding
// Call the cartesian method
}
-
+
///////////////////////////////////////////////////////////////////////////
void TMinosPhysics::BuildSimplePhysicalEvent() {
BuildPhysicalEvent();
@@ -113,736 +113,665 @@ void TMinosPhysics::PreTreat() {
// This method applies thresholds and calibrations
// Isolate 2D tracks with the modified HoughTransformation
// Calculate z for each pad
-
+
data_result.Clear();
fitdata.Clear();
+
+ for(int i = 0; i <2 ; i++){
+ parFit1.push_back(-1000);
+ parFit2.push_back(-1000);
+ parFit3.push_back(-1000);
+ parFit4.push_back(-1000);
+ }
+
+ Xvertex = -1000;
+ Yvertex = -1000;
+ Zvertex = -1000;
+ sumTheta = -1000;
+
+ Theta1 = -1000;
+ Theta2 = -1000;
+ Theta_1= -1000;
+ Theta_2= -1000;
+
+ Phi1 = -1000;
+ Phi2 = -1000;
+ Dmin = -1000;
+
TMinosClust *minosfitdata;
-
+
////////////////////fit variables//////////////////////
-
+
fit_function = new TF1("fit_function",conv_fit, 0, 511, 3);
int fit2DStatus = 0.;
double fit_function_max = 0., fit_function_Tpad = 0.;
double Chi2 = 0.;
double Tshaping = 333.9;
const double TimeBinElec = 30.; // in ns
- double MINOSthresh, VDrift, Z_Shift, DelayTrig;
-
- if(SimulationBool){ // case of simulated data
- MINOSthresh = 1000000000.; //
+ double MINOSthresh, UperFitLim, VDrift, Z_Shift, DelayTrig, ZRot_Minos;
+
+ if(SimulationBool){ // Simulated data
+ UperFitLim = 1000000000.; //
+ MINOSthresh = 899.;
DelayTrig = 0.;
- /* VDrift = 0.036; */
- VDrift = 0.03;
+ VDrift = 0.03475;
Z_Shift = 0;
- }
- else{ // case of experiment data
- MINOSthresh = 100000.;
- DelayTrig = 10170;
- VDrift = 0.03475; // in mm/ns
- /* Z_Shift = 57.16; // */
- Z_Shift = 0. ; //
}
-
+ else{ // Experiment data
+ MINOSthresh = 100;
+ UperFitLim = 100000.;
+ DelayTrig = 1515; // ns, s034 par.
+ ZRot_Minos = 33.81; // Rotation of Minos along Z axis, s034 par.
+ VDrift = 0.03475; // mm/ns, s034 par.
+ Z_Shift = -2.5 ; // mm, s034 par.
+ }
+
///////////////////////////////////////////////////////
ClearPreTreatedData();
/// Read Q(t) signal to fit and extract DriftTime
-
+
unsigned int NbOfPad = m_EventData->GetPadNumberMult();
filled =0;trackNbr=0;
trackNbr_FINAL=0;
- ringsum=0;
- zmax=0.;Iteration=0;
+ ringsum=0;
+ zmax=0.;Iteration=0;
for (unsigned int e = 0; e < NbOfPad; e++){
-
+
x_mm = 0.; y_mm = 0.; z_mm = 0.;maxCharge = 0.; t_pad =0; q_pad=0;
x_mm = m_EventData->GetPadX(e);
y_mm = m_EventData->GetPadY(e);
-
+
if ( !(abs(x_mm)<0.0001 && abs(y_mm)<0.0001) ) {
for (UShort_t o = 0; o < m_EventData->GetTime(e).size() ; ++o){
if(m_EventData->GetCharge(e)[o]> maxCharge){
maxCharge = m_EventData->GetCharge(e)[o];
}
}
-
- if((maxCharge < MINOSthresh && round((sqrt(x_mm*x_mm+y_mm*y_mm)-44.15)/2.1) < 16 && round((sqrt(x_mm*x_mm+y_mm*y_mm)-44.15)/2.1) >4 )) { // to remove rings 1,2,3,4 and 16,17,18
- /* if((maxCharge < MINOSthresh && ( SimulationBool==1 || round((sqrt(x_mm*x_mm+y_mm*y_mm)-44.15)/2.1) < 16 && round((sqrt(x_mm*x_mm+y_mm*y_mm)-43.8)/2.1) >4 )) ) { // to remove rings 1,2,3,4 and 16,17,18 // to not apply the ring condition for simulation*/
- Xpad.push_back(x_mm);
- Ypad.push_back(y_mm);
- Qpad.push_back(maxCharge);
- filled++; //
- }
- }
- } //end of NbOfPad
-
- if(filled>0) { // Nbr of activated pads
- while(Xpad.size()>=10 && Iteration< 20 ) { // !!!!!!!!
- filter_result = 0; // Nbr of pads in the Track
- Iteration++;
- XpadTemp.clear();
- YpadTemp.clear();
- QpadTemp.clear();
- clusterringboolTemp.clear();
-
- filter_result = Tracking_functions->Hough_modified(&Xpad, &Ypad, &Qpad, &XpadTemp, &YpadTemp, &QpadTemp, &clusterringboolTemp);
-
- if(filter_result<0) break;
- if(filter_result>10 && clusterringboolTemp.back()==1){
- /* if(filter_result>4 ){ // modified for the simu w 0 dispersion */
- trackNbr++;
- for(int ik=0; ik<filter_result; ik++) {
- XpadNew.push_back(XpadTemp[ik]);
- YpadNew.push_back(YpadTemp[ik]);
- ZpadNew.push_back(-10000);
- QpadNew.push_back(QpadTemp[ik]);
- clusterringbool.push_back(clusterringboolTemp[ik]);
- clusternbr.push_back(Iteration);
- clusterpads.push_back(filter_result);
- }
+ /* if((maxCharge >= MINOSthresh && round((sqrt(x_mm*x_mm+y_mm*y_mm)-44.15)/2.1) < 14 && round((sqrt(x_mm*x_mm+y_mm*y_mm)-44.15)/2.1) >7 )) { // to remove rings 1,2,3,4 and 16,17,18 */
+ if((maxCharge >= MINOSthresh )) { // to remove rings 1,2,3,4 and 16,17,18
+ /* if((maxCharge >= MINOSthresh && round((sqrt(x_mm*x_mm+y_mm*y_mm)-44.15)/2.1) < 16 && round((sqrt(x_mm*x_mm+y_mm*y_mm)-44.15)/2.1) >4 )) { // to remove rings 1,2,3,4 and 16,17,18 */
+ /* if((maxCharge >= MINOSthresh && ( SimulationBool==1 || (round((sqrt(x_mm*x_mm+y_mm*y_mm)-44.15)/2.1) < 16 && round((sqrt(x_mm*x_mm+y_mm*y_mm)-43.8)/2.1) >4 ))) ) { // to not apply the ring condition for simulation */
+ Xpad.push_back(x_mm);
+ Ypad.push_back(y_mm);
+ Qpad.push_back(maxCharge);
+ filled++; //
}
-
-
- }
- }
-
- for(unsigned int il=0; il<XpadNew.size(); il++){
- minosfitdata = (TMinosClust*)fitdata.ConstructedAt(il);
- minosfitdata->Set(XpadNew[il], YpadNew[il],
- -10000, -10000, QpadNew[il],
- clusternbr[il],clusterpads[il], 0.);
- ZpadNew.push_back(-10000.);
-
- }
-
- if(trackNbr>0 && trackNbr < 5 ){ ///////////
-
- if(filled==0) cerr << "Error !!!" << endl;
- //-------------------------------------------------------
- // STEP 4.1: Fitting taken pads for Qmax and Ttrig info
- //-------------------------------------------------------
-
- for(unsigned int i=0 ; i< NbOfPad; i++ ) {
- hfit->Reset();
- bool fitbool = false;
- x_mm = m_EventData->GetPadX(i);
- y_mm = m_EventData->GetPadY(i);
-
- for(unsigned int jj=0; jj<XpadNew.size(); jj++) {
- if( abs(XpadNew[jj]-x_mm)<0.0001 && abs(YpadNew[jj]-y_mm)<0.0001) {
- fitbool = true;
- indexfill=jj;
- break;
- }
}
-
- if( fitbool==true ) {
-
- /* TGraph gfit( m_EventData->GetTime(i).size() ,&(m_EventData->GetTime(i)[0]),&(m_EventData->GetCharge(i)[0])); */
- /* for (UShort_t o = 0; o < m_EventData->GetTime(i).size() ; ++o){ */
- /* gfit.SetPoint(o,m_EventData->GetTime(i)[o],m_EventData->GetCharge(i)[o]+250); */
- /* if(m_EventData->GetCharge(i)[o]+250> hfit_max){ */
- /* hfit_max_T = m_EventData->GetTime(i)[o]; */
- /* hfit_max=m_EventData->GetCharge(i)[o]+250; */
- /* } */
- /* } */
-
- for(Int_t j=0; j< m_EventData->GetTime(i).size(); j++) {
- if(m_EventData->GetCharge(i)[j]>=0){
- hfit->SetBinContent(hfit->FindBin(m_EventData->GetTime(i)[j]), m_EventData->GetCharge(i)[j]+250);
+ } //end of NbOfPad
+
+ if(filled>0){ // Nbr of activated pads
+ while(Xpad.size()>=7 && Iteration< 20 ) { // !!!!!!!!
+ /* while(Xpad.size()>=10 && Iteration< 20 ) { // !!!!!!!! */
+ filter_result = 0; // Nbr of pads in the Track
+ Iteration++;
+ XpadTemp.clear();
+ YpadTemp.clear();
+ QpadTemp.clear();
+ clusterringboolTemp.clear();
+
+ filter_result = Tracking_functions->Hough_modified(&Xpad, &Ypad, &Qpad, &XpadTemp, &YpadTemp, &QpadTemp, &clusterringboolTemp);
+
+ if(filter_result<0) break;
+ /* if(filter_result>10 && clusterringboolTemp.back()==1){ */
+ /* if(filter_result>4 ){ // modified for the simu w 0 dispersion */
+ if(filter_result>7 ){
+ trackNbr++;
+ for(int ik=0; ik<filter_result; ik++) { // selected pads
+ XpadNew.push_back(XpadTemp[ik]);
+ YpadNew.push_back(YpadTemp[ik]);
+ ZpadNew.push_back(-10000);
+ QpadNew.push_back(QpadTemp[ik]);
+ clusterringbool.push_back(clusterringboolTemp[ik]);
+ clusternbr.push_back(Iteration);
+ clusterpads.push_back(filter_result);
}
}
-
- // Fitting the hfit histogram of last channel if not empty
-
- if(hfit->GetSumOfWeights()>0) {
- hfit->GetXaxis()->SetRange(0,510);
- hfit_max = hfit->GetMaximum();
- hfit_max_T = hfit->GetMaximumBin();
- T_min=-1;
- T_max=-1;
-
- // Find the T_min & T_max limits of the signal non zero
- for(int h=hfit_max_T;h>0;h--) {
- if(T_min == -1 && (hfit->GetBinContent(h))<=250 ) {
- T_min = h;
- break;
- }
- }
- for(int h=hfit_max_T;h<510;h++) {
- if(T_max == -1 && (hfit->GetBinContent(h)) < 1 ) {
- T_max = h;
+ }
+ }
+ for(unsigned int il=0; il<XpadNew.size(); il++){
+ minosfitdata = (TMinosClust*)fitdata.ConstructedAt(il);
+ minosfitdata->Set(XpadNew[il], YpadNew[il],
+ -10000, -10000, QpadNew[il],
+ clusternbr[il],clusterpads[il], 0.);
+ /* ZpadNew.push_back(-10000.); */
+ }
+
+ /* if(trackNbr >= 0 && trackNbr < 5 ){ /////////// */
+ if(trackNbr == 2){ ///////////
+
+ /* if(filled==0) cerr << "Error !!!" << endl; */
+ //-------------------------------------------------------
+ // STEP 4.1: Fitting taken pads for Qmax and Ttrig info
+ //-------------------------------------------------------
+
+ for(unsigned int i=0 ; i< NbOfPad; i++ ) {
+ hfit->Reset();
+ bool fitbool = false;
+ x_mm = m_EventData->GetPadX(i);
+ y_mm = m_EventData->GetPadY(i);
+
+ for(unsigned int jj=0; jj<XpadNew.size(); jj++) {
+ if( abs(XpadNew[jj]-x_mm)<0.0001 && abs(YpadNew[jj]-y_mm)<0.0001) {
+ fitbool = true;
+ indexfill=jj;
break;
}
}
-
- if((hfit_max_T-3.5*(Tshaping/TimeBinElec)) > T_min) T_min = hfit_max_T-2*Tshaping/TimeBinElec;
- if((hfit_max_T+10) < T_max || T_max==-1) T_max = hfit_max_T+10.;
- T_min = max(T_min,0.);
- if(T_max>510) T_max = 510;
-
- // Set fit parameters
- fit_function->SetParameter(0, hfit_max-250.);
- fit_function->SetParameter(1,hfit_max_T - Tshaping/TimeBinElec);
- fit_function->SetParameter(2, Tshaping/TimeBinElec);
- fit_function->SetParLimits(0,0,MINOSthresh);
- fit_function->SetParLimits(1,-20,512);
- fit_function->SetParLimits(2,0,512);
-
- fit2DStatus = hfit->Fit(fit_function,"QN","",T_min,T_max);
- /* fit2DStatus = gfit.Fit(fit_function,"QN","",hfit_max-2*Tshaping,hfit_max+2*Tshaping); */
-
- double fit_function_max = 0., fit_function_Tpad = 0.;
- if(fit2DStatus==0) {
- Chi2 = fit_function->GetChisquare();
- fit_function_max = fit_function->GetMaximum();
- fit_function_Tpad = fit_function->GetParameter(1);
- }
-
- //attribute q_pad and z_mm value
- if(fit2DStatus!=0 || fit_function_max<=20. ||
- fit_function_max >= MINOSthresh || fit_function_Tpad<=0.15 ||
- fit_function_Tpad>=513. || fit_function->GetParameter(2)<=0.15 ||
- fit_function->GetParameter(2)>=513.) {
- q_pad = hfit_max-250.;
- z_mm = -10000;
-
- }
-
- else {
- t_pad = fit_function_Tpad;
- /* z_mm = 300.-(-t_pad*TimeBinElec+DelayTrig)*VDrift[int((sqrt(x_mm*x_mm+y_mm*y_mm)-43.8)/2.1)] */
-
- if(SimulationBool)z_mm = t_pad*TimeBinElec*VDrift; // for simu
- else z_mm = 300.-(-t_pad*TimeBinElec+DelayTrig)*VDrift;
-
- Q_Pad.push_back(fit_function_max-250);
- T_Pad.push_back(t_pad*TimeBinElec);
- X_Pad.push_back(x_mm);
- Y_Pad.push_back(y_mm);
- Z_Pad.push_back(z_mm);
- q_pad = fit_function_max-250.;
-
+ if( fitbool==true ) {
+
+ /* TGraph gfit( m_EventData->GetTime(i).size() ,&(m_EventData->GetTime(i)[0]),&(m_EventData->GetCharge(i)[0])); */
+ /* for (UShort_t o = 0; o < m_EventData->GetTime(i).size() ; ++o){ */
+ /* gfit.SetPoint(o,m_EventData->GetTime(i)[o],m_EventData->GetCharge(i)[o]+250); */
+ /* if(m_EventData->GetCharge(i)[o]+250> hfit_max){ */
+ /* hfit_max_T = m_EventData->GetTime(i)[o]; */
+ /* hfit_max=m_EventData->GetCharge(i)[o]+250; */
+ /* } */
+ /* } */
+
+ for(Int_t j=0; j< m_EventData->GetTime(i).size(); j++) {
+ if(m_EventData->GetCharge(i)[j]>=0){
+ hfit->SetBinContent(hfit->FindBin(m_EventData->GetTime(i)[j]), m_EventData->GetCharge(i)[j]+250);
+ }
}
- ZpadNew[indexfill] = z_mm;
- minosfitdata = (TMinosClust*)fitdata.ConstructedAt(indexfill);
- minosfitdata->Set(XpadNew[indexfill], YpadNew[indexfill], t_pad*TimeBinElec, ZpadNew[indexfill], q_pad, clusternbr[indexfill], clusterpads[indexfill], Chi2);
- QpadNew[indexfill] = q_pad ;
-
- } // end if SumOfWeights > 0
-
- } // end if fitbool==True
+
+ // Fitting the hfit histogram of last channel if not empty
+ if(hfit->GetSumOfWeights()>0) {
+ hfit->GetXaxis()->SetRange(0,510);
+ hfit_max = hfit->GetMaximum();
+ hfit_max_T = hfit->GetMaximumBin();
+ T_min=-1;
+ T_max=-1;
+
+ // Find the T_min & T_max limits of the signal non zero
+ for(int h=hfit_max_T;h>0;h--) {
+ if(T_min == -1 && (hfit->GetBinContent(h))<=250 ) {
+ T_min = h;
+ break;
+ }
+ }
+ for(int h=hfit_max_T;h<510;h++) {
+ if(T_max == -1 && (hfit->GetBinContent(h)) < 1 ) {
+ T_max = h;
+ break;
+ }
+ }
+
+ if((hfit_max_T-3.5*(Tshaping/TimeBinElec)) > T_min) T_min = hfit_max_T-2*Tshaping/TimeBinElec;
+ if((hfit_max_T+10) < T_max || T_max==-1) T_max = hfit_max_T+10.;
+ T_min = max(T_min,0.);
+ if(T_max>510) T_max = 510;
+
+ // Set fit parameters
+ fit_function->SetParameter(0, hfit_max-250.);
+ fit_function->SetParameter(1,hfit_max_T - Tshaping/TimeBinElec);
+ fit_function->SetParameter(2, Tshaping/TimeBinElec);
+ fit_function->SetParLimits(0,0,UperFitLim);
+ fit_function->SetParLimits(1,-20,512);
+ fit_function->SetParLimits(2,0,512);
+
+ fit2DStatus = hfit->Fit(fit_function,"QN","",T_min,T_max);
+ /* fit2DStatus = gfit.Fit(fit_function,"QN","",hfit_max-2*Tshaping,hfit_max+2*Tshaping); */
+
+ double fit_function_max = 0., fit_function_Tpad = 0.;
+ if(fit2DStatus==0) {
+ Chi2 = fit_function->GetChisquare();
+ fit_function_max = fit_function->GetMaximum();
+ fit_function_Tpad = fit_function->GetParameter(1);
+ }
+ //attribute q_pad and z_mm value
+ if(fit2DStatus!=0 || fit_function_max<=20. ||
+ fit_function_max >= UperFitLim || fit_function_Tpad<=0.15 ||
+ fit_function_Tpad>=513. || fit_function->GetParameter(2)<=0.15 ||
+ fit_function->GetParameter(2)>=513.) {
+ q_pad = hfit_max-250.;
+ z_mm = -10000;
+ }
+ else {
+ t_pad = fit_function_Tpad;
+ /* z_mm = 300.-(-t_pad*TimeBinElec+DelayTrig)*VDrift[int((sqrt(x_mm*x_mm+y_mm*y_mm)-43.8)/2.1)] */
+ if(SimulationBool)z_mm = t_pad*TimeBinElec*VDrift; // for simu
+ else z_mm =(t_pad*TimeBinElec-DelayTrig)*VDrift;
+ /* else z_mm = 300.-(-t_pad*TimeBinElec+DelayTrig)*VDrift; */
+ Q_Pad.push_back(fit_function_max-250);
+ T_Pad.push_back(t_pad*TimeBinElec);
+ X_Pad.push_back(x_mm);
+ Y_Pad.push_back(y_mm);
+ Z_Pad.push_back(z_mm);
+ q_pad = fit_function_max-250.;
+ }
+ ZpadNew[indexfill] = z_mm + Z_Shift; // CYRIL Test
+ /* ZpadNew[indexfill] = z_mm; */
+ QpadNew[indexfill] = q_pad;
+ minosfitdata = (TMinosClust*)fitdata.ConstructedAt(indexfill);
+ minosfitdata->Set(XpadNew[indexfill], YpadNew[indexfill], t_pad*TimeBinElec, z_mm, q_pad, clusternbr[indexfill], clusterpads[indexfill], Chi2);
+ } // end if SumOfWeights > 0
+ } // end if fitbool==True
else continue;
- } // end of NbOfPad
-
- //-------------------------------------------------------
- // STEP 3.2: Filtering the tracks off possible noise
- // with Hough3D (3*2D planes)
- //-------------------------------------------------------
-
- cluster_temp = 0;
- int ringtouch[19]={0};
- /* for(int ko=1; ko<19; ko++) ringtouch[ko] = 0; ///// Added by Cyril */
-
- for(unsigned int i=0;i<(XpadNew.size());i++){
- if(xin.size()>0 && ((cluster_temp!=int(clusternbr[i]) && i!=0) || i==(XpadNew.size() - 1))){
-
- Tracking_functions->Hough_3D(&xin, &yin, &zin, &qin, &xout, &yout, &zout, &qout);
-
- for(unsigned int ij=0; ij<xout.size();ij++){
- if(zout[ij]>zmax) {zmax = zout[ij];}
- ringtouch[int(round((sqrt(xout[ij]*xout[ij]+yout[ij]*yout[ij])-44.15)/2.1))]++; // Corr by Cyril
- }
-
- for(int ko=0; ko<19; ko++){
- if(ringtouch[ko]>0) ringsum++;
- }
-
- /* if(zmax> 290 ) ringsum = 16; //comment it if z_mm is not well reconstruct by Vdrift and DelayTrig ! */
-
- // Tracks of interest: >10 pads and >=12 rings hit
-
- /* if(xout.size()>10 &&(SimulationBool==1 || ringsum>=8)){ */
- if(xout.size()>10 && ringsum>=8){
+ } // end of NbOfPad
+
+ //-------------------------------------------------------
+ // STEP 3.2: Filtering the tracks off possible noise
+ // with Hough3D (3*2D planes)
+ //-------------------------------------------------------
+
+ cluster_temp = 0;
+ int ringtouch[19]={0};
+ /* for(int ko=1; ko<19; ko++) ringtouch[ko] = 0; ///// Added by Cyril */
+ for(unsigned int i=0;i<(XpadNew.size());i++){
+ if(xin.size()>0 && ((cluster_temp != int(clusternbr[i]) && i!=0) || i==(XpadNew.size()-1))){ // We fill xin until the next cluster
+ Tracking_functions->Hough_3D(&xin, &yin, &zin, &qin, &xout, &yout, &zout, &qout);
+ for(unsigned int ij=0; ij<xout.size();ij++){
+ if(zout[ij]>zmax) {zmax = zout[ij];}
+ ringtouch[int(round((sqrt(xout[ij]*xout[ij]+yout[ij]*yout[ij])-44.15)/2.1))]++; // Corr by Cyril
+ }
+ for(int ko=0; ko<19; ko++){
+ if(ringtouch[ko]>0) ringsum++;
+ }
+
+ // Tracks of interest: >10 pads and >=12 rings hit
+ if(xout.size()>10 && ringsum>=8){
+
+ npoint=0;
+ trackNbr_FINAL++;
+ double charge_temp=0.;
+ double lenght_temp=0;
+ double zintarget=300;
+ double zouttarget=0;
+ int indexin=0; int indexout=0;
+ for(unsigned int ij=0; ij<xout.size(); ij++){
+
+ point.SetXYZ(xout[ij],yout[ij],zout[ij]);
+ if(!SimulationBool)point.RotateZ(ZRot_Minos*TMath::DegToRad());//15.6*TMath::DegToRad() CHANGE####
+ xoutprime.push_back(point.X());youtprime.push_back(point.Y());zoutprime.push_back(point.Z());
+ TOTxoutprime.push_back(point.X());
+ // Added line to see events used for fit of lines in tree
+ TOTyoutprime.push_back(point.Y());
+ TOTzoutprime.push_back(point.Z());
+ TOTqout.push_back(qout[ij]);
+ // save cluster_temp
+ // save numbtrack
+ trackclusternbr.push_back(clusternbr[i]);
+ tracknbr.push_back( trackNbr_FINAL ) ;
+
+ grxztmp->SetPoint(npoint,zoutprime[ij],xoutprime[ij]);
+ gryztmp->SetPoint(npoint,zoutprime[ij],youtprime[ij]);
+ charge_temp += qout[ij];
+ minosdata_result = (TMinosResult*)data_result.ConstructedAt(array_final);
+ minosdata_result->Set(xoutprime[ij], youtprime[ij], zoutprime[ij], qout[ij], trackNbr_FINAL, xout.size(), zmax);
+ array_final++;
+ npoint++;
+
+ if(zoutprime[ij]<zintarget) {zintarget=zoutprime[ij];indexin=ij;}
+ if(zoutprime[ij]>zouttarget) {zouttarget=zoutprime[ij];indexout=ij;}
+ }
+
+ if(xout.size()>0)lenght_temp=sqrt(pow(zoutprime[indexin]-zoutprime[indexout],2)+pow(youtprime[indexin]-youtprime[indexout],2)+pow(xoutprime[indexin]-xoutprime[indexout],2));
+
+ grxz.push_back(*grxztmp);
+ gryz.push_back(*gryztmp);
+ grxztmp->Set(0);
+ /* if(maxCharge >= MINOSthresh ) { */
+ gryztmp->Set(0);
+
+ chargeTot.push_back(charge_temp);
+ lenght.push_back(lenght_temp);
+ } // end of if(xout.size()>10 &&)
+ //X_mm->Fill();
+ xin.clear();
+ yin.clear();
+ zin.clear();
+ qin.clear();
+ xout.clear();
+ yout.clear();
+ zout.clear();
+ xoutprime.clear();
+ youtprime.clear();
+ zoutprime.clear();
+ qout.clear();
+ npoint_temp=0;
+ ringsum=0;// angle between 1st track and z axis in 3D in degrees
+ zmax=0.;
+ for(int ko=0; ko<18; ko++) ringtouch[ko] = 0;
+
+ } // if(xin.size()>0 && (( cluster_temp .... )
- npoint=0;
- trackNbr_FINAL++;
- double charge_temp=0.;
- double lenght_temp=0;
- double zintarget=300;
- double zouttarget=0;
- int indexin=0; int indexout=0;
- for(unsigned int ij=0; ij<xout.size(); ij++){
-
- point.SetXYZ(xout[ij],yout[ij],zout[ij]);
- point.RotateZ(0.);//15.6*TMath::DegToRad() CHANGE####
- xoutprime.push_back(point.X());youtprime.push_back(point.Y());zoutprime.push_back(point.Z());
- TOTxoutprime.push_back(point.X());
- // Added line to see events used for fit of lines in tree
- TOTyoutprime.push_back(point.Y());
- TOTzoutprime.push_back(point.Z());
- TOTqout.push_back(qout[ij]);
- // save cluster_temp
- // save numbtrack
- trackclusternbr.push_back(clusternbr[i]);
- tracknbr.push_back( trackNbr_FINAL ) ;
-
- grxztmp->SetPoint(npoint,zoutprime[ij],xoutprime[ij]);
- gryztmp->SetPoint(npoint,zoutprime[ij],youtprime[ij]);
- charge_temp += qout[ij];
- minosdata_result = (TMinosResult*)data_result.ConstructedAt(array_final);
- minosdata_result->Set(xoutprime[ij], youtprime[ij], zoutprime[ij], qout[ij], trackNbr_FINAL, xout.size(), zmax);
- array_final++;
- npoint++;
-
- ///////// to exclude data w/ rings from 15to18
- if(sqrt(xout[ij]*xout[ij]+yout[ij]*yout[ij])<70 && sqrt(xout[ij]*xout[ij]+yout[ij]*yout[ij])> 53 ){
- if(zoutprime[ij]<zintarget) {zintarget=zoutprime[ij];indexin=ij;}
- if(zoutprime[ij]>zouttarget) {zouttarget=zoutprime[ij];indexout=ij;}
+ cluster_temp = clusternbr[i]; // Number of the track/cluster
+ /* if(!(clusterpads[i]>=10 && clusterringbool[i]==1 && ZpadNew[i]>-100 && ZpadNew[i]<=310)) continue; // Warning <=310 necessary ? (Cyril) */
+ if(!(clusterpads[i]>=10 && clusterringbool[i]==1 && ZpadNew[i]>-100 && ZpadNew[i]<=310)) continue; // Warning <=310 necessary ? (Cyril)
+ else
+ {
+ xin.push_back(XpadNew[i]);
+ yin.push_back(YpadNew[i]);
+ zin.push_back(ZpadNew[i]);
+ qin.push_back(QpadNew[i]);
+ npoint_temp++;
}
- }
-
- if(xout.size()>0)lenght_temp=sqrt(pow(zoutprime[indexin]-zoutprime[indexout],2)+pow(youtprime[indexin]-youtprime[indexout],2)+pow(xoutprime[indexin]-xoutprime[indexout],2));
-
- grxz.push_back(*grxztmp);
- gryz.push_back(*gryztmp);
- grxztmp->Set(0);
- gryztmp->Set(0);
-
- chargeTot.push_back(charge_temp);
- lenght.push_back(lenght_temp);
- } // end of if(xout.size()>10 &&)
-
- //X_mm->Fill();
-
- xin.clear();
- yin.clear();
- zin.clear();
- qin.clear();
- xout.clear();
- yout.clear();
- zout.clear();
- xoutprime.clear();
- youtprime.clear();
- zoutprime.clear();
- qout.clear();
- npoint_temp=0;
- ringsum=0;// angle between 1st track and z axis in 3D in degrees
- zmax=0.;
- for(int ko=0; ko<18; ko++) ringtouch[ko] = 0;
-
- } // if(xin.size()>0 && (( cluster_temp .... )
-
- cluster_temp = clusternbr[i];
-
- if(!(clusterpads[i]>=10 && clusterringbool[i]==1 && ZpadNew[i]>-100 && ZpadNew[i]<=310)) continue;
- else
- {
- xin.push_back(XpadNew[i]);
- yin.push_back(YpadNew[i]);
- zin.push_back(ZpadNew[i]);
- qin.push_back(QpadNew[i]);
- npoint_temp++;
- }
- }//end of loop on padsNews
-
-
- /* //------------------------------------------------------- */
- /* // STEP 3.2: Fitting the filtered tracks in 3D */
- /* // (weight by charge, TMinuit) */
- /* //------------------------------------------------------- */
-
-
- if(trackNbr_FINAL>= 1){ // !!! Just for 1 or more tracks!
- /* if(trackNbr_FINAL== 2){ */
-
- //////////Minimization in 3D to reconstruct track lines
- allevt_2pfiltered++;
-
- for(int itr= 0 ; itr<trackNbr_FINAL; itr++) {
-
-
- pStart[0]=0; pStart[2]=0; pStart[1]=1; pStart[3]=3;
- min = new TMinuit(4);
- min->SetPrintLevel(-1);
- arglist[0] = 3;
-
- Tracking_functions->FindStart(pStart,chi,fitStatus, &grxz.at(itr), &gryz.at(itr));
-
- NclusterFit = itr+1;
- current_phy=this;
- min->SetFCN(SumDistance);
-
-
- // Set starting values and step sizes for parameters
- min->mnparm(0,"x0",pStart[0],0.1,-500,500,iflag);
- min->mnparm(1,"Ax",pStart[1],0.1,-10,10,iflag);
- min->mnparm(2,"y0",pStart[2],0.1,-500,500,iflag);
- min->mnparm(3,"Ay",pStart[3],0.1,-10,10,iflag);
- arglist[0] = 200; // number of function calls
- arglist[1] = 0.000001; // tolerance
-
- min->mnexcm("MIGRAD",arglist,2,iflag); // minimization with MIGRAD
-
- min->mnstat(amin,edm,errdef,nvpar,nparx,iflag); //returns current status of the minimization
-
-
- // get fit parameters
-
- for(int i = 0; i <4; i++) min->GetParameter(i,parFit_temp[i],err_temp[i]);
-
- /* if( (parFit_temp[0] >-499 && parFit_temp[0]<499) && (parFit_temp[2] >-499 && parFit_temp[2]<499)) { */
- parFit1.push_back(parFit_temp[0]);
- parFit2.push_back(parFit_temp[1]);
- parFit3.push_back(parFit_temp[2]);
- parFit4.push_back(parFit_temp[3]);
- /* } */
- delete min;
- }
-
- if(trackNbr_FINAL==2){
-
- double ParTrack1[4];
- double ParTrack2[4];
-
- ParTrack1[0] = parFit1[0];
- ParTrack1[1] = parFit2[0];
- ParTrack1[2] = parFit3[0];
- ParTrack1[3] = parFit4[0];
-
- ParTrack2[0] = parFit1[1];
- ParTrack2[1] = parFit2[1];
- ParTrack2[2] = parFit3[1];
- ParTrack2[3] = parFit4[1];
-
- Dist_min=0, Theta_tr1=0, Theta_tr2=0;
- Tracking_functions->vertex(ParTrack1, ParTrack2, xv, yv, zv, Dist_min, Theta_tr1, Theta_tr2);
-
- /* if(abs(xv)<20 && abs(yv)<20){ */
- Xvertex.push_back(xv);
- Yvertex.push_back(yv);
- Zvertex.push_back(zv-Z_Shift);
-
- /* if (zv <-50){ */
- /* cout << ParTrack1[0]<< " " << ParTrack1[1]<< " " << ParTrack1[2]<<" " << ParTrack1[3]<<endl; */
- /* cout << xv<< " " << yv<< " " << zv<<endl; */
- /* } */
-
-
- if(Theta_tr1>Theta_tr2){
- Theta1.push_back(Theta_tr1);
- Theta2.push_back(Theta_tr2);
- }
- else if(Theta_tr1<Theta_tr2){
- Theta1.push_back(Theta_tr2);
- Theta2.push_back(Theta_tr1);
- }
+ }//end of PadNews
- sumTheta.push_back(Theta_tr1+Theta_tr2);
- /* sumTheta.push_back(156.9-(Theta_tr1+Theta_tr2)); */
- /* } */
- /* else{ */
-
- /* Xvertex.push_back(-99); */
- /* Yvertex.push_back(-99); */
- /* Zvertex.push_back(-99); */
- /* sumTheta.push_back(1000); */
- /* } */
-
- }
- else
- {
- Xvertex.push_back(-1000);
- Yvertex.push_back(-1000);
- Zvertex.push_back(-1000);
- sumTheta.push_back(1000);
- }
+ /* //------------------------------------------------------- */
+ /* // STEP 3.2: Fitting the filtered tracks in 3D */
+ /* // (weight by charge, TMinuit) */
+ /* //------------------------------------------------------- */
- }// end if trackNbr_FINAL>=1
- else
- {
- Xvertex.push_back(-1000);
- Yvertex.push_back(-1000);
- Zvertex.push_back(-1000);
- sumTheta.push_back(1000);
- }
-
- } // end loop if 0<trackNbr < 5
- else
- {
- Xvertex.push_back(-1000);
- Yvertex.push_back(-1000);
- Zvertex.push_back(-1000);
- sumTheta.push_back(1000);
- }
- // instantiate CalibrationManager
- static CalibrationManager* Cal = CalibrationManager::getInstance();
-}
-///////////////////////////////////////////////////////////////////////////
-void TMinosPhysics::ReadAnalysisConfig() {
- bool ReadingStatus = false;
+ /* if(trackNbr_FINAL>= 1){ // !!! Just for 1 or more tracks! */
+ if(trackNbr_FINAL== 2){
- // path to file
- string FileName = "./configs/ConfigMinos.dat";
+ //////////Minimization in 2D to reconstruct track lines
+ allevt_2pfiltered++;
- // open analysis config file
- ifstream AnalysisConfigFile;
- AnalysisConfigFile.open(FileName.c_str());
+ for(int itr= 0 ; itr < trackNbr_FINAL; itr++) {
- if (!AnalysisConfigFile.is_open()) {
- cout << " No ConfigMinos.dat found: Default parameter loaded for Analayis " << FileName << endl;
- return;
- }
- cout << " Loading user parameter for Analysis from ConfigMinos.dat " << endl;
-
- // Save it in a TAsciiFile
- TAsciiFile* asciiConfig = RootOutput::getInstance()->GetAsciiFileAnalysisConfig();
- asciiConfig->AppendLine("%%% ConfigMinos.dat %%%");
- asciiConfig->Append(FileName.c_str());
- asciiConfig->AppendLine("");
- // read analysis config file
- string LineBuffer,DataBuffer,whatToDo;
- while (!AnalysisConfigFile.eof()) {
- // Pick-up next line
- getline(AnalysisConfigFile, LineBuffer);
-
- // search for "header"
- string name = "ConfigMinos";
- if (LineBuffer.compare(0, name.length(), name) == 0)
- ReadingStatus = true;
-
- // loop on tokens and data
- while (ReadingStatus ) {
- whatToDo="";
- AnalysisConfigFile >> whatToDo;
-
- // Search for comment symbol (%)
- if (whatToDo.compare(0, 1, "%") == 0) {
- AnalysisConfigFile.ignore(numeric_limits<streamsize>::max(), '\n' );
- }
+ pStart[0]=0; pStart[2]=0; pStart[1]=1; pStart[3]=3;
+ min = new TMinuit(4);
+ min->SetPrintLevel(-1);
+ arglist[0] = 3;
- else if (whatToDo=="E_RAW_THRESHOLD") {
- AnalysisConfigFile >> DataBuffer;
- m_E_RAW_Threshold = atof(DataBuffer.c_str());
- cout << whatToDo << " " << m_E_RAW_Threshold << endl;
- }
+ Tracking_functions->FindStart(pStart,chi,fitStatus, &grxz.at(itr), &gryz.at(itr));
- else if (whatToDo=="E_THRESHOLD") {
- AnalysisConfigFile >> DataBuffer;
- m_E_Threshold = atof(DataBuffer.c_str());
- cout << whatToDo << " " << m_E_Threshold << endl;
- }
+ NclusterFit = itr+1;
+ current_phy=this;
+ min->SetFCN(SumDistance);
- else {
- ReadingStatus = false;
- }
- }
- }
+ // Set starting values and step sizes for parameters
+ min->mnparm(0,"x0",pStart[0],0.1,-500,500,iflag);
+ min->mnparm(1,"Ax",pStart[1],0.1,-10,10,iflag);
+ min->mnparm(2,"y0",pStart[2],0.1,-500,500,iflag);
+ min->mnparm(3,"Ay",pStart[3],0.1,-10,10,iflag);
+ arglist[0] = 200; // number of function calls
+ arglist[1] = 0.000001; // tolerance
-}
+ min->mnexcm("MIGRAD",arglist,2,iflag); // minimization with MIGRAD
-double TMinosPhysics::distance2(double x,double y,double z, double *p) {
- // distance line point is D= | (xp-x0) cross ux |
- // where ux is direction of line and x0 is a point in the line (like t = 0)
- ROOT::Math::XYZVector xp(x,y,z); //point of the track
- ROOT::Math:: XYZVector x0(p[0], p[2], 0. );
- ROOT::Math::XYZVector x1(p[0] + p[1], p[2] + p[3], 1. ); //line
- ROOT::Math::XYZVector u = (x1-x0).Unit();
- double d2 = ((xp-x0).Cross(u)) .Mag2();
- return d2;
-}
+ min->mnstat(amin,edm,errdef,nvpar,nparx,iflag); //returns current status of the minimization
-void TMinosPhysics::SumDistance(int &, double *, double & sum, double * par, int) {
- TMinosPhysics* phy = current_phy;
- int nused=0;
- double qtot=0;
- sum = 0;
-
- for(int i=0; i < phy->data_result.GetEntriesFast(); i++)
- {
- phy->minosdata_result = (TMinosResult*)phy->data_result.At(i);
- if(phy->minosdata_result->n_Cluster==NclusterFit)
- {
- float x=phy->minosdata_result->x_mm;
- float y=phy->minosdata_result->y_mm;
- float z=phy->minosdata_result->z_mm;
- float q=phy->minosdata_result->Chargemax;
- //if(nused<2)cout<<minosdata_result->n_Cluster<<" "<<x<<" "<<y<<" "<<z<<" "<<q<<endl;
- /* double d = Tracking_functions->distance2(x, y, z, par); */
- double d = TMinosPhysics::distance2(x, y, z, par);
- sum += d*q;
- nused++;
- qtot+=q;
- }
- }
- //sum/=nused;
- sum/=qtot;
- return;
-}
+ // get fit parameters
+ for(int i = 0; i <4; i++) min->GetParameter(i,parFit_temp[i],err_temp[i]);
-///////////////////////////////////////////////////////////////////////////
-void TMinosPhysics::Clear() {
-
- Xpad.clear();
- Ypad.clear();
- Qpad.clear();
- XpadNew.clear();
- YpadNew.clear();
- ZpadNew.clear();
- QpadNew.clear();
-
- clusterringboolTemp.clear();
- clusterringbool.clear();
- clusternbr.clear();
- clusterpads.clear();
- hfit->Reset();
-
- xin.clear();
- yin.clear();
- zin.clear();
- qin.clear();
- xout.clear();
- yout.clear();
- zout.clear();
- xoutprime.clear();
- youtprime.clear();
- zoutprime.clear();
- qout.clear();
-
- trackclusternbr.clear();
- tracknbr.clear();
- TOTxoutprime.clear();
- TOTyoutprime.clear();
- TOTzoutprime.clear();
- TOTqout.clear();
-
- lenght.clear();
- chargeTot.clear();
- parFit1.clear();
- parFit2.clear();
- parFit3.clear();
- parFit4.clear();
-
- Xvertex.clear();
- Yvertex.clear();
- Zvertex.clear();
- sumTheta.clear();
-
- grxz.clear();
- gryz.clear();
-
- Theta1.clear();
- Theta2.clear();
-
- hfit->Reset();
-
-
- filled=0;
- indexfill = 0;
- ChargeBin = 0.;
- maxCharge = 0.;
- Iteration=0;
- filter_result=0;
- fit2DStatus=0;
- trackNbr=0;
- trackNbr_FINAL=0;
- x_mm = 0.; y_mm = 0.; z_mm = 0.; q_pad = 0.; t_pad = 0.;
- array_final=0;
- ringsum=0;
- zmax=0.;
+ /* if( (parFit_temp[0] >-499 && parFit_temp[0]<499) && (parFit_temp[2] >-499 && parFit_temp[2]<499)) { */
+ /* parFit1[itr] = push_back(parFit_temp[0]); */
+ /* parFit2.push_back(parFit_temp[1]); */
+ /* parFit3.push_back(parFit_temp[2]); */
+ /* parFit4.push_back(parFit_temp[3]); */
+ parFit1[itr]=parFit_temp[0];
+ parFit2[itr]=parFit_temp[1];
+ parFit3[itr]=parFit_temp[2];
+ parFit4[itr]=parFit_temp[3];
-}
+ /* } */
+ delete min;
+ }
-///////////////////////////////////////////////////////////////////////////
-void TMinosPhysics::ReadConfiguration(NPL::InputParser parser) {
- /* vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("Minos"); */
- /* if(NPOptionManager::getInstance()->GetVerboseLevel()) */
- /* cout << "//// " << blocks.size() << " detectors found " << endl; */
-
- /* vector<string> cart = {"POS","Shape"}; */
- /* vector<string> sphe = {"R","THETA","POS","Phi","TargetLength"}; */
-
- /* for(unsigned int i = 0 ; i < blocks.size() ; i++){ */
- /* if(blocks[i]->HasTokenList(cart)){ */
- /* if(NPOptionManager::getInstance()->GetVerboseLevel()) */
- /* cout << endl << "//// Minos " << i+1 << endl; */
-
- /* TVector3 Pos = blocks[i]->GetTVector3("POS","mm"); */
- /* string Shape = blocks[i]->GetString("Shape"); */
- /* AddDetector(Pos,Shape); */
- /* } */
- /* else if(blocks[i]->HasTokenList(sphe)){ */
- /* if(NPOptionManager::getInstance()->GetVerboseLevel()) */
- /* cout << endl << "//// Minos " << i+1 << endl; */
- /* double R = blocks[i]->GetDouble("R","mm"); */
- /* double Theta = blocks[i]->GetDouble("THETA","deg"); */
- /* TVector3 POS = blocks[i]->GetTVector3("POS","cm"); */
- /* double Phi = blocks[i]->GetDouble("PHI","deg"); */
- /* double TargetLenght = blocks[i]->GetDouble("TargetLength","mm"); */
- /* AddDetector(R,Theta,POS,Phi,TargetLenght); */
- /* } */
- /* else{ */
- /* cout << "ERROR: check your input file formatting " << endl; */
- /* exit(1); */
- /* } */
- /* } */
-}
+ if(trackNbr_FINAL==2){
+ /* if(trackNbr_FINAL>=0){ */
+ double ParTrack1[4];
+ double ParTrack2[4];
-///////////////////////////////////////////////////////////////////////////
-void TMinosPhysics::AddParameterToCalibrationManager() {
- CalibrationManager* Cal = CalibrationManager::getInstance();
- for (int i = 0; i < m_NumberOfDetectors; ++i) {
- Cal->AddParameter("Minos", "D"+ NPL::itoa(i+1)+"_ENERGY","Minos_D"+ NPL::itoa(i+1)+"_ENERGY");
- Cal->AddParameter("Minos", "D"+ NPL::itoa(i+1)+"_TIME","Minos_D"+ NPL::itoa(i+1)+"_TIME");
- }
-}
+ ParTrack1[0] = parFit1[0];
+ ParTrack1[1] = parFit2[0];
+ ParTrack1[2] = parFit3[0];
+ ParTrack1[3] = parFit4[0];
-///////////////////////////////////////////////////////////////////////////
-void TMinosPhysics::InitializeRootInputRaw() {
- TChain* inputChain = RootInput::getInstance()->GetChain();
- inputChain->SetBranchStatus("Minos", true );
- inputChain->SetBranchAddress("Minos", &m_EventData );
- if( strlen(inputChain->GetTree()->GetName())==13) SimulationBool= true; // Test if TreeName==SimulatedTree
-}
+ ParTrack2[0] = parFit1[1];
+ ParTrack2[1] = parFit2[1];
+ ParTrack2[2] = parFit3[1];
+ ParTrack2[3] = parFit4[1];
-///////////////////////////////////////////////////////////////////////////
-void TMinosPhysics::InitializeRootInputPhysics() {
- TChain* inputChain = RootInput::getInstance()->GetChain();
- inputChain->SetBranchAddress("Minos", &m_EventPhysics);
-}
+ Dist_min=0, Theta_tr1=0, Theta_tr2=0;
-///////////////////////////////////////////////////////////////////////////
-void TMinosPhysics::InitializeRootOutput() {
- TTree* outputTree = RootOutput::getInstance()->GetTree();
- outputTree->Branch("Minos", "TMinosPhysics", &m_EventPhysics);
-}
+ Tracking_functions->vertex(ParTrack1, ParTrack2, xv, yv, zv, Dist_min, Theta_tr1, Theta_tr2, Phi1, Phi2);
+
+ Xvertex=xv;
+ Yvertex=yv;
+ Zvertex = zv;
+ /* Zvertex = zv+Z_Shift; */
+
+ Dmin = Dist_min;
+ Theta_1 = Theta_tr1;
+ Theta_2 = Theta_tr2 ;
+
+ Theta1 = Theta_tr1;
+ Theta2 = Theta_tr2;
-////////////////////////////////////////////////////////////////////////////////
-// Construct Method to be pass to the DetectorFactory //
-////////////////////////////////////////////////////////////////////////////////
-NPL::VDetector* TMinosPhysics::Construct() {
- return (NPL::VDetector*) new TMinosPhysics();
-}
+ sumTheta = Theta_tr1+Theta_tr2;
+ }
-////////////////////////////////////////////////////////////////////////////////
-// Registering the construct method to the factory //
-////////////////////////////////////////////////////////////////////////////////
-extern "C"{
-class proxy_Minos{
- public:
- proxy_Minos(){
- NPL::DetectorFactory::getInstance()->AddToken("Minos","Minos");
- NPL::DetectorFactory::getInstance()->AddDetector("Minos",TMinosPhysics::Construct);
- }
-};
-
-proxy_Minos p_Minos;
-}
+ }// end if trackNbr_FINAL>=1
+
+ } // end loop if 0<trackNbr < 5
+
+ // instantiate CalibrationManager
+ static CalibrationManager* Cal = CalibrationManager::getInstance();
+ }
+
+ ///////////////////////////////////////////////////////////////////////////
+ void TMinosPhysics::ReadAnalysisConfig() {
+ bool ReadingStatus = false;
+
+ // path to file
+ string FileName = "./configs/ConfigMinos.dat";
+
+ // open analysis config file
+ ifstream AnalysisConfigFile;
+ AnalysisConfigFile.open(FileName.c_str());
+
+ if (!AnalysisConfigFile.is_open()) {
+ cout << " No ConfigMinos.dat found: Default parameter loaded for Analayis " << FileName << endl;
+ return;
+ }
+ cout << " Loading user parameter for Analysis from ConfigMinos.dat " << endl;
+
+ // Save it in a TAsciiFile
+ TAsciiFile* asciiConfig = RootOutput::getInstance()->GetAsciiFileAnalysisConfig();
+ asciiConfig->AppendLine("%%% ConfigMinos.dat %%%");
+ asciiConfig->Append(FileName.c_str());
+ asciiConfig->AppendLine("");
+ // read analysis config file
+ string LineBuffer,DataBuffer,whatToDo;
+ while (!AnalysisConfigFile.eof()) {
+ // Pick-up next line
+ getline(AnalysisConfigFile, LineBuffer);
+
+ // search for "header"
+ string name = "ConfigMinos";
+ if (LineBuffer.compare(0, name.length(), name) == 0)
+ ReadingStatus = true;
+
+ // loop on tokens and data
+ while (ReadingStatus ) {
+ whatToDo="";
+ AnalysisConfigFile >> whatToDo;
+
+ // Search for comment symbol (%)
+ if (whatToDo.compare(0, 1, "%") == 0) {
+ AnalysisConfigFile.ignore(numeric_limits<streamsize>::max(), '\n' );
+ }
+
+ else if (whatToDo=="E_RAW_THRESHOLD") {
+ AnalysisConfigFile >> DataBuffer;
+ m_E_RAW_Threshold = atof(DataBuffer.c_str());
+ cout << whatToDo << " " << m_E_RAW_Threshold << endl;
+ }
+
+ else if (whatToDo=="E_THRESHOLD") {
+ AnalysisConfigFile >> DataBuffer;
+ m_E_Threshold = atof(DataBuffer.c_str());
+ cout << whatToDo << " " << m_E_Threshold << endl;
+ }
+
+ else {
+ ReadingStatus = false;
+ }
+ }
+ }
+
+ }
+
+ double TMinosPhysics::distance2(double x,double y,double z, double *p) {
+ // distance line point is D= | (xp-x0) cross ux |
+ // where ux is direction of line and x0 is a point in the line (like t = 0)
+ ROOT::Math::XYZVector xp(x,y,z); //point of the track
+ ROOT::Math:: XYZVector x0(p[0], p[2], 0. );
+ ROOT::Math::XYZVector x1(p[0] + p[1], p[2] + p[3], 1. ); //line
+ ROOT::Math::XYZVector u = (x1-x0).Unit();
+ double d2 = ((xp-x0).Cross(u)) .Mag2();
+ return d2;
+ }
+
+ void TMinosPhysics::SumDistance(int &, double *, double & sum, double * par, int) {
+ TMinosPhysics* phy = current_phy;
+ int nused=0;
+ double qtot=0;
+ sum = 0;
+
+ for(int i=0; i < phy->data_result.GetEntriesFast(); i++)
+ {
+ phy->minosdata_result = (TMinosResult*)phy->data_result.At(i);
+ if(phy->minosdata_result->n_Cluster==NclusterFit)
+ {
+ float x=phy->minosdata_result->x_mm;
+ float y=phy->minosdata_result->y_mm;
+ float z=phy->minosdata_result->z_mm;
+ float q=phy->minosdata_result->Chargemax;
+ //if(nused<2)cout<<minosdata_result->n_Cluster<<" "<<x<<" "<<y<<" "<<z<<" "<<q<<endl;
+ /* double d = Tracking_functions->distance2(x, y, z, par); */
+ double d = TMinosPhysics::distance2(x, y, z, par);
+ sum += d*q;
+ nused++;
+ qtot+=q;
+ }
+ }
+ //sum/=nused;
+ sum/=qtot;
+ return;
+ }
+
+ ///////////////////////////////////////////////////////////////////////////
+ void TMinosPhysics::Clear() {
+
+ Xpad.clear();
+ Ypad.clear();
+ Qpad.clear();
+ XpadNew.clear();
+ YpadNew.clear();
+ ZpadNew.clear();
+ QpadNew.clear();
+
+ clusterringboolTemp.clear();
+ clusterringbool.clear();
+ clusternbr.clear();
+ clusterpads.clear();
+ hfit->Reset();
+
+ xin.clear();
+ yin.clear();
+ zin.clear();
+ qin.clear();
+ xout.clear();
+ yout.clear();
+ zout.clear();
+ xoutprime.clear();
+ youtprime.clear();
+ zoutprime.clear();
+ qout.clear();
+
+ trackclusternbr.clear();
+ tracknbr.clear();
+ TOTxoutprime.clear();
+ TOTyoutprime.clear();
+ TOTzoutprime.clear();
+ TOTqout.clear();
+
+ lenght.clear();
+ chargeTot.clear();
+ parFit1.clear();
+ parFit2.clear();
+ parFit3.clear();
+ parFit4.clear();
+
+ grxz.clear();
+ gryz.clear();
+
+ hfit->Reset();
+
+
+ filled=0;
+ indexfill = 0;
+ ChargeBin = 0.;
+ maxCharge = 0.;
+ Iteration=0;
+ filter_result=0;
+ fit2DStatus=0;
+ trackNbr=0;
+ trackNbr_FINAL=0;
+ x_mm = 0.; y_mm = 0.; z_mm = 0.; q_pad = 0.; t_pad = 0.;
+ array_final=0;
+ ringsum=0;
+ zmax=0.;
+
+
+ }
+
+ ///////////////////////////////////////////////////////////////////////////
+ void TMinosPhysics::ReadConfiguration(NPL::InputParser parser) {
+ }
+
+
+ ///////////////////////////////////////////////////////////////////////////
+ void TMinosPhysics::AddParameterToCalibrationManager() {
+ CalibrationManager* Cal = CalibrationManager::getInstance();
+ for (int i = 0; i < m_NumberOfDetectors; ++i) {
+ Cal->AddParameter("Minos", "D"+ NPL::itoa(i+1)+"_ENERGY","Minos_D"+ NPL::itoa(i+1)+"_ENERGY");
+ Cal->AddParameter("Minos", "D"+ NPL::itoa(i+1)+"_TIME","Minos_D"+ NPL::itoa(i+1)+"_TIME");
+ }
+ }
+
+ ///////////////////////////////////////////////////////////////////////////
+ void TMinosPhysics::InitializeRootInputRaw() {
+ TChain* inputChain = RootInput::getInstance()->GetChain();
+ inputChain->SetBranchStatus("Minos", true );
+ inputChain->SetBranchAddress("Minos", &m_EventData );
+ if( strlen(inputChain->GetTree()->GetName())==13) SimulationBool= true; // Test if TreeName==SimulatedTree
+ }
+
+ ///////////////////////////////////////////////////////////////////////////
+ void TMinosPhysics::InitializeRootInputPhysics() {
+ TChain* inputChain = RootInput::getInstance()->GetChain();
+ inputChain->SetBranchAddress("Minos", &m_EventPhysics);
+ }
+
+ ///////////////////////////////////////////////////////////////////////////
+ void TMinosPhysics::InitializeRootOutput() {
+ TTree* outputTree = RootOutput::getInstance()->GetTree();
+ outputTree->Branch("Minos", "TMinosPhysics", &m_EventPhysics);
+ }
+
+ ////////////////////////////////////////////////////////////////////////////////
+ // Construct Method to be pass to the DetectorFactory //
+ ////////////////////////////////////////////////////////////////////////////////
+ NPL::VDetector* TMinosPhysics::Construct() {
+ return (NPL::VDetector*) new TMinosPhysics();
+ }
+
+ ////////////////////////////////////////////////////////////////////////////////
+ // Registering the construct method to the factory //
+ ////////////////////////////////////////////////////////////////////////////////
+ extern "C"{
+ class proxy_Minos{
+ public:
+ proxy_Minos(){
+ NPL::DetectorFactory::getInstance()->AddToken("Minos","Minos");
+ NPL::DetectorFactory::getInstance()->AddDetector("Minos",TMinosPhysics::Construct);
+ }
+ };
+
+ proxy_Minos p_Minos;
+ }
diff --git a/NPLib/Detectors/Minos/TMinosPhysics.h b/NPLib/Detectors/Minos/TMinosPhysics.h
index 23d9d57e5..2ac1d0cb9 100644
--- a/NPLib/Detectors/Minos/TMinosPhysics.h
+++ b/NPLib/Detectors/Minos/TMinosPhysics.h
@@ -63,8 +63,6 @@ class TMinosPhysics : public TObject, public NPL::VDetector {
void Clear(const Option_t*) {};
- /* static void SumDistance(int &, double *, double & sum, double * par, int); */
-
//////////////////////////////////////////////////////////////
// data obtained after BuildPhysicalEvent() and stored in
// output ROOT file
@@ -116,12 +114,10 @@ class TMinosPhysics : public TObject, public NPL::VDetector {
NPL::Tracking* Tracking_functions; //!
- /* int NclusterFit; //! */
// fit function for the Q(t) signals at each pad
static double conv_fit(double *x, double *p);
static double distance2(double x,double y,double z, double *p);
static void SumDistance(int &, double *, double & sum, double * par, int);
-
// remove bad channels, calibrate the data and apply thresholds
void PreTreat();
@@ -164,16 +160,22 @@ TMinosResult *minosdata_result;//!
vector<double> GetPad_X() {return X_Pad;}
vector<double> GetPad_Y() {return Y_Pad;}
vector<double> GetPad_Z() {return Z_Pad;}
+ vector<double> GetPad_T() {return T_Pad;}
vector<double> GetParFit1() {return parFit1;}
vector<double> GetParFit2() {return parFit2;}
vector<double> GetParFit3() {return parFit3;}
vector<double> GetParFit4() {return parFit4;}
- double GetVertexZ() {return Zvertex[0];}
+ double GetVertexZ() {return Zvertex;}
+ double GetVertexX() {return Xvertex;}
+ double GetVertexY() {return Yvertex;}
- /* TClonesArray fitdata; */
-
- /* fitdata.SetClass("TMinosClust"); */
-
+ double GetTheta1() {return Theta1 ;}
+ double GetTheta2() {return Theta2 ;}
+ double GetTheta_1() {return Theta_1 ;}
+ double GetTheta_2() {return Theta_2 ;}
+
+ double GetPhi1(){return Phi1;}
+ double GetPhi2(){return Phi2;}
// parameters used in the analysis
private:
@@ -256,13 +258,19 @@ TMinosResult *minosdata_result;//!
vector<double> TOTzoutprime;
vector<double> TOTqout;
- vector<double> Xvertex;
- vector<double> Yvertex;
- vector<double> Zvertex;
+ double Xvertex;
+ double Yvertex;
+ double Zvertex;
+ double Dmin;
- vector<double> sumTheta;
- vector<double> Theta1;
- vector<double> Theta2;
+ double sumTheta;
+ double Theta1;
+ double Theta2;
+ double Theta_1;
+ double Theta_2;
+
+ double Phi1;
+ double Phi2;
vector<int> trackclusternbr;//!
vector<int> tracknbr;//!
@@ -302,12 +310,11 @@ TMinosResult *minosdata_result;//!
vector<double> errFit3_local;
vector<double> errFit4_local;
-
double xv;//!
double yv;//!
double zv;//!
- double Dist_min;//!
+ double Dist_min;
double Theta_tr1;//!
double Theta_tr2;//!
// number of detectors
diff --git a/NPLib/TrackReconstruction/Tracking.cxx b/NPLib/TrackReconstruction/Tracking.cxx
index 02bda352a..3edd675d5 100644
--- a/NPLib/TrackReconstruction/Tracking.cxx
+++ b/NPLib/TrackReconstruction/Tracking.cxx
@@ -48,427 +48,446 @@ NPL::Tracking::~Tracking()
/* } */
int NPL::Tracking::Hough_modified(vector<double> *x,vector<double> *y,vector<double> *q, vector<double> *x_out,vector<double> *y_out, vector<double> *q_out, vector<int> *ringbool) {
- double bint1=2.;
- double bint2=2.;
- int maxt = 360.;
- int mint = 0.;
- int nt1=(maxt-mint)/bint1;
- int nt2=(maxt-mint)/bint1;
-
- double PI = TMath::Pi();
- double Rint = 45.2;
- double Rext = 45.2 + 18*2.1;
- int filter_result = 0;
-
- TH2F *hp_xy = new TH2F("hp_xy","hp_xy",nt1,mint,maxt,nt2,mint,maxt);
- TH2F *hpDiag_xy = new TH2F("hpDiag_xy","hpDiag_xy",nt1,mint,maxt,nt2,mint,maxt);
-
- double max_xy;
-
- vector<double> xTemp, yTemp, qTemp;
-
- double theta1, theta2, xt1, yt1, xt2, yt2;
- double line0=0., line1=0.;
- double delta=0., AA=0., BB=0., CC=0.;
- double maxtheta1=0., maxtheta2=0., xmax1=0., ymax1=0., xmax2=0., ymax2=0.;
- double par0=0., par1=0.;
- double r_mm=0.;
- int ringsum=0;
- bool maxfound = false;
-
- for(unsigned int i=0;i<x->size();i++){
- xTemp.push_back(x->at(i));
- yTemp.push_back(y->at(i));
- qTemp.push_back(q->at(i));
-
+ double bint1=2.;
+ double bint2=2.;
+ int maxt = 360.;
+ int mint = 0.;
+ int nt1=(maxt-mint)/bint1;
+ int nt2=(maxt-mint)/bint1;
+
+ double PI = TMath::Pi();
+ double Rint = 45.2;
+ double Rext = 45.2 + 18*2.1;
+ int filter_result = 0;
+
+ TH2F *hp_xy = new TH2F("hp_xy","hp_xy",nt1,mint,maxt,nt2,mint,maxt);
+ TH2F *hpDiag_xy = new TH2F("hpDiag_xy","hpDiag_xy",nt1,mint,maxt,nt2,mint,maxt);
+
+ double max_xy;
+
+ vector<double> xTemp, yTemp, qTemp;
+
+ double theta1, theta2, xt1, yt1, xt2, yt2;
+ double line0=0., line1=0.;
+ double delta=0., AA=0., BB=0., CC=0.;
+ double maxtheta1=0., maxtheta2=0., xmax1=0., ymax1=0., xmax2=0., ymax2=0.;
+ double par0=0., par1=0.;
+ double r_mm=0.;
+ int ringsum=0;
+ bool maxfound = false;
+
+ for(unsigned int i=0;i<x->size();i++){
+ xTemp.push_back(x->at(i));
+ yTemp.push_back(y->at(i));
+ qTemp.push_back(q->at(i));
+
//Loop of indices
- for(int j=0;j<nt1;j++){
-
+ for(int j=0;j<nt1;j++){
+
theta1 = (j+0.5)*bint1 + mint;
- xt1 = Rint * TMath::Cos(theta1*PI/180.);
- yt1 = Rint * TMath::Sin(theta1*PI/180.);
- line1 = (yt1 - y->at(i))/(xt1 - x->at(i));
- line0 = yt1 - xt1 * line1;
- AA = 1 + line1*line1;
- BB = 2*line0*line1;
- CC = line0*line0 - Rext*Rext;
-
- delta = BB*BB - 4*AA*CC;
-
+ xt1 = Rint * TMath::Cos(theta1*PI/180.);
+ yt1 = Rint * TMath::Sin(theta1*PI/180.);
+ line1 = (yt1 - y->at(i))/(xt1 - x->at(i));
+ line0 = yt1 - xt1 * line1;
+ AA = 1 + line1*line1;
+ BB = 2*line0*line1;
+ /* CC = line0*line0 - Rext; // Warning likely Error : Rext and not Rext^2 (Cyril) */
+ CC = line0*line0 - Rext*Rext; // Warning likely Error : Rext and not Rext^2 (Cyril)
+
+ delta = BB*BB - 4*AA*CC;
+
if(delta>=0){
- xt2 = (-BB - sqrt(delta))/(2*AA);
- yt2 = line0 + line1*xt2;
-
+ xt2 = (-BB - sqrt(delta))/(2*AA);
+ yt2 = line0 + line1*xt2;
if(xt2<=0) theta2= 180 - asin(yt2/Rext)*180/PI;
-
else if(xt2>0){
- if(yt2>0) theta2= asin(yt2/Rext)*180/PI;
+ if(yt2>0) theta2= asin(yt2/Rext)*180/PI;
else if(yt2<=0) theta2= 360 + asin(yt2/Rext)*180/PI;
- }
-
- if( (xt1*x->at(i) + yt1*y->at(i))>=0 && (xt2*x->at(i) + yt2*y->at(i))>=0 && (xt1*xt2+yt1*yt2)>=0){
- hp_xy->Fill(theta1,theta2);
- if(abs(theta1-theta2)<=10) hpDiag_xy->Fill(theta1,theta2);
- }
- else{
- if(delta!=0){
- xt2 = (-BB + sqrt(delta))/(2*AA);
- yt2 = line0 + line1*xt2;
-
+ }
+ if( (xt1*x->at(i) + yt1*y->at(i))>=0 && (xt2*x->at(i) + yt2*y->at(i))>=0 && (xt1*xt2+yt1*yt2)>=0){
+ hp_xy->Fill(theta1,theta2);
+ if(abs(theta1-theta2)<=10) hpDiag_xy->Fill(theta1,theta2);
+ }
+ else{
+ if(delta!=0){
+ xt2 = (-BB + sqrt(delta))/(2*AA);
+ yt2 = line0 + line1*xt2;
if(xt2<=0) theta2= 180 - asin(yt2/Rext)*180/PI;
-
else if(xt2>0){
if(yt2>0) theta2= asin(yt2/Rext)*180/PI;
else if(yt2<=0) theta2= 360 + asin(yt2/Rext)*180/PI;
- }
-
+ }
if( (xt1*x->at(i) + yt1*y->at(i))>=0 && (xt2*x->at(i) + yt2*y->at(i))>=0 && (xt1*xt2+yt1*yt2)>=0){
- hp_xy->Fill(theta1,theta2);
+ hp_xy->Fill(theta1,theta2);
if(abs(theta1-theta2)<=10) hpDiag_xy->Fill(theta1,theta2);
- }
- }
- }
- } // end if delta>=0
- } // end loop on indices
- } // end loop on pads
-
+ }
+ }
+ }
+ } // end if delta>=0
+ } // end loop on indices
+ } // end loop on pads
+
x->clear();
- y->clear();
- q->clear();
-
- if(hpDiag_xy->GetMaximum()>=10) max_xy = hpDiag_xy->GetMaximum();
-// cout << "Max taken in diag... withh value=" << max_xy << endl;
- else max_xy = hp_xy->GetMaximum();
-
- for(int ii=0; ii<nt1; ii++){
- if(maxfound ==true) break;
- for(int jj=0; jj<nt2; jj++){
- if(hp_xy->GetBinContent(ii+1, jj+1) == max_xy){
- maxtheta1 = (ii+0.5)*bint1 + mint;
- maxtheta2 = (jj+0.5)*bint2 + mint;
- maxfound = true;
- //cout << "xy: theta max are " << maxtheta1 << " , " << maxtheta2 << endl;
- }
- if(maxfound ==true) break;
- }
- }
-
- xmax1 = Rint * TMath::Cos(maxtheta1*PI/180.);
- ymax1 = Rint * TMath::Sin(maxtheta1*PI/180.);
- xmax2 = Rext * TMath::Cos(maxtheta2*PI/180.);
- ymax2 = Rext * TMath::Sin(maxtheta2*PI/180.);
-
- // xy PEAK
- par1 = (ymax2-ymax1)/(xmax2-xmax1);
- par0 = (ymax1 - xmax1*par1);
- //Selection of x,y points IN the maxmean+/-1 found in Obertelli transform of xy plane
- for(unsigned int i=0;i<xTemp.size();i++){
- if( (abs(par1*xTemp[i]-yTemp[i]+par0)/sqrt(1+par1*par1))<= 6 && ((xmax1*xTemp[i] + ymax1*yTemp[i]) >= 0) && ((xmax2*xTemp[i] + ymax2*yTemp[i]) >= 0) && ((xmax1*xmax2 + ymax1*ymax2) >= 0)){
-// couti<< "Taken points= " << xTemp[i] << " , " << yTemp[i] << " , " << zTemp[i] << endl;
-// hcnew_xy->Fill(xTemp[i],yTemp[i],qTemp[i]);
- x_out->push_back(xTemp[i]);
- y_out->push_back(yTemp[i]);
- q_out->push_back(qTemp[i]);
- filter_result++;
- /* r_mm = sqrt(xTemp[i]*xTemp[i]+yTemp[i]*yTemp[i]); */
+ y->clear();
+ q->clear();
+
+ if(hpDiag_xy->GetMaximum()>=10) max_xy = hpDiag_xy->GetMaximum();
+ // cout << "Max taken in diag... withh value=" << max_xy << endl;
+ else max_xy = hp_xy->GetMaximum();
+
+ for(int ii=0; ii<nt1; ii++){
+ if(maxfound ==true) break;
+ for(int jj=0; jj<nt2; jj++){
+ if(hp_xy->GetBinContent(ii+1, jj+1) == max_xy){
+ maxtheta1 = (ii+0.5)*bint1 + mint;
+ maxtheta2 = (jj+0.5)*bint2 + mint;
+ maxfound = true;
+ //cout << "xy: theta max are " << maxtheta1 << " , " << maxtheta2 << endl;
+ }
+ if(maxfound ==true) break;
+ }
+ }
+
+ xmax1 = Rint * TMath::Cos(maxtheta1*PI/180.);
+ ymax1 = Rint * TMath::Sin(maxtheta1*PI/180.);
+ xmax2 = Rext * TMath::Cos(maxtheta2*PI/180.);
+ ymax2 = Rext * TMath::Sin(maxtheta2*PI/180.);
+
+ // xy PEAK
+ par1 = (ymax2-ymax1)/(xmax2-xmax1);
+ par0 = (ymax1 - xmax1*par1);
+ //Selection of x,y points IN the maxmean+/-1 found in Obertelli transform of xy plane
+ for(unsigned int i=0;i<xTemp.size();i++){
+ if( (abs(par1*xTemp[i]-yTemp[i]+par0)/sqrt(1+par1*par1))<= 6 && ((xmax1*xTemp[i] + ymax1*yTemp[i]) >= 0) && ((xmax2*xTemp[i] + ymax2*yTemp[i]) >= 0) && ((xmax1*xmax2 + ymax1*ymax2) >= 0)){
+ // couti<< "Taken points= " << xTemp[i] << " , " << yTemp[i] << " , " << zTemp[i] << endl;
+ // hcnew_xy->Fill(xTemp[i],yTemp[i],qTemp[i]);
+ x_out->push_back(xTemp[i]);
+ y_out->push_back(yTemp[i]);
+ q_out->push_back(qTemp[i]);
+ filter_result++;
+ /* r_mm = sqrt(xTemp[i]*xTemp[i]+yTemp[i]*yTemp[i]); */
/* if(r_mm<(45.2+5*2.1)) ringsum++; */ // commented by Cyril
-
- }
- else{
- x->push_back(xTemp[i]);
- y->push_back(yTemp[i]);
- q->push_back(qTemp[i]);
-
- }
- }
-
+ }
+ else{
+ x->push_back(xTemp[i]);
+ y->push_back(yTemp[i]);
+ q->push_back(qTemp[i]);
+ }
+ }
/* cout << "ringsum = " << ringsum << endl; */
/* cout << "Filter_result = " << filter_result << endl; */
- for(int ip=0; ip<filter_result; ip++){
- /* if(ringsum>2) ringbool->push_back(1); */
- /* else ringbool->push_back(0); */
- ringbool->push_back(1);
- }
-
- delete hp_xy;
- delete hpDiag_xy;
-
- return filter_result;
-
+ for(int ip=0; ip<filter_result; ip++){
+ /* if(ringsum>2) ringbool->push_back(1); */
+ /* else ringbool->push_back(0); */
+ ringbool->push_back(1);
+ }
+ delete hp_xy;
+ delete hpDiag_xy;
+ return filter_result;
}
double NPL::Tracking::FitFunction(double *x, double *p) {
- double val=p[0]+p[1]*x[0];
- return(val);
+ double val=p[0]+p[1]*x[0];
+ return(val);
}
void NPL::Tracking::FindStart(double pStart[4], double chi[2], int fitStatus[2],TGraph *grxz, TGraph *gryz) {
- double par1D[2];
- TF1 *myfit1 = new TF1("myfit1","[0]+[1]*x", -100,500);
- myfit1->SetParameter(0,0);
- myfit1->SetParameter(1,10);
- fitStatus[0] =0;
- grxz->Fit(myfit1,"RQM");
- chi[0]=myfit1->GetChisquare();
- par1D[0]=myfit1->GetParameter(0);
- par1D[1]=myfit1->GetParameter(1);
- pStart[0]=par1D[0];
- pStart[1]=par1D[1];
- fitStatus[1] =0;
- gryz->Fit(myfit1,"RQM");
- chi[1]=myfit1->GetChisquare();
- par1D[0]=myfit1->GetParameter(0);
- par1D[1]=myfit1->GetParameter(1);
- pStart[2]=par1D[0];
- pStart[3]=par1D[1];
- //AC 07/12/14
- delete myfit1;
+ double par1D[2];
+ TF1 *myfit1 = new TF1("myfit1","[0]+[1]*x", -100,500);
+ myfit1->SetParameter(0,0);
+ myfit1->SetParameter(1,10);
+ fitStatus[0] =0;
+ grxz->Fit(myfit1,"RQM");
+ chi[0]=myfit1->GetChisquare();
+ par1D[0]=myfit1->GetParameter(0);
+ par1D[1]=myfit1->GetParameter(1);
+ pStart[0]=par1D[0];
+ pStart[1]=par1D[1];
+ fitStatus[1] =0;
+ gryz->Fit(myfit1,"RQM");
+ chi[1]=myfit1->GetChisquare();
+ par1D[0]=myfit1->GetParameter(0);
+ par1D[1]=myfit1->GetParameter(1);
+ pStart[2]=par1D[0];
+ pStart[3]=par1D[1];
+ //AC 07/12/14
+ delete myfit1;
}
// Calculation of the distance line-point
double NPL::Tracking::distance2(double x,double y,double z, double *p) {
- // distance line point is D= | (xp-x0) cross ux |
- // where ux is direction of line and x0 is a point in the line (like t = 0)
- ROOT::Math::XYZVector xp(x,y,z); //point of the track
- ROOT::Math:: XYZVector x0(p[0], p[2], 0. );
- ROOT::Math::XYZVector x1(p[0] + p[1], p[2] + p[3], 1. ); //line
- ROOT::Math::XYZVector u = (x1-x0).Unit();
- double d2 = ((xp-x0).Cross(u)) .Mag2();
- return d2;
+ // distance line point is D= | (xp-x0) cross ux |
+ // where ux is direction of line and x0 is a point in the line (like t = 0)
+ ROOT::Math::XYZVector xp(x,y,z); //point of the track
+ ROOT::Math:: XYZVector x0(p[0], p[2], 0. );
+ ROOT::Math::XYZVector x1(p[0] + p[1], p[2] + p[3], 1. ); //line
+ ROOT::Math::XYZVector u = (x1-x0).Unit();
+ double d2 = ((xp-x0).Cross(u)) .Mag2();
+ return d2;
}
//void Hough_3D(TCanvas *c1, vector<double> *x,vector<double> *y,vector<double> *z,vector<double> *q, vector<double> *x_out,vector<double> *y_out,vector<double> *z_out,vector<double> *q_out) {
void NPL::Tracking::Hough_3D(vector<double> *x,vector<double> *y,vector<double> *z,vector<double> *q,vector<double> *x_out,vector<double> *y_out,vector<double> *z_out,vector<double> *q_out) {
int nt_xy=180;
- int nt_xz=180;
- int nt_yz=180;
- int nr_xy=45;
- int nr_xz=300;
- int nr_yz=300;
- double bint_xy=2.;
- double bint_xz=2.;
- double bint_yz=2.;
- double binr_xy=3.;
- double binr_xz=3.;
- double binr_yz=3.;
- int nt,nr;
- double PI = TMath::Pi();
-
- double rho_xy,rho_xz,rho_yz;
- double theta_xy,theta_xz,theta_yz;
-
- TH2F *hp_xy = new TH2F("hp_xy","hp_xy",nt_xy,0,180,nr_xy,-1*nr_xy,nr_xy);
- TH2F *hp_xz = new TH2F("hp_xz","hp_xz",nt_xz,0,180,nr_xz,-1*nr_xz,nr_xz);
- TH2F *hp_yz = new TH2F("hp_yz","hp_yz",nt_yz,0,180,nr_yz,-1*nr_yz,nr_yz);
-
+ int nt_xz=180;
+ int nt_yz=180;
+ int nr_xy=45;
+ int nr_xz=300;
+ int nr_yz=300;
+ double bint_xy=2.;
+ double bint_xz=2.;
+ double bint_yz=2.;
+ double binr_xy=3.;
+ double binr_xz=3.;
+ double binr_yz=3.;
+ int nt,nr;
+ double PI = TMath::Pi();
+
+ double rho_xy,rho_xz,rho_yz;
+ double theta_xy,theta_xz,theta_yz;
+
+ TH2F *hp_xy = new TH2F("hp_xy","hp_xy",nt_xy,0,180,nr_xy,-1*nr_xy,nr_xy);
+ TH2F *hp_xz = new TH2F("hp_xz","hp_xz",nt_xz,0,180,nr_xz,-1*nr_xz,nr_xz);
+ TH2F *hp_yz = new TH2F("hp_yz","hp_yz",nt_yz,0,180,nr_yz,-1*nr_yz,nr_yz);
+
// int npeaks_xy, npeaks_xz, npeaks_yz;
vector<double> thetapeaks_xy, rpeaks_xy, thetapeaks_xz, rpeaks_xz, thetapeaks_yz, rpeaks_yz;
- double max_xy, max_xz, max_yz;
- double rmean_xy=0, thetamean_xy=0, rmean_xz=0, thetamean_xz=0, rmean_yz=0, thetamean_yz=0;
-
- double r0_xy=0., r0_xz=0., r0_yz=0., rmin_xy=0., rmin_xz=0., rmin_yz=0., rmax_xy=0., rmax_xz=0., rmax_yz=0.;
- double tmin=0., tmax=0.;
- double rinf=0., rsup=0.;
-
- nt=nt_xy;
- nr=nr_xy;
- if(nt<nt_xz)nt=nt_xz;
- if(nr<nr_xz)nr=nr_xz;
- if(nt<nt_yz)nt=nt_yz;
- if(nr<nr_yz)nr=nr_yz;
-
+ double max_xy, max_xz, max_yz;
+ double rmean_xy=0, thetamean_xy=0, rmean_xz=0, thetamean_xz=0, rmean_yz=0, thetamean_yz=0;
+
+ double r0_xy=0., r0_xz=0., r0_yz=0., rmin_xy=0., rmin_xz=0., rmin_yz=0., rmax_xy=0., rmax_xz=0., rmax_yz=0.;
+ double tmin=0., tmax=0.;
+ double rinf=0., rsup=0.;
+
+ nt=nt_xy;
+ nr=nr_xy;
+ if(nt<nt_xz)nt=nt_xz;
+ if(nr<nr_xz)nr=nr_xz;
+ if(nt<nt_yz)nt=nt_yz;
+ if(nr<nr_yz)nr=nr_yz;
+
for(unsigned int i=0;i<x->size();i++)
- {
- //Fill coordinate space histograms for plots
- //Loop of indices and fill Histograms
+ {
+ //Fill coordinate space histograms for plots
+ //Loop of indices and fill Histograms
for(int j=0;j<nt;j++)
- {
-
+ {
//xy
- theta_xy = j*180./nt_xy;
-
- /* cout <<y->at(i) << endl; */
+ theta_xy = j*180./nt_xy;
rho_xy = x->at(i)*TMath::Cos(theta_xy*PI/180.)+y->at(i)*TMath::Sin(theta_xy*PI/180.);
if(abs(theta_xy)<180. && abs(rho_xy)<nr_xy)
- {
- hp_xy->Fill(theta_xy,rho_xy);
- }
-
- //xz
- theta_xz = j*180./nt_xz;
- rho_xz = z->at(i)*TMath::Cos(theta_xz*PI/180.)+x->at(i)*TMath::Sin(theta_xz*PI/180.);
- if(abs(theta_xz)<180. && abs(rho_xz)<nr_xz)
- {
- hp_xz->Fill(theta_xz,rho_xz);
- }
-
- //yz
- theta_yz = j*180./nt_yz;
- rho_yz = z->at(i)*TMath::Cos(theta_yz*PI/180.)+y->at(i)*TMath::Sin(theta_yz*PI/180.);
- if(abs(theta_yz)<180. && abs(rho_yz)<nr_yz)
- {
- hp_yz->Fill(theta_yz,rho_yz);
- }
- }
- }
-
- max_xy = hp_xy->GetMaximum();
- max_xz = hp_xz->GetMaximum();
- max_yz = hp_yz->GetMaximum();
-
- for(int ii=0; ii<nt; ii++)
- {
- for(int jj=0; jj<nr; jj++)
- {
- if(hp_xy->GetBinContent(ii+1, jj+1) == max_xy && jj<nr_xy)
- {
- thetapeaks_xy.push_back((ii+0.5)*nt_xy/nt);
- rpeaks_xy.push_back((jj+0.5)*2 - nr_xy);
- rmean_xy += rpeaks_xy.back();
- thetamean_xy += thetapeaks_xy.back();
- }
- if(hp_xz->GetBinContent(ii+1, jj+1) == max_xz)
- {
- thetapeaks_xz.push_back((ii+0.5)*nt_xz/nt);
- rpeaks_xz.push_back((jj+0.5)*2 - nr_xz);
- rmean_xz += rpeaks_xz.back();
- thetamean_xz += thetapeaks_xz.back();
- }
- if(hp_yz->GetBinContent(ii+1, jj+1) == max_yz)
- {
- thetapeaks_yz.push_back((ii+0.5)*nt_yz/nt);
- rpeaks_yz.push_back((jj+0.5)*2 - nr_yz);
- rmean_yz += rpeaks_yz.back();
- thetamean_yz += thetapeaks_yz.back();
- }
- }
- }
-
- // xy PEAK
- rmean_xy = rmean_xy / rpeaks_xy.size();
- thetamean_xy = thetamean_xy / thetapeaks_xy.size();
-
- // xz PEAK
- rmean_xz = rmean_xz / rpeaks_xz.size();
- thetamean_xz = thetamean_xz / thetapeaks_xz.size();
-
- // yz PEAK
- rmean_yz = rmean_yz / rpeaks_yz.size();
- thetamean_yz = thetamean_yz / thetapeaks_yz.size();
-
- rmean_xy = rpeaks_xy[0];
- thetamean_xy = thetapeaks_xy[0];
- rmean_xz = rpeaks_xz[0];
- thetamean_xz = thetapeaks_xz[0];
- rmean_yz = rpeaks_yz[0];
- thetamean_yz = thetapeaks_yz[0];
-
- //Selection of x,y,z points COMMON to the 3 maxmean+/-1 found in Hough spaces for xy, xz and yz spaces
- for(unsigned int i=0;i<x->size();i++)
- {
- r0_xy = x->at(i)*TMath::Cos(thetamean_xy*PI/180.)+y->at(i)*TMath::Sin(thetamean_xy*PI/180.);
- tmin = thetamean_xy-bint_xy;
- tmax = thetamean_xy+bint_xy;
- if((tmin)<0) tmin = tmin + 180.;
- if((tmax)>180) tmax = tmax - 180.;
- rmin_xy = x->at(i)*TMath::Cos(tmin*PI/180.)+y->at(i)*TMath::Sin(tmin*PI/180.);
- rmax_xy = x->at(i)*TMath::Cos(tmax*PI/180.)+y->at(i)*TMath::Sin(tmax*PI/180.);
-
- rinf = min( rmean_xy - binr_xy, rmean_xy + binr_xy);
- rsup = max( rmean_xy - binr_xy, rmean_xy + binr_xy);
- if((r0_xy>=rinf || rmin_xy>=rinf || rmax_xy>=rinf) && (r0_xy<=rsup || rmin_xy<=rsup || rmax_xy<=rsup))
- {
- r0_xz = z->at(i)*TMath::Cos(thetamean_xz*PI/180.)+x->at(i)*TMath::Sin(thetamean_xz*PI/180.);
- tmin = thetamean_xz-bint_xz;
- tmax = thetamean_xz+bint_xz;
- if((tmin)<0) tmin = tmin + 180.;
- if((tmax)>180) tmax = tmax - 180.;
- rmin_xz = z->at(i)*TMath::Cos(tmin*PI/180.)+x->at(i)*TMath::Sin(tmin*PI/180.);
- rmax_xz = z->at(i)*TMath::Cos(tmax*PI/180.)+x->at(i)*TMath::Sin(tmax*PI/180.);
-
- rinf = min( rmean_xz - binr_xz, rmean_xz + binr_xz);
- rsup = max( rmean_xz - binr_xz, rmean_xz + binr_xz);
-
- if((r0_xz>=rinf || rmin_xz>=rinf || rmax_xz>=rinf) && (r0_xz<=rsup || rmin_xz<=rsup || rmax_xz<=rsup))
- {
- r0_yz = z->at(i)*TMath::Cos(thetamean_yz*PI/180.)+y->at(i)*TMath::Sin(thetamean_yz*PI/180.);
- tmin = thetamean_yz-bint_yz;
- tmax = thetamean_yz+bint_yz;
- if((tmin)<0) tmin = tmin + 180.;
- if((tmax)>180) tmax = tmax - 180.;
- rmin_yz = z->at(i)*TMath::Cos(tmin*PI/180.)+y->at(i)*TMath::Sin(tmin*PI/180.);
- rmax_yz = z->at(i)*TMath::Cos(tmax*PI/180.)+y->at(i)*TMath::Sin(tmax*PI/180.);
-
- rinf = min( rmean_yz - binr_yz, rmean_yz + binr_yz);
- rsup = max( rmean_yz - binr_yz, rmean_yz + binr_yz);
-
- if((r0_yz>=rinf || rmin_yz>=rinf || rmax_yz>=rinf) && (r0_yz<=rsup || rmin_yz<=rsup || rmax_yz<=rsup))
- {
- x_out->push_back(x->at(i));
- y_out->push_back(y->at(i));
- z_out->push_back(z->at(i));
- q_out->push_back(q->at(i));
- }
- }
- }
-
- }
- delete hp_xy;
- delete hp_xz;
- delete hp_yz;
-
+ {
+ hp_xy->Fill(theta_xy,rho_xy);
+ }
+
+ //xz
+ theta_xz = j*180./nt_xz;
+ rho_xz = z->at(i)*TMath::Cos(theta_xz*PI/180.)+x->at(i)*TMath::Sin(theta_xz*PI/180.);
+ if(abs(theta_xz)<180. && abs(rho_xz)<nr_xz)
+ {
+ hp_xz->Fill(theta_xz,rho_xz);
+ }
+
+ //yz
+ theta_yz = j*180./nt_yz;
+ rho_yz = z->at(i)*TMath::Cos(theta_yz*PI/180.)+y->at(i)*TMath::Sin(theta_yz*PI/180.);
+ if(abs(theta_yz)<180. && abs(rho_yz)<nr_yz)
+ {
+ hp_yz->Fill(theta_yz,rho_yz);
+ }
+ }
+ }
+
+ max_xy = hp_xy->GetMaximum();
+ max_xz = hp_xz->GetMaximum();
+ max_yz = hp_yz->GetMaximum();
+
+ for(int ii=0; ii<nt; ii++)
+ {
+ for(int jj=0; jj<nr; jj++)
+ {
+ if(hp_xy->GetBinContent(ii+1, jj+1) == max_xy && jj<nr_xy)
+ {
+ thetapeaks_xy.push_back((ii+0.5)*nt_xy/nt);
+ rpeaks_xy.push_back((jj+0.5)*2 - nr_xy);
+ rmean_xy += rpeaks_xy.back();
+ thetamean_xy += thetapeaks_xy.back();
+ }
+ if(hp_xz->GetBinContent(ii+1, jj+1) == max_xz)
+ {
+ thetapeaks_xz.push_back((ii+0.5)*nt_xz/nt);
+ rpeaks_xz.push_back((jj+0.5)*2 - nr_xz);
+ rmean_xz += rpeaks_xz.back();
+ thetamean_xz += thetapeaks_xz.back();
+ }
+ if(hp_yz->GetBinContent(ii+1, jj+1) == max_yz)
+ {
+ thetapeaks_yz.push_back((ii+0.5)*nt_yz/nt);
+ rpeaks_yz.push_back((jj+0.5)*2 - nr_yz);
+ rmean_yz += rpeaks_yz.back();
+ thetamean_yz += thetapeaks_yz.back();
+ }
+ }
+ }
+
+ // xy PEAK
+ rmean_xy = rmean_xy / rpeaks_xy.size();
+ thetamean_xy = thetamean_xy / thetapeaks_xy.size();
+
+ // xz PEAK
+ rmean_xz = rmean_xz / rpeaks_xz.size();
+ thetamean_xz = thetamean_xz / thetapeaks_xz.size();
+
+ // yz PEAK
+ rmean_yz = rmean_yz / rpeaks_yz.size();
+ thetamean_yz = thetamean_yz / thetapeaks_yz.size();
+
+ rmean_xy = rpeaks_xy[0];
+ thetamean_xy = thetapeaks_xy[0];
+ rmean_xz = rpeaks_xz[0];
+ thetamean_xz = thetapeaks_xz[0];
+ rmean_yz = rpeaks_yz[0];
+ thetamean_yz = thetapeaks_yz[0];
+
+ //Selection of x,y,z points COMMON to the 3 maxmean+/-1 found in Hough spaces for xy, xz and yz spaces
+ for(unsigned int i=0;i<x->size();i++)
+ {
+ r0_xy = x->at(i)*TMath::Cos(thetamean_xy*PI/180.)+y->at(i)*TMath::Sin(thetamean_xy*PI/180.);
+ tmin = thetamean_xy-bint_xy;
+ tmax = thetamean_xy+bint_xy;
+ if((tmin)<0) tmin = tmin + 180.;
+ if((tmax)>180) tmax = tmax - 180.;
+ rmin_xy = x->at(i)*TMath::Cos(tmin*PI/180.)+y->at(i)*TMath::Sin(tmin*PI/180.);
+ rmax_xy = x->at(i)*TMath::Cos(tmax*PI/180.)+y->at(i)*TMath::Sin(tmax*PI/180.);
+
+ rinf = min( rmean_xy - binr_xy, rmean_xy + binr_xy);
+ rsup = max( rmean_xy - binr_xy, rmean_xy + binr_xy);
+ if((r0_xy>=rinf || rmin_xy>=rinf || rmax_xy>=rinf) && (r0_xy<=rsup || rmin_xy<=rsup || rmax_xy<=rsup))
+ {
+ r0_xz = z->at(i)*TMath::Cos(thetamean_xz*PI/180.)+x->at(i)*TMath::Sin(thetamean_xz*PI/180.);
+ tmin = thetamean_xz-bint_xz;
+ tmax = thetamean_xz+bint_xz;
+ if((tmin)<0) tmin = tmin + 180.;
+ if((tmax)>180) tmax = tmax - 180.;
+ rmin_xz = z->at(i)*TMath::Cos(tmin*PI/180.)+x->at(i)*TMath::Sin(tmin*PI/180.);
+ rmax_xz = z->at(i)*TMath::Cos(tmax*PI/180.)+x->at(i)*TMath::Sin(tmax*PI/180.);
+
+ rinf = min( rmean_xz - binr_xz, rmean_xz + binr_xz);
+ rsup = max( rmean_xz - binr_xz, rmean_xz + binr_xz);
+
+ if((r0_xz>=rinf || rmin_xz>=rinf || rmax_xz>=rinf) && (r0_xz<=rsup || rmin_xz<=rsup || rmax_xz<=rsup))
+ {
+ r0_yz = z->at(i)*TMath::Cos(thetamean_yz*PI/180.)+y->at(i)*TMath::Sin(thetamean_yz*PI/180.);
+ tmin = thetamean_yz-bint_yz;
+ tmax = thetamean_yz+bint_yz;
+ if((tmin)<0) tmin = tmin + 180.;
+ if((tmax)>180) tmax = tmax - 180.;
+ rmin_yz = z->at(i)*TMath::Cos(tmin*PI/180.)+y->at(i)*TMath::Sin(tmin*PI/180.);
+ rmax_yz = z->at(i)*TMath::Cos(tmax*PI/180.)+y->at(i)*TMath::Sin(tmax*PI/180.);
+
+ rinf = min( rmean_yz - binr_yz, rmean_yz + binr_yz);
+ rsup = max( rmean_yz - binr_yz, rmean_yz + binr_yz);
+
+ if((r0_yz>=rinf || rmin_yz>=rinf || rmax_yz>=rinf) && (r0_yz<=rsup || rmin_yz<=rsup || rmax_yz<=rsup))
+ {
+ x_out->push_back(x->at(i));
+ y_out->push_back(y->at(i));
+ z_out->push_back(z->at(i));
+ q_out->push_back(q->at(i));
+ }
+ }
+ }
+
+ }
+ delete hp_xy;
+ delete hp_xz;
+ delete hp_yz;
+
}
// Calculation of the minimal distance between 2 lines in 3D space & calculation of mid-point=>vertex of interaction
-void NPL::Tracking::vertex(double *p, double *pp, double &xv,double &yv,double &zv,double &min_dist, double &Theta_tr1, double &Theta_tr2) {
- double a1 = p[0];
- double a2 = p[2];
- double b1 = p[1];
- double b2 = p[3];
- double ap1 = pp[0];
- double ap2 = pp[2];
- double bp1 = pp[1];
- double bp2 = pp[3];
-
- double alpha, beta, A, B, C;
-
- alpha = (bp1*(a1-ap1)+bp2*(a2-ap2))/(bp1*bp1 + bp2*bp2 + 1);
- beta = (bp1*b1+bp2*b2+1)/(bp1*bp1 + bp2*bp2 + 1);
-
- A = beta*(bp1*bp1 + bp2*bp2 + 1) - (bp1*b1 + bp2*b2 + 1);
- B = (b1*b1 + b2*b2 + 1) - beta*(bp1*b1+bp2*b2+1);
- C = beta*(bp1*(ap1-a1) + bp2*(ap2-a2)) - (b1*(ap1-a1) + b2*(ap2-a2));
-
- double sol1, solf1;
- double x,y,z,xp,yp,zp;
-
-
- sol1 = -(A*alpha + C)/(A*beta + B);
- solf1 = alpha + beta* sol1;
-
- x = a1 + b1*sol1;
- y = a2 + b2*sol1;
- z = sol1;
- xp = ap1 + bp1*solf1;
- yp = ap2 + bp2*solf1;
- zp = solf1;
-
- xv = (x+xp)/2.;
- yv = (y+yp)/2.;
- zv = (z+zp)/2.;
-
- Theta_tr1 = TMath::ATan(TMath::Sqrt((x-a1)*(x-a1)+(y-a2)*(y-a2))/TMath::Abs(sol1));
- Theta_tr2 = TMath::ATan(TMath::Sqrt((xp-ap1)*(xp-ap1)+(yp-ap2)*(yp-ap2))/TMath::Abs(solf1));
- Theta_tr1 = Theta_tr1*180./TMath::Pi();
- Theta_tr2 = Theta_tr2*180./TMath::Pi();
-
- //cout << "Vertex 1st :" << x << "," << y << "," << z << endl;
- //cout << "Vertex 2nd :" << xp << "," << yp << "," << zp << endl;
- //cout << "Vertex middle :" << xv << "," << yv << "," << zv << endl;
-
- //cout << "min dist " << sqrt(pow((x-xp),2) + pow((y-yp),2) + pow((z-zp),2)) << endl;
-
+void NPL::Tracking::vertex(double *p, double *pp, double &xv,double &yv,double &zv,double &min_dist, double &Theta_tr1, double &Theta_tr2, double &Phi1, double &Phi2) {
+ double a1 = p[0];
+ double a2 = p[2];
+ double b1 = p[1];
+ double b2 = p[3];
+ double ap1 = pp[0];
+ double ap2 = pp[2];
+ double bp1 = pp[1];
+ double bp2 = pp[3];
+
+ // calcul the closest pointis between track1 && track2
+ double alpha, beta, A, B, C;
+
+ alpha = (bp1*(a1-ap1)+bp2*(a2-ap2))/(bp1*bp1 + bp2*bp2 + 1);
+ beta = (bp1*b1+bp2*b2+1)/(bp1*bp1 + bp2*bp2 + 1);
+
+ A = beta*(bp1*bp1 + bp2*bp2 + 1) - (bp1*b1 + bp2*b2 + 1);
+ B = (b1*b1 + b2*b2 + 1) - beta*(bp1*b1+bp2*b2+1);
+ C = beta*(bp1*(ap1-a1) + bp2*(ap2-a2)) - (b1*(ap1-a1) + b2*(ap2-a2));
+
+ double sol1, solf1;
+ double x,y,z,xp,yp,zp;
+
+ sol1 = -(A*alpha + C)/(A*beta + B);
+ solf1 = alpha + beta* sol1;
+
+ // point from track1
+ x = a1 + b1*sol1;
+ y = a2 + b2*sol1;
+ z = sol1;
+ // point from track2
+ xp = ap1 + bp1*solf1;
+ yp = ap2 + bp2*solf1;
+ zp = solf1;
+ // vertex (mid-ditance between the 2 points)
+ xv = (x+xp)/2.;
+ yv = (y+yp)/2.;
+ zv = (z+zp)/2.;
+
+ // calulate Theta and Phi
+ double xa,ya,za,zap,xap,yap;
+ double zpoint = 3000;
+ xa = a1 + b1*zpoint;
+ ya = a2 + b2*zpoint;
+
+ xap = ap1 + bp1*zpoint;
+ yap = ap2 + bp2*zpoint;
+
+ /* Theta_tr1 = TMath::ACos((zpoint-z)/TMath::Sqrt((xa-x)*(xa-x)+(ya-y)*(ya-y)+(zpoint-z)*(zpoint-z))); */
+ /* Theta_tr2 = TMath::ACos((zpoint-zp)/TMath::Sqrt((xap-xp)*(xap-xp)+(yap-yp)*(yap-yp)+(zpoint-zp)*(zpoint-zp))); */
+ // Aldric Revel version :
+ Theta_tr1 = TMath::ATan(TMath::Sqrt((x-a1)*(x-a1)+(y-a2)*(y-a2))/TMath::Abs(sol1));
+ Theta_tr2 = TMath::ATan(TMath::Sqrt((xp-ap1)*(xp-ap1)+(yp-ap2)*(yp-ap2))/TMath::Abs(solf1));
+ Phi1 = TMath::ATan2(b2,b1);
+ Phi2 = TMath::ATan2(bp2,bp1);
+
+ /* Phi1 = TMath::RadToDeg()*TMath::ATan2(parFit4->at(0),parFit2->at(0)); */
+ /* if( xa-x > 0 && ya-y >= 0) */
+ /* Phi1 = TMath::ATan((ya-y)/(xa-x)); */
+ /* else if(xa-x > 0 && ya-y < 0) */
+ /* Phi1 = TMath::ATan((ya-y)/(xa-x)) + 2*TMath::Pi(); */
+ /* else if(xa-x < 0) */
+ /* Phi1 = TMath::ATan((ya-y)/(xa-x)) + TMath::Pi(); */
+
+ /* if(xap-xp > 0 && yap-yp >= 0) */
+ /* Phi2 = TMath::ATan((yap-yp)/(xap-xp)); */
+ /* else if(xap-xp > 0 && yap-yp < 0) */
+ /* Phi2 = TMath::ATan((yap-yp)/(xap-xp)) + 2*TMath::Pi(); */
+ /* else if(xap-xp < 0) */
+ /* Phi2 = TMath::ATan((yap-yp)/(xap-xp)) + TMath::Pi(); */
+
+ /* if(Phi1>TMath::Pi()) Phi1 = Phi1-2*TMath::Pi(); */
+ /* if(Phi2>TMath::Pi()) Phi2 = Phi2-2*TMath::Pi(); */
+
+ Phi1 = 180*Phi1/TMath::Pi();
+ Phi2 = 180*Phi2/TMath::Pi();
+ Theta_tr1 = Theta_tr1*180./TMath::Pi();
+ Theta_tr2 = Theta_tr2*180./TMath::Pi();
+
+ min_dist = sqrt(pow((x-xp),2) + pow((y-yp),2) + pow((z-zp),2));
+
}
void NPL::Tracking::ParFor_Vertex(double *a, double *b, double *parFit) {
@@ -478,7 +497,7 @@ void NPL::Tracking::ParFor_Vertex(double *a, double *b, double *parFit) {
APX=a[0]-b[0]; APY=a[1]-b[1]; APZ=a[2]-b[2];
APX0=APX/APZ; APY0=APY/APZ;
X0=a[0]-a[2]*APX0; Y0=a[1]-a[2]*APY0; //Z0=0
-
+
parFit[0]=X0;
parFit[1]=APX0;
parFit[2]=Y0;
diff --git a/NPLib/TrackReconstruction/Tracking.h b/NPLib/TrackReconstruction/Tracking.h
index ee158a5bb..a198da8a9 100644
--- a/NPLib/TrackReconstruction/Tracking.h
+++ b/NPLib/TrackReconstruction/Tracking.h
@@ -23,7 +23,7 @@ class Tracking {
void FindStart(double pStart[4], double chi[2], int fitStatus[2],TGraph *grxz, TGraph *gryz);
double distance2(double x,double y,double z, double *p);
void Hough_3D(vector<double> *x,vector<double> *y,vector<double> *z,vector<double> *q,vector<double> *x_out,vector<double> *y_out,vector<double> *z_out,vector<double> *q_out);
- void vertex(double *p, double *pp, double &xv,double &yv,double &zv, double &min_dist, double &Theta_tr1, double &Theta_tr2);
+ void vertex(double *p, double *pp, double &xv,double &yv,double &zv, double &min_dist, double &Theta_tr1, double &Theta_tr2, double &Phi1, double &Phi2);
void ParFor_Vertex(double *a, double *b, double *parFit);
ClassDef(Tracking,1);
diff --git a/NPSimulation/Core/MaterialManager.cc b/NPSimulation/Core/MaterialManager.cc
index 7a94a28d6..0158e7b50 100644
--- a/NPSimulation/Core/MaterialManager.cc
+++ b/NPSimulation/Core/MaterialManager.cc
@@ -63,11 +63,11 @@ void MaterialManager::Destroy() {}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void MaterialManager::ClearMaterialLibrary() {
- // map<string, G4Material*>::iterator it;
- // for (it = m_Material.begin(); it != m_Material.end(); it++) {
- // delete it->second;
- // }
-
+ // map<string, G4Material*>::iterator it;
+ // for (it = m_Material.begin(); it != m_Material.end(); it++) {
+ // delete it->second;
+ // }
+
// Geant4 own pointer to the material
// we can forget about them but not delete it
// as they can be deleted by the kernel e.g. when Cleaning the PVPStore
@@ -80,7 +80,7 @@ void MaterialManager::ClearMaterialLibrary() {
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4Material* MaterialManager::GetMaterialFromLibrary(string Name,
- double density) {
+ double density) {
// Search if the material is already instantiate
map<string, G4Material*>::iterator it;
it = m_Material.find(Name);
@@ -147,11 +147,12 @@ G4Material* MaterialManager::GetMaterialFromLibrary(string Name,
else if (Name == "Kapton") {
if (!density)
- density = 1.39 * g / cm3;
- G4Material* material = new G4Material("NPS_" + Name, density, 3);
- material->AddElement(GetElementFromLibrary("H"), 4);
- material->AddElement(GetElementFromLibrary("C"), 10);
- material->AddElement(GetElementFromLibrary("O"), 2);
+ density = 1.42 * g / cm3;
+ G4Material* material = new G4Material("NPS_" + Name, density, 4);
+ material->AddElement(GetElementFromLibrary("H"), 0.026);
+ material->AddElement(GetElementFromLibrary("C"), 0.69);
+ material->AddElement(GetElementFromLibrary("O"), 0.21);
+ material->AddElement(GetElementFromLibrary("N"), 0.074);
m_Material[Name] = material;
return material;
}
@@ -197,7 +198,7 @@ G4Material* MaterialManager::GetMaterialFromLibrary(string Name,
if (!density)
density = 0.808 * g / cm3;
G4Material* material = new G4Material("NPS_" + Name, 7, 14.01 * g / mole,
- density, kStateLiquid, 77 * kelvin);
+ density, kStateLiquid, 77 * kelvin);
m_Material[Name] = material;
return material;
}
@@ -315,6 +316,57 @@ G4Material* MaterialManager::GetMaterialFromLibrary(string Name,
m_Material[Name] = material;
return material;
}
+
+ else if (Name=="MgO"){ //cyril
+ if (!density)
+ density = 3.6* g / cm3;
+ G4Material* material = new G4Material("NPS_" + Name, density, 2);
+ material->AddElement(GetElementFromLibrary("Mg"), 1);
+ material->AddElement(GetElementFromLibrary("O"), 1);
+ m_Material[Name] = material;
+ return material;
+ }
+ else if (Name=="mixMINOS"){ //cyril
+ if (!density)
+ density = 0.0019836* g / cm3;
+ G4Material* material = new G4Material("NPS_" + Name, density, 3);
+ material->AddMaterial(GetMaterialFromLibrary("CF4"), .15);
+ material->AddMaterial(GetMaterialFromLibrary("isobutane"), .03);
+ material->AddElement(GetElementFromLibrary("Ar"), .82);
+ m_Material[Name] = material;
+ return material;
+ }
+ else if (Name=="mumetal"){ //cyril
+ if (!density)
+ density = 8.7* g / cm3;
+ G4Material* material = new G4Material("NPS_" + Name, density, 3);
+ material->AddElement(GetElementFromLibrary("Ni"), .8);
+ material->AddElement(GetElementFromLibrary("Fe"), .15);
+ material->AddElement(GetElementFromLibrary("Mo"), .05);
+ m_Material[Name] = material;
+ return material;
+ }
+ else if (Name=="LH2"){ //cyril
+ if (!density)
+ density = 0.07293* g / cm3;
+ G4Material* material = new G4Material("NPS_" + Name, density, 1);
+ material->AddElement(GetElementFromLibrary("H"), 2);
+ m_Material[Name] = material;
+ return material;
+ }
+ else if (Name=="Rohacell"){ //cyril
+ if (!density)
+ density = 0.075* g / cm3;
+ G4Material* material = new G4Material("NPS_" + Name, density, 4);
+ material->AddElement(GetElementFromLibrary("H"), 0.0805);
+ material->AddElement(GetElementFromLibrary("C"), 0.6014);
+ material->AddElement(GetElementFromLibrary("O"), 0.3154);
+ material->AddElement(GetElementFromLibrary("N"), 0.00276);
+ m_Material[Name] = material;
+ return material;
+ }
+
+
// Usual detector material
else if (Name == "Si") {
if (!density)
@@ -456,13 +508,13 @@ G4Material* MaterialManager::GetMaterialFromLibrary(string Name,
wl = wlmin + i * step;
// Formula from www.refractiveindex.info
rindex[i] = sqrt(1 + 0.27587 + 0.68689 / (1 - pow(0.130 / wl, 2))
- + 0.26090 / (1 - pow(0.147 / wl, 2))
- + 0.06256 / (1 - pow(0.163 / wl, 2))
- + 0.06527 / (1 - pow(0.177 / wl, 2))
- + 0.14991 / (1 - pow(0.185 / wl, 2))
- + 0.51818 / (1 - pow(0.206 / wl, 2))
- + 0.01918 / (1 - pow(0.218 / wl, 2))
- + 3.38229 / (1 - pow(161.29 / wl, 2)));
+ + 0.26090 / (1 - pow(0.147 / wl, 2))
+ + 0.06256 / (1 - pow(0.163 / wl, 2))
+ + 0.06527 / (1 - pow(0.177 / wl, 2))
+ + 0.14991 / (1 - pow(0.185 / wl, 2))
+ + 0.51818 / (1 - pow(0.206 / wl, 2))
+ + 0.01918 / (1 - pow(0.218 / wl, 2))
+ + 3.38229 / (1 - pow(161.29 / wl, 2)));
// check below block
energy_r[i] = h_Planck * c_light / wl;
// To be defined properly
@@ -476,7 +528,7 @@ G4Material* MaterialManager::GetMaterialFromLibrary(string Name,
MPT->AddProperty("SCINTILLATION", energy_e, scint, 5); // check
MPT->AddProperty("RINDEX", energy_r, rindex, NumberOfPoints); // check
MPT->AddProperty("ABSLENGTH", energy_r, absorption,
- NumberOfPoints); // check
+ NumberOfPoints); // check
MPT->AddProperty("FASTCOMPONENT", energy_e, fast, 2.1); // good
MPT->AddProperty("SLOWCOMPONENT", energy_e, slow, 22.6); // good
MPT->AddConstProperty("RESOLUTIONSCALE", 1.0); // check
@@ -570,13 +622,13 @@ G4Material* MaterialManager::GetMaterialFromLibrary(string Name,
wl = wlmin + i * step;
// Formula from www.refractiveindex.info
rindex[i] = sqrt(1 + 0.27587 + 0.68689 / (1 - pow(0.130 / wl, 2))
- + 0.26090 / (1 - pow(0.147 / wl, 2))
- + 0.06256 / (1 - pow(0.163 / wl, 2))
- + 0.06527 / (1 - pow(0.177 / wl, 2))
- + 0.14991 / (1 - pow(0.185 / wl, 2))
- + 0.51818 / (1 - pow(0.206 / wl, 2))
- + 0.01918 / (1 - pow(0.218 / wl, 2))
- + 3.38229 / (1 - pow(161.29 / wl, 2)));
+ + 0.26090 / (1 - pow(0.147 / wl, 2))
+ + 0.06256 / (1 - pow(0.163 / wl, 2))
+ + 0.06527 / (1 - pow(0.177 / wl, 2))
+ + 0.14991 / (1 - pow(0.185 / wl, 2))
+ + 0.51818 / (1 - pow(0.206 / wl, 2))
+ + 0.01918 / (1 - pow(0.218 / wl, 2))
+ + 3.38229 / (1 - pow(161.29 / wl, 2)));
energy_r[i] = h_Planck * c_light / wl;
// To be defined properly
@@ -674,7 +726,7 @@ G4Material* MaterialManager::GetMaterialFromLibrary(string Name,
if (!density)
density = 1.74 * mg / cm3;
G4Material* material
- = new G4Material("NPS_" + Name, density, 3); //@ 0K, 1 atm
+ = new G4Material("NPS_" + Name, density, 3); //@ 0K, 1 atm
material->AddElement(GetElementFromLibrary("Ar"), 0.9222);
material->AddElement(GetElementFromLibrary("C"), 0.0623);
material->AddElement(GetElementFromLibrary("H"), 0.0155);
@@ -686,7 +738,7 @@ G4Material* MaterialManager::GetMaterialFromLibrary(string Name,
if (!density)
density = 0.57 * mg / cm3;
G4Material* material
- = new G4Material("NPS_" + Name, density, 3); //@ 0K, 1/3 atm
+ = new G4Material("NPS_" + Name, density, 3); //@ 0K, 1/3 atm
material->AddElement(GetElementFromLibrary("Ar"), 0.9222);
material->AddElement(GetElementFromLibrary("C"), 0.0623);
material->AddElement(GetElementFromLibrary("H"), 0.0155);
@@ -704,13 +756,23 @@ G4Material* MaterialManager::GetMaterialFromLibrary(string Name,
return material;
}
+ else if(Name == "iC4H10" || Name == "Isobutane" || Name == "isobutane"){
+ density = 0.002506*g/cm3;
+ G4Material* material = new G4Material("NPS_"+Name,density,2);
+ material->AddElement(GetElementFromLibrary("C"), 4);
+ material->AddElement(GetElementFromLibrary("H"), 10);
+ m_Material[Name]=material;
+ return material;
+ }
+
else if (Name == "CF4") { // 52 torr
if (!density)
density = 3.78 * mg / cm3;
G4Material* material = new G4Material("NPS_" + Name, density, 2,
- kStateGas, 300, 0.0693276 * bar);
+ kStateGas, 300, 0.0693276 * bar);
material->AddElement(GetElementFromLibrary("C"), 1);
material->AddElement(GetElementFromLibrary("F"), 4);
+ material->GetIonisation()->SetMeanExcitationEnergy(20.0*eV);
m_Material[Name] = material;
return material;
}
@@ -740,19 +802,19 @@ G4Material* MaterialManager::GetMaterialFromLibrary(string Name,
const G4int NUMENTRIES = 15;
G4double PMMA_PP[NUMENTRIES]
- = {10 * eV, 3.25 * eV, 3.099 * eV, 2.88 * eV, 2.695 * eV,
- 2.53 * eV, 2.38 * eV, 2.254 * eV, 2.138 * eV, 2.033 * eV,
- 1.937 * eV, 1.859 * eV, 1.771 * eV, 1.6 * eV, 0 * eV};
+ = {10 * eV, 3.25 * eV, 3.099 * eV, 2.88 * eV, 2.695 * eV,
+ 2.53 * eV, 2.38 * eV, 2.254 * eV, 2.138 * eV, 2.033 * eV,
+ 1.937 * eV, 1.859 * eV, 1.771 * eV, 1.6 * eV, 0 * eV};
G4double PMMA_RIND[NUMENTRIES]
- = {1.47, 1.47, 1.47, 1.47, 1.47, 1.47, 1.47, 1.47,
- 1.47, 1.47, 1.47, 1.47, 1.47, 1.47, 1.47};
+ = {1.47, 1.47, 1.47, 1.47, 1.47, 1.47, 1.47, 1.47,
+ 1.47, 1.47, 1.47, 1.47, 1.47, 1.47, 1.47};
G4double PMMA_ABSL[NUMENTRIES]
- = {35. * cm, 35. * cm, 35. * cm, 35. * cm, 35. * cm,
- 35. * cm, 35. * cm, 35. * cm, 35. * cm, 35. * cm,
- 35. * cm, 35. * cm, 35. * cm, 35. * cm, 35. * cm};
+ = {35. * cm, 35. * cm, 35. * cm, 35. * cm, 35. * cm,
+ 35. * cm, 35. * cm, 35. * cm, 35. * cm, 35. * cm,
+ 35. * cm, 35. * cm, 35. * cm, 35. * cm, 35. * cm};
/*G4double THICK_ABSL[NUMENTRIES] =
- {0.01*mm,0.01*mm,0.01*mm,0.01*mm,0.01*mm,0.01*mm,0.01*mm,0.01*mm,
- 0.01*mm,0.01*mm,0.01*mm,0.01*mm,0.01*mm,0.01*mm,0.01*mm};*/
+ {0.01*mm,0.01*mm,0.01*mm,0.01*mm,0.01*mm,0.01*mm,0.01*mm,0.01*mm,
+ 0.01*mm,0.01*mm,0.01*mm,0.01*mm,0.01*mm,0.01*mm,0.01*mm};*/
G4MaterialPropertiesTable* myMPT2 = new G4MaterialPropertiesTable();
myMPT2->AddProperty("RINDEX", PMMA_PP, PMMA_RIND, NUMENTRIES);
@@ -774,30 +836,30 @@ G4Material* MaterialManager::GetMaterialFromLibrary(string Name,
//--------------------- Optical Properties ---------------------//
const G4int NUMENTRIES = 15;
G4double CsI_PP[NUMENTRIES]
- = {10 * eV, 3.5 * eV, 3.25 * eV, 3.2 * eV, 3.15 * eV,
- 3.099 * eV, 3.0 * eV, 2.95 * eV, 2.88 * eV, 2.75 * eV,
- 2.695 * eV, 2.53 * eV, 2.38 * eV, 2.30 * eV, 0 * eV};
+ = {10 * eV, 3.5 * eV, 3.25 * eV, 3.2 * eV, 3.15 * eV,
+ 3.099 * eV, 3.0 * eV, 2.95 * eV, 2.88 * eV, 2.75 * eV,
+ 2.695 * eV, 2.53 * eV, 2.38 * eV, 2.30 * eV, 0 * eV};
G4double CsI_SCINT[NUMENTRIES]
- = {0.0, 0.0, 0.1, 0.2, 0.4, 0.65, 0.8, 0.95,
- 0.82, 0.7, 0.5, 0.21, 0.05, 0, 0};
+ = {0.0, 0.0, 0.1, 0.2, 0.4, 0.65, 0.8, 0.95,
+ 0.82, 0.7, 0.5, 0.21, 0.05, 0, 0};
G4double CsI_RIND[NUMENTRIES]
- = {1.505, 1.505, 1.505, 1.505, 1.505, 1.505, 1.505, 1.505,
- 1.505, 1.505, 1.505, 1.505, 1.505, 1.505, 1.505};
+ = {1.505, 1.505, 1.505, 1.505, 1.505, 1.505, 1.505, 1.505,
+ 1.505, 1.505, 1.505, 1.505, 1.505, 1.505, 1.505};
G4double CsI_ABSL[NUMENTRIES]
- = {1.5 * m, 1.5 * m, 1.5 * m, 1.5 * m, 1.5 * m,
- 1.5 * m, 1.5 * m, 1.5 * m, 1.5 * m, 1.5 * m,
- 1.5 * m, 1.5 * m, 1.5 * m, 1.5 * m, 1.5 * m};
+ = {1.5 * m, 1.5 * m, 1.5 * m, 1.5 * m, 1.5 * m,
+ 1.5 * m, 1.5 * m, 1.5 * m, 1.5 * m, 1.5 * m,
+ 1.5 * m, 1.5 * m, 1.5 * m, 1.5 * m, 1.5 * m};
G4MaterialPropertiesTable* myMPT1 = new G4MaterialPropertiesTable();
myMPT1->AddProperty("RINDEX", CsI_PP, CsI_RIND, NUMENTRIES); /// Constant?
myMPT1->AddProperty("ABSLENGTH", CsI_PP, CsI_ABSL,
- NUMENTRIES); // Constant?
+ NUMENTRIES); // Constant?
myMPT1->AddProperty("FASTCOMPONENT", CsI_PP, CsI_SCINT,
- NUMENTRIES); // Spectrum
+ NUMENTRIES); // Spectrum
myMPT1->AddProperty("SLOWCOMPONENT", CsI_PP, CsI_SCINT,
- NUMENTRIES); // Spectrum
+ NUMENTRIES); // Spectrum
myMPT1->AddConstProperty("SCINTILLATIONYIELD", 13000. / MeV);
myMPT1->AddConstProperty("RESOLUTIONSCALE", 1.0);
@@ -813,8 +875,8 @@ G4Material* MaterialManager::GetMaterialFromLibrary(string Name,
else {
G4cout << "ERROR: Material requested \"" << Name
- << "\" is not available in the Material Library, trying with NIST"
- << G4endl;
+ << "\" is not available in the Material Library, trying with NIST"
+ << G4endl;
G4NistManager* man = G4NistManager::Instance();
return man->FindOrBuildMaterial(Name.c_str());
}
@@ -837,7 +899,7 @@ G4Element* MaterialManager::GetElementFromLibrary(string Name) {
if (!m_D) {
m_D = new G4Element(Name.c_str(), Name.c_str(), 1);
G4Isotope* isotope
- = new G4Isotope(Name.c_str(), 1, 2, 2.01410178 * g / mole);
+ = new G4Isotope(Name.c_str(), 1, 2, 2.01410178 * g / mole);
m_D->AddIsotope(isotope, 1);
}
return m_D;
@@ -847,7 +909,7 @@ G4Element* MaterialManager::GetElementFromLibrary(string Name) {
if (!m_T) {
m_T = new G4Element(Name.c_str(), Name.c_str(), 1);
G4Isotope* isotope
- = new G4Isotope(Name.c_str(), 1, 3, 3.0160492 * g / mole);
+ = new G4Isotope(Name.c_str(), 1, 3, 3.0160492 * g / mole);
m_T->AddIsotope(isotope, 1);
}
return m_T;
@@ -857,7 +919,7 @@ G4Element* MaterialManager::GetElementFromLibrary(string Name) {
if (!m_He3) {
m_He3 = new G4Element(Name.c_str(), Name.c_str(), 1);
G4Isotope* isotope
- = new G4Isotope(Name.c_str(), 2, 1, 3.0160293 * g / mole);
+ = new G4Isotope(Name.c_str(), 2, 1, 3.0160293 * g / mole);
m_He3->AddIsotope(isotope, 1);
}
return m_He3;
@@ -870,95 +932,104 @@ G4Element* MaterialManager::GetElementFromLibrary(string Name) {
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//
G4Material* MaterialManager::GetGasFromLibrary(string Name, double Pressure, double Temperature){
- ostringstream oss;
- oss << Name<< "_"<<Pressure<<"_"<<Temperature;
- string newName= oss.str();
- map<string,G4Material*>::iterator it;
- it = m_Material.find(Name);
-
- double density = 0 ;
-
- G4double Vm=0.08206*Temperature*atmosphere/(Pressure*kelvin);
-
- // The element is not found
- if(it==m_Material.end()){
- if(Name == "CF4"){ // 52 torr
- density = 3.72*kg/m3;
- double refTemp= (273.15+15)*kelvin;
- double refPres= 1.01325*bar;
- density = density*(refTemp/Temperature)/(refPres/Pressure);
- G4Material* material = new G4Material("NPS_"+newName,density,2,kStateGas,Temperature,Pressure);
- material->AddElement(GetElementFromLibrary("C"), 1);
- material->AddElement(GetElementFromLibrary("F"), 4);
- m_Material[newName]=material;
- return material;
- }
-
- if(Name == "He"){
- density = (4.0026/Vm)*mg/cm3;
- G4Material* material = new G4Material("NPS_"+newName,density,1,kStateGas,Temperature,Pressure);
- material->AddElement(GetElementFromLibrary("He"), 1);
- m_Material[newName]=material;
- return material;
- }
-
- if(Name == "iC4H10" || Name == "Isobutane" || Name == "isobutane"){
- density = ((4*12.0107+10*1.00794)/Vm)*mg/cm3;
- G4Material* material = new G4Material("NPS_"+newName,density,2,kStateGas,Temperature,Pressure);
- material->AddElement(GetElementFromLibrary("C"), 4);
- material->AddElement(GetElementFromLibrary("H"), 10);
- m_Material[newName]=material;
- return material;
- }
-
- if(Name == "CH4"){
- density = ((12.0107+4*1.00794)/Vm)*mg/cm3;
- G4Material* material = new G4Material("NPS_"+newName,density,2,kStateGas,Temperature,Pressure);
- material->AddElement(GetElementFromLibrary("C"), 1);
- material->AddElement(GetElementFromLibrary("H"), 4);
- m_Material[newName]=material;
- return material;
- }
-
- if(Name == "CO2"){
- density = ((12.0107+2*16)/Vm)*mg/cm3;
- G4Material* material = new G4Material("NPS_"+newName,density,2,kStateGas,Temperature,Pressure);
- material->AddElement(GetElementFromLibrary("C"), 1);
- material->AddElement(GetElementFromLibrary("O"), 2);
- m_Material[newName]=material;
- return material;
- }
-
- if(Name == "H2"){
- density = (2*1.00794/Vm)*mg/cm3;
- G4Material* material = new G4Material("NPS_"+newName,density,1,kStateGas,Temperature,Pressure);
- material->AddElement(GetElementFromLibrary("H"), 2);
- //material->AddElement(GetElementFromLibrary("H"), 1);
- m_Material[newName]=material;
- return material;
- }
-
- if(Name == "D2"){
- density = (2*2.0140/Vm)*mg/cm3;
- G4Material* material = new G4Material("NPS_"+newName,density,1,kStateGas,Temperature,Pressure);
- material->AddElement(GetElementFromLibrary("D"), 2);
- //material->AddElement(GetElementFromLibrary("D"), 1);
- m_Material[newName]=material;
- return material;
- }
-
-
- else{
- exit(1);
- }
- }
+ ostringstream oss;
+ oss << Name<< "_"<<Pressure<<"_"<<Temperature;
+ string newName= oss.str();
+ map<string,G4Material*>::iterator it;
+ it = m_Material.find(Name);
+
+ double density = 0 ;
+
+ G4double Vm=0.08206*Temperature*atmosphere/(Pressure*kelvin);
+
+ // The element is not found
+ if(it==m_Material.end()){
+ if(Name == "CF4"){ // 52 torr
+ density = 3.72*kg/m3;
+ double refTemp= (273.15+15)*kelvin;
+ double refPres= 1.01325*bar;
+ density = density*(refTemp/Temperature)/(refPres/Pressure);
+ G4Material* material = new G4Material("NPS_"+newName,density,2,kStateGas,Temperature,Pressure);
+ material->AddElement(GetElementFromLibrary("C"), 1);
+ material->AddElement(GetElementFromLibrary("F"), 4);
+ m_Material[newName]=material;
+ return material;
+ }
+
+ if(Name == "He"){
+ density = (4.0026/Vm)*mg/cm3;
+ G4Material* material = new G4Material("NPS_"+newName,density,1,kStateGas,Temperature,Pressure);
+ material->AddElement(GetElementFromLibrary("He"), 1);
+ m_Material[newName]=material;
+ return material;
+ }
+
+ if(Name == "iC4H10" || Name == "Isobutane" || Name == "isobutane"){
+ density = ((4*12.0107+10*1.00794)/Vm)*mg/cm3;
+ G4Material* material = new G4Material("NPS_"+newName,density,2,kStateGas,Temperature,Pressure);
+ material->AddElement(GetElementFromLibrary("C"), 4);
+ material->AddElement(GetElementFromLibrary("H"), 10);
+ m_Material[newName]=material;
+ return material;
+ }
+
+ if(Name == "CH4"){
+ density = ((12.0107+4*1.00794)/Vm)*mg/cm3;
+ G4Material* material = new G4Material("NPS_"+newName,density,2,kStateGas,Temperature,Pressure);
+ material->AddElement(GetElementFromLibrary("C"), 1);
+ material->AddElement(GetElementFromLibrary("H"), 4);
+ m_Material[newName]=material;
+ return material;
+ }
+
+ if(Name == "CO2"){
+ density = ((12.0107+2*16)/Vm)*mg/cm3;
+ G4Material* material = new G4Material("NPS_"+newName,density,2,kStateGas,Temperature,Pressure);
+ material->AddElement(GetElementFromLibrary("C"), 1);
+ material->AddElement(GetElementFromLibrary("O"), 2);
+ m_Material[newName]=material;
+ return material;
+ }
+
+ if(Name == "H2"){
+ density = (2*1.00794/Vm)*mg/cm3;
+ G4Material* material = new G4Material("NPS_"+newName,density,1,kStateGas,Temperature,Pressure);
+ material->AddElement(GetElementFromLibrary("H"), 2);
+ //material->AddElement(GetElementFromLibrary("H"), 1);
+ m_Material[newName]=material;
+ return material;
+ }
+
+ if(Name == "D2"){
+ density = (2*2.0140/Vm)*mg/cm3;
+ G4Material* material = new G4Material("NPS_"+newName,density,1,kStateGas,Temperature,Pressure);
+ material->AddElement(GetElementFromLibrary("D"), 2);
+ //material->AddElement(GetElementFromLibrary("D"), 1);
+ m_Material[newName]=material;
+ return material;
+ }
+
+ /* if(Name == "mixMINOS"){ */
+ /* density = (2*2.0140/Vm)*mg/cm3; */
+ /* G4Material* material = new G4Material("NPS_"+newName,density,1,kStateGas,Temperature,Pressure); */
+ /* material->AddElement(GetElementFromLibrary("D"), 2); */
+ /* //material->AddElement(GetElementFromLibrary("D"), 1); */
+ /* m_Material[newName]=material; */
+ /* return material; */
+ /* } */
+
+
+ else{
+ exit(1);
+ }
+ }
return NULL;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Generate a DEDX file table using the material used in the geometry
void MaterialManager::WriteDEDXTable(G4ParticleDefinition* Particle,
- G4double Emin, G4double Emax) {
+ G4double Emin, G4double Emax) {
map<string, G4Material*>::iterator it;
if (Particle->GetPDGCharge() == 0)
return;
@@ -969,7 +1040,7 @@ void MaterialManager::WriteDEDXTable(G4ParticleDefinition* Particle,
// Remove NPS name
Name.erase(0, 4);
G4String Path = GlobalPath + "/Inputs/EnergyLoss/"
- + Particle->GetParticleName() + "_" + Name + ".G4table";
+ + Particle->GetParticleName() + "_" + Name + ".G4table";
ofstream File;
File.open(Path);
@@ -978,8 +1049,8 @@ void MaterialManager::WriteDEDXTable(G4ParticleDefinition* Particle,
return;
File << "Table from Geant4 generate using NPSimulation \t"
- << "Particle: " << Particle->GetParticleName()
- << "\tMaterial: " << it->second->GetName() << G4endl;
+ << "Particle: " << Particle->GetParticleName()
+ << "\tMaterial: " << it->second->GetName() << G4endl;
// G4cout << Particle->GetParticleName() << "\tMaterial: " <<
// it->second->GetName() <<endl;
@@ -993,7 +1064,7 @@ void MaterialManager::WriteDEDXTable(G4ParticleDefinition* Particle,
for (G4double E = Emin; E < Emax; E += step) {
dedx = emCalculator.ComputeTotalDEDX(E, Particle, it->second)
- / (MeV / micrometer);
+ / (MeV / micrometer);
if (before) {
if (abs(before - dedx) / abs(before) < 0.01)
step *= 2;
@@ -1009,11 +1080,11 @@ void MaterialManager::WriteDEDXTable(G4ParticleDefinition* Particle,
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Generate a DEDX file table using the material used in the geometry
void MaterialManager::WriteDEDXTable(std::set<string> Particle, G4double Emin,
- G4double Emax) {
+ G4double Emax) {
std::set<string>::iterator it;
for (it = Particle.begin(); it != Particle.end(); it++) {
G4ParticleDefinition* p
- = G4ParticleTable::GetParticleTable()->FindParticle((*it));
+ = G4ParticleTable::GetParticleTable()->FindParticle((*it));
WriteDEDXTable(p, Emin, Emax);
}
}
@@ -1024,18 +1095,18 @@ void MaterialManager::CreateSampleVolumes(G4LogicalVolume* world_log) {
G4double SampleSize = 1 * um;
G4double WorldSize = 10.0 * m;
G4Box* sample_box
- = new G4Box("sample_box", SampleSize, SampleSize, SampleSize);
+ = new G4Box("sample_box", SampleSize, SampleSize, SampleSize);
G4int i = 1;
G4double Coord1 = WorldSize - SampleSize;
G4double Coord2 = 0;
map<string, G4Material*>::iterator it;
for (it = m_Material.begin(); it != m_Material.end(); it++) {
G4LogicalVolume* sample_log
- = new G4LogicalVolume(sample_box, it->second, "sample_log", 0, 0, 0);
+ = new G4LogicalVolume(sample_box, it->second, "sample_log", 0, 0, 0);
sample_log->SetVisAttributes(G4VisAttributes::Invisible);
Coord2 = WorldSize - i * 2 * SampleSize;
i++;
new G4PVPlacement(0, G4ThreeVector(Coord1, Coord2, -Coord1), sample_log,
- "sample", world_log, false, 0);
+ "sample", world_log, false, 0);
}
}
diff --git a/NPSimulation/Detectors/Dali/Dali.cc b/NPSimulation/Detectors/Dali/Dali.cc
index 84c1eb247..c8848e871 100644
--- a/NPSimulation/Detectors/Dali/Dali.cc
+++ b/NPSimulation/Detectors/Dali/Dali.cc
@@ -73,13 +73,13 @@ namespace Dali_NS{
// Energy and time Resolution
const double EnergyThreshold = 0*MeV;
const double ResoTime = 0.0*ns; //4.5*ns ;
- const double ResoEnergy = 1.36*MeV ; // mean Resolution(FWHM) 1.7% of 80MeV from slides 20170214-SAMURAI34-setup-DALI.pdf if 1.7% of 30MeV = 0.51 MeV // 0.001*MeV ;
+ const double ResoEnergy = 0.122; //Relative resolution DeltaE = 0.122*Sqrt(E)
+ /* const double ResoEnergy = 1.36*MeV ; // mean Resolution(FWHM) 1.7% of 80MeV from slides 20170214-SAMURAI34-setup-DALI.pdf if 1.7% of 30MeV = 0.51 MeV // 0.001*MeV ; */
const double Radius = 50*mm ;
const double Width = 49.76*mm ;
const double Hight = 84.81*mm ;
const double Thickness = 164.82*mm ;
const double LengthPMT = 152.62*mm ;
- const string Material = "NaI";
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -87,287 +87,194 @@ namespace Dali_NS{
// Dali Specific Method
Dali::Dali(){
m_Event = new TDaliData() ;
- m_DaliScorer = 0;
m_SquareDetector = 0;
- m_SquareDetector_Can = 0;
- m_SquareDetector_CanMgO =0;
- m_CylindricalDetector = 0;
- m_SquareDetector_Crystal = 0;
- Logic_ArrayDali_1 =0;
-
- // RGB Color + Transparency
- m_VisSquare = new G4VisAttributes(G4Colour(0, 1, 1/*, 0.3*/));
- m_VisCylinder = new G4VisAttributes(G4Colour(0, 0, 1/*, 0.3*/));
-
-}
+ m_DaliScorer = 0;
+ Log_Dali_1Volume = 0;
+ Log_Al_Cryst_can = 0;
+ Log_MgO_Cryst_can =0;
+ Log_Crystal = 0;
+ Log_Dali_3Volume=0;
+ }
Dali::~Dali(){
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void Dali::AddDetector(G4ThreeVector POS, string Shape){
+void Dali::AddDetector(G4ThreeVector POS ){
// Convert the POS value to R theta Phi as Cylindrical coordinate is easier in G4
m_R.push_back(POS.perp());
m_Alpha.push_back(POS.phi());
m_Zeta.push_back(POS.y());
- m_Shape.push_back(Shape);
}
-
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void Dali::AddDetector(double R, double Theta, double Phi, string Shape){
-
+void Dali::AddDetector(double R, double Theta, double Phi){
double m_r, m_alpha, m_zeta;
-
m_r = R*cos(Phi);
m_alpha = Theta;
m_zeta = R*sin(Phi);
-
m_R.push_back(m_r);
m_Alpha.push_back(m_alpha);
m_Zeta.push_back(m_zeta);
- m_Shape.push_back(Shape);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void Dali::AddDetector2(double R, double Alpha, double Zeta, string Shape){
+void Dali::AddDetector(double R, double Alpha, double Zeta, int Ring){
m_R.push_back(R);
m_Alpha.push_back(Alpha);
m_Zeta.push_back(Zeta);
- m_Shape.push_back(Shape);
+ m_Ring.push_back(Ring);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Definition Materials MgO and NaI(Tl)
void Dali::DefinitionMaterials()
{
- //G4Element* H = new G4Element("Hydrogen","H" , 1., 1.01*g/mole);
-
- G4Isotope* Mg24 = new G4Isotope ("Mg24", 12, 24, 23.985041*g/mole);
- G4Isotope* Mg25 = new G4Isotope ("Mg25", 12, 25, 24.985836*g/mole);
- G4Isotope* Mg26 = new G4Isotope ("Mg26", 12, 26, 25.982592*g/mole);
- G4Element* Mg= new G4Element("elMagnesium","Mg",3);
- Mg->AddIsotope(Mg24, 78.99*perCent);
- Mg->AddIsotope(Mg25, 10*perCent);
- Mg->AddIsotope(Mg26, 11.01*perCent);
-
- G4Isotope* O16 = new G4Isotope ("O16", 8, 16, 15.99*g/mole);
- G4Isotope* O17 = new G4Isotope ("O17", 8, 17, 17.00*g/mole);
- G4Isotope* O18 = new G4Isotope ("O18", 8, 18, 18.00*g/mole);
- G4Element* O = new G4Element("elOxygen","O",3);
- O->AddIsotope(O16, 99.76*perCent);
- O->AddIsotope(O17, 0.04*perCent);
- O->AddIsotope(O18, 0.20*perCent);
+
+ G4Element *Tl = new G4Element("Thallium","Tl",81.,204.383*g/mole );
-
- MgO = new G4Material("MgO",3.6*g/cm3, 2 );
- MgO->AddElement(Mg,1);
- MgO->AddElement(O, 1);
-
- G4Element *elTl = new G4Element("Thallium","Tl",81.,204.383*g/mole );
-
NaI_Tl = new G4Material("NaI_Tl",3.6667*g/cm3, 2);
NaI_Tl->AddMaterial(MaterialManager::getInstance()->GetMaterialFromLibrary("NaI"), 99.6*perCent);
- NaI_Tl->AddElement(elTl, 0.4*perCent);
-
+ NaI_Tl->AddElement(Tl, 0.4*perCent);
+
}
-
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
G4LogicalVolume* Dali::BuildSquareDetector(){
if(!m_SquareDetector){
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-
- G4Box* box_3can = new G4Box("Dali_3BoxCan", Dali_NS::Hight*0.5,
- Dali_NS::Width*0.5*3, Dali_NS::Thickness*0.5 + Dali_NS::LengthPMT/2.+11.5/2.*mm /*last part is PMTVolume*/ );
- G4Material* Aria = MaterialManager::getInstance()->GetMaterialFromLibrary("Air");
- Logic_ArrayDali_1 = new G4LogicalVolume(box_3can,Aria,"logic_ArrayDali",0,0,0);
-
- Logic_ArrayDali_1->SetVisAttributes(G4VisAttributes(G4Colour(1,1,1, 0)));
-
- G4Box* box_can = new G4Box("Dali_BoxCan", Dali_NS::Hight*0.5,
- Dali_NS::Width*0.5, Dali_NS::Thickness*0.5);
-
- G4Box* box_canandPMT = new G4Box("Dali_BoxCan", Dali_NS::Hight*0.5,
- Dali_NS::Width*0.5, Dali_NS::Thickness*0.5 + Dali_NS::LengthPMT/2.+11.5/2.*mm /*last part is PMTVolume*/ );
-
std::vector<G4TwoVector> polygon;
- polygon.push_back(G4TwoVector(Dali_NS::Hight*0.5, Dali_NS::Width*0.5*3. ) ) ;
- polygon.push_back(G4TwoVector(Dali_NS::Hight*0.5, -Dali_NS::Width*0.5*3. ) ) ;
- polygon.push_back(G4TwoVector(-Dali_NS::Hight*0.5, -Dali_NS::Width*0.5*3. ) ) ;
- polygon.push_back(G4TwoVector(-Dali_NS::Hight*0.5, Dali_NS::Width*0.5*3. ) ) ;
-
- // std::vector<ZSection> zsection;
- // zsection.push_back(ZSection (Dali_NS::Thickness*0.5, {0,0}, 1. ) );
- // zsection.push_back(ZSection (-Dali_NS::Thickness*0.5-19.5*2.*mm , {0,0}, 1. ) );
-
+ polygon.push_back(G4TwoVector(Dali_NS::Hight*0.5, Dali_NS::Width*0.5*3.) ) ;
+ polygon.push_back(G4TwoVector(Dali_NS::Hight*0.5, -Dali_NS::Width*0.5*3.) ) ;
+ polygon.push_back(G4TwoVector(-Dali_NS::Hight*0.5, -Dali_NS::Width*0.5*3.) ) ;
+ polygon.push_back(G4TwoVector(-Dali_NS::Hight*0.5, Dali_NS::Width*0.5*3.) ) ;
- G4Box* Extrudedbox_can = new G4Box("Dali_BoxCan", Dali_NS::Hight*0.5,Dali_NS::Width*0.5, Dali_NS::LengthPMT/2.+11.5/2.*mm);
-
+ G4Material* Aria = MaterialManager::getInstance()->GetMaterialFromLibrary("Air");
+ G4Material* Al = MaterialManager::getInstance()->GetMaterialFromLibrary("Al");
+ G4Material* MgO = MaterialManager::getInstance()->GetMaterialFromLibrary("MgO");
+ G4Material* muMetal = MaterialManager::getInstance()->GetMaterialFromLibrary("mumetal");
+ G4Material* Vacuum = MaterialManager::getInstance()->GetMaterialFromLibrary("Vacuum");
+ G4Material* BoroSili_Glass = MaterialManager::getInstance()->GetMaterialFromLibrary("Borosillicate_Glass");
+
+ ///// CONTENAIR VOLUMES
+ G4Box* Dali_3Volume = new G4Box("Dali_3Volume",
+ Dali_NS::Hight*0.5,
+ Dali_NS::Width*0.5*3,
+ Dali_NS::Thickness*0.5 + Dali_NS::LengthPMT/2.+11.5/2.*mm );
+ Log_Dali_3Volume = new G4LogicalVolume(Dali_3Volume,
+ Aria,"log_Dali_3Volume",0,0,0);
+ G4Box* Dali_1Volume = new G4Box("Dali_1Volume",
+ Dali_NS::Hight*0.5,
+ Dali_NS::Width*0.5,
+ Dali_NS::Thickness*0.5 + Dali_NS::LengthPMT/2.+11.5/2.*mm );
+ Log_Dali_1Volume = new G4LogicalVolume(Dali_1Volume,
+ Aria,"logic_1DaliVolume",0,0,0);
+
+ G4Box* Extrudedbox_can = new G4Box("extrude_box", Dali_NS::Hight*0.5,
+ Dali_NS::Width*0.5,
+ Dali_NS::LengthPMT/2.+11.5/2.*mm);
AriaExtrude = new G4LogicalVolume(Extrudedbox_can,Aria, "logic_Ariaextrude",0,0,0);
- G4Material* DetectorCanMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary("Al");
- m_SquareDetector_Can = new G4LogicalVolume(box_can,DetectorCanMaterial,"logic_Dali_Can",0,0,0);
- m_Square2Detector_Can = new G4LogicalVolume(box_canandPMT, Aria,"logic_Dali_CanandPMT",0,0,0);
-
- //THE PMT
+ /////// PMTs objets
G4Tubs* AlPMT = new G4Tubs("AlPMT",16.5*mm, 19.5*mm,Dali_NS::LengthPMT/2.,0*deg,360*deg);
+ lAlPMT = new G4LogicalVolume(AlPMT, Vacuum ,"lAlPMT",0,0,0);
+
G4Tubs* MuPMT = new G4Tubs("MuPMT",16.5*mm,20.*mm,Dali_NS::LengthPMT/2.,0*deg,360*deg);
+ lMuPMT = new G4LogicalVolume(MuPMT, Vacuum ,"lMuPMT",0,0,0);
+
G4Box* TopPlatePMT = new G4Box("TopPlatePMT", Dali_NS::Hight*0.5-1*mm,
- Dali_NS::Width*0.5-1*mm, 11.5/2.*mm );
+ Dali_NS::Width*0.5-1*mm,
+ 11.5/2.*mm );
+ lTopPlatePMT = new G4LogicalVolume(TopPlatePMT, Vacuum ,"lTopPlatePMT",0,0,0);
+
G4Tubs* GlassPMT = new G4Tubs("GlassPMT", 0. , 16.5*mm , 11.5/2.*mm ,0*deg,360*deg);
+ lGlassPMT = new G4LogicalVolume(GlassPMT , Vacuum ,"lGlassPMT",0,0,0); // Cyril
- lAlPMT = new G4LogicalVolume(AlPMT, DetectorCanMaterial ,"lAlPMT",0,0,0);
- lMuPMT = new G4LogicalVolume(MuPMT, DetectorCanMaterial ,"lMuPMT",0,0,0);
- lTopPlatePMT = new G4LogicalVolume(TopPlatePMT, DetectorCanMaterial ,"lTopPlatePMT",0,0,0);
- lGlassPMT = new G4LogicalVolume(GlassPMT , MaterialManager::getInstance()->GetMaterialFromLibrary("Borosillicate_Glass") ,"lGlassPMT",0,0,0);
-
-
- G4VisAttributes* Can_Attributes = new G4VisAttributes(G4Colour(0.5,0.5,0.5, .3));
- m_SquareDetector_Can->SetVisAttributes(Can_Attributes);
- m_Square2Detector_Can->SetVisAttributes(G4VisAttributes(G4Colour(1,1,1,0)));
-
- AriaExtrude->SetVisAttributes(G4VisAttributes(G4Colour(1,1,1,0)));
- lAlPMT->SetVisAttributes(Can_Attributes);
- lMuPMT->SetVisAttributes(Can_Attributes);
- lTopPlatePMT->SetVisAttributes(Can_Attributes);
-
- G4Box* box_MgO = new G4Box("Dali_BoxMgO", Dali_NS::Hight*0.5-1*mm,
- Dali_NS::Width*0.5-1*mm, Dali_NS::Thickness*0.5-1*mm); // Size of Al Can but w/o thickness of AlCan
-
- m_SquareDetector_CanMgO = new G4LogicalVolume(box_MgO,MgO,"logic_Dali_CanMg0",0,0,0);
- G4Box* box_crystal = new G4Box("Dali_BoxNaI", Dali_NS::Hight*0.5-2.4*mm,
- Dali_NS::Width*0.5-2.4*mm, Dali_NS::Thickness*0.5-2.4*mm); // Size of AlCan but w/o thickness of AlCan and MgO
-// G4Material* DetectorMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary(Dali_NS::Material);
+ //////////////// CRYSTAL part
+ G4Box* Al_Cryst_can = new G4Box("Al_Cryst_can",
+ Dali_NS::Hight*0.5,
+ Dali_NS::Width*0.5,
+ Dali_NS::Thickness*0.5);
+ Log_Al_Cryst_can = new G4LogicalVolume(Al_Cryst_can,
+ Al,"log_Al_Cryst_can",0,0,0);
+
+ G4Box* MgO_Cryst_can = new G4Box("MgO_Cryst_can",
+ Dali_NS::Hight*0.5-1*mm,
+ Dali_NS::Width*0.5-1*mm,
+ Dali_NS::Thickness*0.5-1*mm);
+ Log_MgO_Cryst_can = new G4LogicalVolume(MgO_Cryst_can,
+ MgO, "logic_Dali_CanMg0",0,0,0);
+
+ G4Box* Crystal = new G4Box("Crystal",
+ Dali_NS::Hight*0.5-2.4*mm,
+ Dali_NS::Width*0.5-2.4*mm,
+ Dali_NS::Thickness*0.5-2.4*mm);
+ Log_Crystal = new G4LogicalVolume(Crystal,
+ NaI_Tl,"logic_Dali_Box",0,0,0);
- m_SquareDetector_Crystal = new G4LogicalVolume(box_crystal,
- NaI_Tl,
- // DetectorMaterial,
- "logic_Dali_Box",0,0,0);
-
G4ThreeVector positionnull = G4ThreeVector(0,0,0);
-
- // PMT Volume -
+ // PMT part -
new G4PVPlacement(0, positionnull,
- lAlPMT ,
- "AlPMT",
- lMuPMT,
- false,
- 0);
+ lAlPMT ,"AlPMT",lMuPMT,false,0);
new G4PVPlacement(0, G4ThreeVector(0,0, -11.5/2.*mm ),
- lMuPMT ,
- "MuPMT",
- AriaExtrude,
- false,
- 0);
-
-
- new G4PVPlacement(0, positionnull,
- lGlassPMT,
- "GlassPMT",
- lTopPlatePMT,
- false,
- 0);
+ lMuPMT ,"MuPMT",AriaExtrude,false,0);
+ new G4PVPlacement(0, positionnull,lGlassPMT,"GlassPMT",
+ lTopPlatePMT,false,0);
new G4PVPlacement(0, G4ThreeVector(0,0, Dali_NS::LengthPMT/2. ),
- lTopPlatePMT,
- "TopPlatePMT",
- AriaExtrude,
- false,
- 0);
-
-
-
+ lTopPlatePMT,"TopPlatePMT",AriaExtrude,false,0);
new G4PVPlacement(0, G4ThreeVector(0,0, -Dali_NS::Thickness*0.5 ),
- AriaExtrude,
- "PMTVolume",
- m_Square2Detector_Can,
- false,
- 0);
-
-
-
+ AriaExtrude,"PMTVolume",Log_Dali_1Volume,false,0);
+
- new G4PVPlacement(0, G4ThreeVector(0,0, Dali_NS::LengthPMT/2.+11.5/2.*mm ),
- m_SquareDetector_Can,
+ // Cryst Part -
+ new G4PVPlacement(0, G4ThreeVector(0,0,
+ Dali_NS::LengthPMT/2.+11.5/2.*mm ),
+ Log_Al_Cryst_can,
"DetectorVolume",
- m_Square2Detector_Can,
- false,
- 0);
-
-
- // MgO Volume -
+ Log_Dali_1Volume,false,0);
+
new G4PVPlacement(0, positionnull,
- m_SquareDetector_CanMgO,
- "MgO",
- m_SquareDetector_Can,
- false,
- 0);
- G4VisAttributes* MgO_Attributes = new G4VisAttributes(G4Colour(1,1,1, .3));
- m_SquareDetector_CanMgO->SetVisAttributes(MgO_Attributes);
-
-
- // NaI Volume -
+ Log_MgO_Cryst_can,
+ "MgO_Can",
+ Log_Al_Cryst_can,false,0);
new G4PVPlacement(0, positionnull,
- m_SquareDetector_Crystal,
- "CrystalNaI",
- m_SquareDetector_CanMgO,
- false,
- 0);
-
- m_SquareDetector_Crystal->SetVisAttributes(m_VisSquare);
- m_SquareDetector_Crystal->SetSensitiveDetector(m_DaliScorer);
-
-
- // new G4PVPlacement(0, positionnull,
- // m_SquareDetector_Crystal,
- // "CrystalNaI",
- // m_SquareDetector_CanMgO,
- // false,
- // 0);
-
-
-
+ Log_Crystal,
+ "CrystalNaI",
+ Log_MgO_Cryst_can,false,0);
new G4PVReplica("DaliArrayElement",
- m_Square2Detector_Can,
- Logic_ArrayDali_1 ,
- kYAxis,
- 3,
- Dali_NS::Width, //?????????
- 0);
-
-
-
- }
-
- return Logic_ArrayDali_1;
-}
+ Log_Dali_1Volume,
+ Log_Dali_3Volume ,
+ kYAxis,3,Dali_NS::Width,0);
+
+ G4VisAttributes* MgO_Color = new G4VisAttributes(G4Colour(1,1,1, .4));
+ G4VisAttributes* Al_Color = new G4VisAttributes(G4Colour(0.5,0.5,0.5, .3));
+ G4VisAttributes* Crystal_Color = new G4VisAttributes(G4Colour(0, 1, 1));
+ G4VisAttributes* mumetal_Color = new G4VisAttributes(G4Colour(0, 0.5, 1, .3));
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-G4LogicalVolume* Dali::BuildCylindricalDetector(){
- if(!m_CylindricalDetector){
- G4Tubs* tub = new G4Tubs("Dali_Cyl",0,Dali_NS::Radius, Dali_NS::Thickness*0.5,0,360*deg);
- G4Material* DetectorMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary(Dali_NS::Material);
- m_CylindricalDetector = new G4LogicalVolume(tub,DetectorMaterial,"logic_Dali_tub",0,0,0);
- m_CylindricalDetector->SetVisAttributes(m_VisSquare);
- m_CylindricalDetector->SetSensitiveDetector(m_DaliScorer);
+ Log_Dali_3Volume->SetVisAttributes(G4VisAttributes(G4Colour(1,1,1, 0)));
+ Log_Dali_1Volume->SetVisAttributes(G4VisAttributes(G4Colour(1,1,1,0)));
+ AriaExtrude->SetVisAttributes(G4VisAttributes(G4Colour(1,1,1,0)));
+
+ Log_MgO_Cryst_can->SetVisAttributes(MgO_Color);
+ Log_Crystal->SetVisAttributes(Crystal_Color);
+ Log_Al_Cryst_can->SetVisAttributes(Al_Color);
+
+ lAlPMT->SetVisAttributes(Al_Color);
+ lMuPMT->SetVisAttributes(Al_Color);
+ lTopPlatePMT->SetVisAttributes(Al_Color);
+ Log_Crystal->SetSensitiveDetector(m_DaliScorer);
}
- return m_CylindricalDetector;
-}
+ return Log_Dali_3Volume;
+} // end BuildSquareDetector
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -381,36 +288,25 @@ void Dali::ReadConfiguration(NPL::InputParser parser){
if(NPOptionManager::getInstance()->GetVerboseLevel())
cout << "//// " << blocks.size() << " detectors found " << endl;
- vector<string> cart = {"POS","Shape"};
- vector<string> sphe = {"R","Theta","Phi","Shape"};
- vector<string> cyli = {"R","Alpha","Zeta","Shape"};
+ vector<string> cart = {"POS"};
+ vector<string> cyli = {"R","Alpha","Zeta"};
for(unsigned int i = 0 ; i < blocks.size() ; i++){
if(blocks[i]->HasTokenList(cart)){
if(NPOptionManager::getInstance()->GetVerboseLevel())
cout << endl << "//// Dali " << i+1 << endl;
-
G4ThreeVector Pos = NPS::ConvertVector(blocks[i]->GetTVector3("POS","mm"));
- string Shape = blocks[i]->GetString("Shape");
- AddDetector(Pos,Shape);
- }
- else if(blocks[i]->HasTokenList(sphe)){
- if(NPOptionManager::getInstance()->GetVerboseLevel())
- cout << endl << "//// Dali " << i+1 << endl;
- double R = blocks[i]->GetDouble("R","mm");
- double Theta = blocks[i]->GetDouble("Theta","deg");
- double Phi = blocks[i]->GetDouble("Phi","deg");
- string Shape = blocks[i]->GetString("Shape");
- AddDetector(R,Theta,Phi,Shape);
+ AddDetector(Pos);
}
+
else if(blocks[i]->HasTokenList(cyli)){
if(NPOptionManager::getInstance()->GetVerboseLevel())
cout << endl << "//// Dali " << i+1 << endl;
double R = blocks[i]->GetDouble("R","mm");
double Alpha = blocks[i]->GetDouble("Alpha","deg");
double Zeta = blocks[i]->GetDouble("Zeta","mm");
- string Shape = blocks[i]->GetString("Shape");
- AddDetector2(R,Alpha,Zeta,Shape);
+ int Ring = blocks[i]->GetInt("Ring");
+ AddDetector(R,Alpha,Zeta,Ring);
}
else{
cout << "ERROR: check your input file formatting " << endl;
@@ -426,44 +322,25 @@ void Dali::ReadConfiguration(NPL::InputParser parser){
// Called After DetecorConstruction::AddDetector Method
void Dali::ConstructDetector(G4LogicalVolume* world){
-
DefinitionMaterials();
-
for (unsigned short i = 0 ; i < m_R.size() ; i++) {
-
G4double wX = m_R[i] * cos(m_Alpha[i] ) ;
G4double wY = m_R[i] * sin(m_Alpha[i] ) ;
G4double wZ = m_Zeta[i];
- if(m_Zeta[i]<0) wZ = wZ - Dali_NS::LengthPMT/2.+11.5/2.*mm;
- else wZ = wZ + Dali_NS::LengthPMT/2.+11.5/2.*mm;
+ if(m_Ring[i]==1) wZ = wZ - (Dali_NS::Thickness+11.5*mm+Dali_NS::LengthPMT)/2. + Dali_NS::Thickness/2.;
+ else wZ = wZ + (Dali_NS::Thickness+11.5*mm+Dali_NS::LengthPMT)/2. - Dali_NS::Thickness/2.;
G4ThreeVector Det_pos = G4ThreeVector(wX, wY, wZ) ;
-
-
G4RotationMatrix* Rot = new G4RotationMatrix();
-
-
-
Rot->rotateX(180*deg);
- if(m_Zeta[i]<0){
+ if(m_Ring[i]==1){
Rot->rotateY(180*deg); Rot->rotateZ(m_Alpha[i]);
} else{Rot->rotateZ(m_Alpha[i]);}
-
- if(m_Shape[i] == "Cylindrical"){
- new G4PVPlacement(G4Transform3D(*Rot,Det_pos),
- BuildCylindricalDetector(),
- "Dali",world,false,i+1);
- }
-
- else if(m_Shape[i] == "Square"){
- new G4PVPlacement(G4Transform3D(*Rot,Det_pos),
- BuildSquareDetector(),
- "Dali",world,false,i+1);
-
-
- }
+ new G4PVPlacement(G4Transform3D(*Rot,Det_pos),
+ BuildSquareDetector(),
+ "Dali",world,false,i+1);
}
}
@@ -483,20 +360,17 @@ void Dali::InitializeRootOutput(){
// Read sensitive part and fill the Root tree.
// Called at in the EventAction::EndOfEventAvtion
void Dali::ReadSensitive(const G4Event* ){
- m_Event->Clear();
+ m_Event->Clear();
///////////
-
// Calorimeter scorer
CalorimeterScorers::PS_Calorimeter* Scorer= (CalorimeterScorers::PS_Calorimeter*) m_DaliScorer->GetPrimitive(0);
//cout << m_DaliScorer->GetNumberOfPrimitives()<<endl;
-
-
unsigned int size = Scorer->GetMult();
for(unsigned int i = 0 ; i < size ; i++){
vector<unsigned int> level = Scorer->GetLevel(i);
- double Energy = RandGauss::shoot(Scorer->GetEnergy(i),Dali_NS::ResoEnergy);
+ double Energy = RandGauss::shoot(Scorer->GetEnergy(i),Dali_NS::ResoEnergy*std::sqrt(Scorer->GetEnergy(i)));
// cout << Energy << endl;
if(Energy>Dali_NS::EnergyThreshold){
double Time = RandGauss::shoot(Scorer->GetTime(i),Dali_NS::ResoTime);
@@ -505,6 +379,7 @@ void Dali::ReadSensitive(const G4Event* ){
int DetectorNbr = (ArrayNbr-1)*3+DetectinsArrayNbr;
m_Event->SetEnergy(DetectorNbr,Energy);
m_Event->SetTime(DetectorNbr,Time);
+ /* m_Event->SetParticleID(Scorer->GetParticleID(i)); */
}
}
}
@@ -517,14 +392,13 @@ void Dali::InitializeScorers() {
vector<int> NestingLevel;
NestingLevel.push_back(3);
NestingLevel.push_back(4);
-
+
m_DaliScorer = CheckScorer("DaliScorer",already_exist) ;
if(already_exist) //Necessary?
return ; //Necessary?
-
// Otherwise the scorer is initialised
-// vector<int> level; level.push_back(0);
+ // vector<int> level; level.push_back(0);
G4VPrimitiveScorer* Calorimeter= new CalorimeterScorers::PS_Calorimeter("Calorimeter", NestingLevel) ;
G4VPrimitiveScorer* Interaction= new InteractionScorers::PS_Interactions("Interaction",ms_InterCoord, 0) ;
//and register it to the multifunctionnal detector
diff --git a/NPSimulation/Detectors/Dali/Dali.hh b/NPSimulation/Detectors/Dali/Dali.hh
index 02fb8c6dd..9b265abce 100644
--- a/NPSimulation/Detectors/Dali/Dali.hh
+++ b/NPSimulation/Detectors/Dali/Dali.hh
@@ -51,35 +51,33 @@ class Dali : public NPS::VDetector{
////////////////////////////////////////////////////
public:
// Cartesian
- void AddDetector(G4ThreeVector POS, string Shape);
+ void AddDetector(G4ThreeVector POS);
// Spherical
- void AddDetector(double R,double Theta,double Phi,string Shape);
+ void AddDetector(double R,double Theta,double Phi);
//Cylindrical
- void AddDetector2(double R,double Alpha,double Zeta,string Shape);
+ void AddDetector(double R,double Alpha,double Zeta, int Ring);
G4LogicalVolume* BuildSquareDetector();
G4LogicalVolume* BuildCylindricalDetector();
private:
+
G4LogicalVolume* m_SquareDetector;
- G4LogicalVolume* m_SquareDetector_Can;
- G4LogicalVolume* m_Square2Detector_Can;
- G4LogicalVolume* m_CylindricalDetector;
- G4LogicalVolume* m_SquareDetector_CanMgO;
- G4LogicalVolume* m_SquareDetector_Crystal;
+ G4LogicalVolume* Log_Dali_box;
+ G4LogicalVolume* Log_Al_Cryst_can;
+ G4LogicalVolume* Log_MgO_Cryst_can;
+ G4LogicalVolume* Log_Crystal;
G4LogicalVolume* lAlPMT;
+ G4LogicalVolume* Log_Dali_1Volume;
+ G4LogicalVolume* Log_Dali_3Volume;
+ G4LogicalVolume* m_CylindricalDetector;
G4LogicalVolume* lMuPMT;
G4LogicalVolume* lTopPlatePMT;
G4LogicalVolume* lGlassPMT;
G4LogicalVolume* AriaExtrude;
- G4LogicalVolume* Logic_ArrayDali_1;
- G4Material* MgO;
G4Material* NaI_Tl;
-
-
-
////////////////////////////////////////////////////
////// Inherite from NPS::VDetector class /////////
////////////////////////////////////////////////////
@@ -91,7 +89,6 @@ class Dali : public NPS::VDetector{
// Definition of materials
// Called in ConstructDetector()
void DefinitionMaterials();
-
// Construct detector and inialise sensitive part.
// Called After DetecorConstruction::AddDetector Method
@@ -125,13 +122,10 @@ class Dali : public NPS::VDetector{
vector<double> m_Zeta;
vector<double> m_R;
vector<double> m_Alpha;
+ vector<int> m_Ring;
- // Shape type
- vector<string> m_Shape ;
-
// Visualisation Attribute
G4VisAttributes* m_VisSquare;
- G4VisAttributes* m_VisCylinder;
// Needed for dynamic loading of the library
public:
diff --git a/NPSimulation/Detectors/Minos/Minos.cc b/NPSimulation/Detectors/Minos/Minos.cc
index 121b390ae..a9cc583fb 100644
--- a/NPSimulation/Detectors/Minos/Minos.cc
+++ b/NPSimulation/Detectors/Minos/Minos.cc
@@ -80,18 +80,19 @@ using namespace CLHEP;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
namespace Minos_NS{
- // TPC
- const G4double ChamberInnerRadius = 37.*mm;
- const G4double ChamberThickness = 2.*mm;
- const G4double ChamberLength = 300.*mm;
- const G4double InnerRohacellThickness = 1.*mm;
- const G4double KaptonThickness = 0.125*mm;
- const G4double OuterRohacellThickness = 2.*mm;
- const G4double TPCRadiusExt = 91.525*mm;
-
- // MINOS
- const G4double TargetRadius = 28.*mm;
- const G4double WindowThickness = 0.150*mm;
+// TPC
+const G4double ChamberInnerRadius = 37.*mm;
+/* const G4double ChamberInnerRadius = 29*mm; */ //big TPC
+const G4double ChamberThickness = 2.*mm;
+const G4double ChamberLength = 300.*mm;
+const G4double KaptonThickness = 0.125*mm;
+const G4double RohacellThickness = 2.*mm;
+const G4double TPCRadiusExt = 91.525*mm;
+/* const G4double TPCRadiusExt = 150*mm; //big TPC */
+
+// MINOS
+const G4double TargetRadius = 28.*mm;
+const G4double WindowThickness = 0.150*mm;
}
@@ -101,49 +102,41 @@ using namespace Minos_NS;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Minos Specific Method
Minos::Minos(){
- m_Event = new TMinosData() ;
- m_MinosPadScorer = 0;
- m_ReactionRegion=NULL;
-
- // RGB Color + Transparency
- m_VisTarget= new G4VisAttributes(G4Colour(0.6,1.,1., .4));
- m_VissimpleBox= new G4VisAttributes(G4Colour(0,1,0,.6));
- m_VisTPC= new G4VisAttributes(G4Colour(1.,0.5,0.6,0.3));
- m_VisInnerRohacell= new G4VisAttributes(G4Colour(1.,1.,1., .3));
- m_VisOuterRohacell= new G4VisAttributes(G4Colour(1.,1.,1., .4));
- m_VisOuterOuterRohacell= new G4VisAttributes(G4Colour(1.,1.,1., .8));
- m_VisKapton = new G4VisAttributes(G4Colour(1.,1.,0.6,0.4));
- m_VisTargetCell = new G4VisAttributes(G4Colour(0,0,1, .4));
- m_VisTargetCell->SetForceSolid(true);
- m_VisOuterKapton = new G4VisAttributes(G4Colour(1.,1.,0.6,0.8));
-
- Raw_Signal = new TH1F("raw_Signal","raw_Signal",350,0,350); // 30ns per bin
- Elec_Signal = new TH1F("Elec_Signal","Elec_Signal",350,0,350);
-
+m_Event = new TMinosData() ;
+m_MinosPadScorer = 0;
+m_ReactionRegion=NULL;
+
+// RGB Color + Transparency
+m_VisTarget= new G4VisAttributes(G4Colour(0.6,1.,1., .4));
+m_VissimpleBox= new G4VisAttributes(G4Colour(0,1,0,.6));
+m_VisTPC= new G4VisAttributes(G4Colour(1.,0.5,0.6,0.3));
+m_VisRohacell= new G4VisAttributes(G4Colour(1.,1.,1., .8));
+m_VisKapton = new G4VisAttributes(G4Colour(1.,1.,0.6,0.4));
+m_VisTargetCell = new G4VisAttributes(G4Colour(0,0,1, .4));
+m_VisTargetCell->SetForceSolid(true);
+m_VisOuterKapton = new G4VisAttributes(G4Colour(1.,1.,0.6,0.8));
+
+Raw_Signal = new TH1F("raw_Signal","raw_Signal",350,0,350); // 30ns per bin
+Elec_Signal = new TH1F("Elec_Signal","Elec_Signal",350,0,350);
+
/* Raw_Signal = new TH1F;*/
/* Elec_Signal = new TH1F;*/
-
- fa1 = new TF1("fa1","[0]*exp(-3.*(x-[1])/[2]) * sin((x-[1])/[2]) * pow((x-[1])/[2], 3)",0,1000);
-
- solidTarget=0;
- logicTarget=0;
- solidChamber=0;
- logicChamber=0;
- solidTPC=0;
- logicTPC=0;
- solidWindow0=0;
- logicWindow0=0;
- solidWindow1=0;
- logicWindow1=0;
- solidWindow2=0;
- logicWindow2=0;
- solidInnerRohacell=0;
- logicInnerRohacell=0;
- solidOuterRohacell=0;
- logicOuterRohacell=0;
- solidKapton=0;
- logicKapton=0;
+
+fa1 = new TF1("fa1","[0]*exp(-3.*(x-[1])/[2]) * sin((x-[1])/[2]) * pow((x-[1])/[2], 3)",0,1000);
+
+solidTarget=0;
+logicTarget=0;
+solidChamber=0;
+logicChamber=0;
+solidTPC=0;
+logicTPC=0;
+solidWindow0=0;
+logicWindow0=0;
+solidRohacell=0;
+logicRohacell=0;
+solidKapton=0;
+logicKapton=0;
}
Minos::~Minos(){
@@ -154,208 +147,19 @@ Minos::~Minos(){
/* void Minos::AddDetector(G4ThreeVector POS, double LengthOfTarget, int PresenceOfMinos){ */
void Minos::AddDetector(G4ThreeVector POS, double LengthOfTarget, G4String MaterialOfTarget,G4String MaterialOfCell, int PresenceOfMinos){
- m_POS.push_back(POS);
- m_TargetLength.push_back(LengthOfTarget);
- m_TargetMaterial.push_back(MaterialOfTarget);
- m_CellMaterial.push_back(MaterialOfCell);
- m_TPCOnly.push_back(PresenceOfMinos);
+m_POS.push_back(POS);
+m_TargetLength.push_back(LengthOfTarget);
+m_TargetMaterial.push_back(MaterialOfTarget);
+m_CellMaterial.push_back(MaterialOfCell);
+m_TPCOnly.push_back(PresenceOfMinos);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-void Minos::DefineMaterials()
-{
- //This function illustrates the possible ways to define materials
-
- G4String symbol; //a=mass of a mole;
- G4double a, z, density; //z=mean number of protons;
-
- G4int ncomponents, natoms;
- G4double fractionmass;
-
- //
- // define Elements
- //
-
- G4Element* H = new G4Element("Hydrogen",symbol="H" , z= 1., a= 1.01*g/mole);
- G4Element* C = new G4Element("Carbon" ,symbol="C" , z= 6., a= 12.01*g/mole);
- G4Element* F = new G4Element("Fluorin" ,symbol="F" , z= 9., a= 19.0*g/mole);
- G4Element* N = new G4Element("Nitrogen",symbol="N" , z= 7., a= 14.01*g/mole);
- G4Element* O = new G4Element("Oxygen" ,symbol="O" , z= 8., a= 16.00*g/mole);
- G4Element* Ar = new G4Element("Argon", symbol="Ar", z=18, a=39.948*g/mole);
- G4Element* Fe = new G4Element("Fer",symbol="Fe" , z= 26., a= 55.9*g/mole);
- G4Element* Cr = new G4Element("Chrome",symbol="Cr" , z= 24., a= 52.*g/mole);
-
- new G4Material("Aluminium", z=13., a=26.98*g/mole, density=2.700*g/cm3);
- new G4Material("liquidArgon", z=18., a= 39.95*g/mole, density= 1.390*g/cm3);
- new G4Material("Lead" , z=82., a= 207.19*g/mole, density= 11.35*g/cm3);
- new G4Material("Silicium", z=14., a=28.09*g/mole, density=2.330*g/cm3);
- new G4Material("Titanium", z=22., a=47.87*g/mole, density=4.510*g/cm3);
-
- G4Material* iso = new G4Material("isobutane", density=0.002506*g/cm3, ncomponents=2);
- iso->AddElement(C, natoms=4);
- iso->AddElement(H, natoms=10);
- G4Material* CF4 =
- new G4Material("CF4", density= 0.0036586*g/cm3, ncomponents=2);
- CF4->AddElement(C, natoms=1);
- CF4->AddElement(F, natoms=4);
- // overwrite computed meanExcitationEnergy with ICRU recommended value
- CF4->GetIonisation()->SetMeanExcitationEnergy(20.0*eV);
- G4Material* mix =
- new G4Material("mix", density= 0.0019836*g/cm3, ncomponents=3);
- mix->AddMaterial(CF4, fractionmass=15.*perCent);
- mix->AddMaterial(iso, fractionmass=3.*perCent);
- mix->AddElement(Ar, fractionmass=82.*perCent);
- // overwrite computed meanExcitationEnergy with ICRU recommended value
- mix->GetIonisation()->SetMeanExcitationEnergy(25.0*eV);
-
- ///////// Drift properties
-
- G4MaterialPropertiesTable* MPT = new G4MaterialPropertiesTable();
- MPT->AddConstProperty("DE_PAIRENERGY",30*eV);
- MPT->AddConstProperty("DE_ABSLENGTH",10*pc);
- MPT->AddConstProperty("DE_DRIFTSPEED",3.6*cm/microsecond);
- MPT->AddConstProperty("DE_TRANSVERSALSPREAD",7e-5*mm2/ns);
- MPT->AddConstProperty("DE_LONGITUDINALSPREAD",7e-5*mm2/ns);
- /* MPT->AddConstProperty("DE_TRANSVERSALSPREAD",7e-15*mm2/ns); */
- /* MPT->AddConstProperty("DE_LONGITUDINALSPREAD",7e-15*mm2/ns); */
- mix->SetMaterialPropertiesTable(MPT);
-
- G4Material* Ar_CF4_95_5 =
- new G4Material("Ar_CF4_95_5", density= 0.0017611*g/cm3, ncomponents=2);
- Ar_CF4_95_5->AddMaterial(CF4, fractionmass=5.*perCent);
- Ar_CF4_95_5->AddElement(Ar, fractionmass=95.*perCent);
- Ar_CF4_95_5->GetIonisation()->SetMeanExcitationEnergy(30.0*eV);
-
- G4Material* Ar_CF4_90_10 =
- new G4Material("Ar_CF4_90_10", density= 0.0018610*g/cm3, ncomponents=2);
- Ar_CF4_90_10->AddMaterial(CF4, fractionmass=10.*perCent);
- Ar_CF4_90_10->AddElement(Ar, fractionmass=90.*perCent);
- Ar_CF4_90_10->GetIonisation()->SetMeanExcitationEnergy(25.0*eV);
-
- G4Material* Ar_iso_97_3 =
- new G4Material("Ar_iso_97_3", density= 0.0016838*g/cm3, ncomponents=2);
- Ar_iso_97_3->AddMaterial(iso, fractionmass=3.*perCent);
- Ar_iso_97_3->AddElement(Ar, fractionmass=97.*perCent);
- Ar_iso_97_3->GetIonisation()->SetMeanExcitationEnergy(25.0*eV);
-
- G4Material* Ar_iso_95_5 =
- new G4Material("Ar_iso_95_5", density= 0.0016990*g/cm3, ncomponents=2);
- Ar_iso_95_5->AddMaterial(iso, fractionmass=5.*perCent);
- Ar_iso_95_5->AddElement(Ar, fractionmass=95.*perCent);
- Ar_iso_95_5->GetIonisation()->SetMeanExcitationEnergy(25.0*eV);
-
- G4Material* LH2 =
- new G4Material("LH2", density= 0.07293*g/cm3, ncomponents=1);
- LH2->AddElement(H, natoms=2);
- G4Material* Myl =
- new G4Material("Mylar", density= 1.397*g/cm3, ncomponents=3);
- Myl->AddElement(C, natoms=10);
- Myl->AddElement(H, natoms= 8);
- Myl->AddElement(O, natoms= 4);
-
- G4Material* Epo =
- new G4Material("Epoxy", density= 1.85*g/cm3, ncomponents=3); //density of FR4 (Wikipedia)
- Epo->AddElement(C, natoms=18);
- Epo->AddElement(H, natoms= 20);
- Epo->AddElement(O, natoms= 3);
-
- G4Material* Air =
- new G4Material("Air" , density= 1.290*mg/cm3, ncomponents=3);
- Air->AddElement(N, fractionmass=0.781);
- Air->AddElement(O, fractionmass=0.21);
- Air->AddElement(Ar, fractionmass=0.009);
-
- /* G4Material* Vacuum = */
- /* new G4Material("Vacuum", z=1., a=1.01*g/mole,density= universe_mean_density, */
- /* kStateGas, 2.73*kelvin, 3.e-18*pascal); */
-
- G4Material* Vacuum =
- new G4Material("Vacuum", density= universe_mean_density, ncomponents=2);
- Vacuum-> AddElement(N,.7);
- Vacuum-> AddElement(O,.3);
-
- G4Material* beam =
- new G4Material("Beam", density= 1.e-5*g/cm3, ncomponents=1,
- kStateGas, STP_Temperature, 2.e-2*bar);
- beam->AddMaterial(Air, fractionmass=1.);
-
- G4Material* Inox =
- new G4Material("Inox", density= 8.02*g/cm3, ncomponents=3);
- Inox->AddElement(C, fractionmass=0.001);
- Inox->AddElement(Fe, fractionmass=0.829);
- Inox->AddElement(Cr, fractionmass=0.17);
-
- G4Material* Kapton =
- new G4Material("Kapton", density= 1.42*g/cm3, ncomponents=4);
- Kapton->AddElement(C, fractionmass=0.691133);
- Kapton->AddElement(H, fractionmass=0.026362);
- Kapton->AddElement(O, fractionmass=0.209235);
- Kapton->AddElement(N, fractionmass=0.073270);
-
- G4Material* Rohacell =
- new G4Material("Rohacell", density= 0.075*g/cm3, ncomponents=4);
- Rohacell->AddElement(C, fractionmass=0.6014);
- Rohacell->AddElement(H, fractionmass=0.0805);
- Rohacell->AddElement(O, fractionmass=0.3154);
- Rohacell->AddElement(N, fractionmass=0.00276);
-
- //default materials of the World
- defaultMaterial = Vacuum;
-}
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-//SET MATERIALS
-void Minos::SetTargetMaterial(G4String materialChoice)
-{
- // search the material by its name
- G4Material* pttoMaterial = G4Material::GetMaterial(materialChoice);
- if (pttoMaterial) TargetMaterial = pttoMaterial;
-}
-void Minos::SetChamberMaterial(G4String materialChoice)
-{
- // search the material by its name
- G4Material* pttoMaterial = G4Material::GetMaterial(materialChoice);
- if (pttoMaterial) ChamberMaterial = pttoMaterial;
-}
-
-void Minos::SetTPCMaterial(G4String materialChoice)
-{
- // search the material by its name
- G4Material* pttoMaterial = G4Material::GetMaterial(materialChoice);
- if (pttoMaterial) TPCMaterial = pttoMaterial;
-}
-
-void Minos::SetWindowMaterial(G4String materialChoice)
-{
- // search the material by its name
- G4Material* pttoMaterial = G4Material::GetMaterial(materialChoice);
- if (pttoMaterial) WindowMaterial = pttoMaterial;
-}
-void Minos::SetInnerRohacellMaterial(G4String materialChoice)
-{
- // search the material by its name
- G4Material* pttoMaterial = G4Material::GetMaterial(materialChoice);
- if (pttoMaterial) InnerRohacellMaterial = pttoMaterial;
-}
-void Minos::SetOuterRohacellMaterial(G4String materialChoice)
-{
- // search the material by its name
- G4Material* pttoMaterial = G4Material::GetMaterial(materialChoice);
- if (pttoMaterial) OuterRohacellMaterial = pttoMaterial;
-}
-void Minos::SetKaptonMaterial(G4String materialChoice)
-{
- // search the material by its name
- G4Material* pttoMaterial = G4Material::GetMaterial(materialChoice);
- if (pttoMaterial) KaptonMaterial = pttoMaterial;
-}
-
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-// Below are vis attributes that permits someone to test / play
// with the interactive expansion / contraction geometry system of the
// vis/OpenInventor driver :
@@ -376,122 +180,90 @@ G4LogicalVolume* Minos::BuildTarget(){
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-G4LogicalVolume* Minos::BuildChamber(){
- if(!logicChamber){
- //
- // Chamber
+G4LogicalVolume* Minos::BuildTPC(){
+ if(!logicTPC){
+ //
+ // TPC
//
- solidChamber = new G4Tubs("Chamber", //its name
- ChamberInnerRadius,ChamberInnerRadius+ChamberThickness,ChamberLength/2.,0,360.); //size
- logicChamber = new G4LogicalVolume(solidChamber, //its solid
- ChamberMaterial, //its material
- "Chamber"); //its name
- m_VissimpleBox->SetVisibility(true);
- logicChamber->SetVisAttributes(m_VissimpleBox);
+ solidTPC = new G4Tubs("TPC",
+ ChamberInnerRadius ,TPCRadiusExt,ChamberLength/2.,0,360.);
+ logicTPC = new G4LogicalVolume(solidTPC, //its solid
+ TPCMaterial, //its material
+ "TPC"); //name
+ logicTPC->SetVisAttributes(m_VisTPC);
}
- return logicChamber;
+ return logicTPC;
}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+/* G4LogicalVolume* Minos::BuildChamber(){ */
+/* if(!logicChamber){ */
+/* solidChamber = new G4Tubs("Chamber", //its name */
+/* ChamberInnerRadius,ChamberInnerRadius+ChamberThickness,ChamberLength/2.,0,360.); //size */
+/* logicChamber = new G4LogicalVolume(solidChamber, //its solid */
+/* ChamberMaterial, //its material */
+/* "Chamber"); //its name */
+/* m_VissimpleBox->SetVisibility(true); */
+/* logicChamber->SetVisAttributes(m_VissimpleBox); */
+/* } */
+/* return logicChamber; */
+/* } */
+
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4LogicalVolume* Minos::BuildInnerRohacell(){
- if(!logicInnerRohacell){
- //
- // Inner Rohacell
- //
- solidInnerRohacell = new G4Tubs("InnerRohacell", //its name
- ChamberInnerRadius /*+ ChamberThickness*/,ChamberInnerRadius + ChamberThickness+InnerRohacellThickness,ChamberLength/2.,0,360.); //size
- logicInnerRohacell = new G4LogicalVolume(solidInnerRohacell, //its solid
- InnerRohacellMaterial, //its material
+ if(!logicRohacell){
+ solidRohacell = new G4Tubs("InnerRohacell", //its name
+ ChamberInnerRadius ,ChamberInnerRadius + RohacellThickness ,ChamberLength/2.,0,360.); //size
+ logicRohacell = new G4LogicalVolume(solidRohacell, //its solid
+ RohacellMaterial, //its material
"InnerRohacell"); //its name
- logicInnerRohacell->SetVisAttributes(m_VisInnerRohacell);
+ logicRohacell->SetVisAttributes(m_VisRohacell);
}
- return logicInnerRohacell;
+ return logicRohacell;
}
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-G4LogicalVolume* Minos::BuildOuterRohacell(){
- if(!logicOuterRohacell){
- //
- // Outer Rohacell
- //
- solidOuterRohacell = new G4Tubs("OuterRohacell", //its name
- ChamberInnerRadius /*+ ChamberThickness + InnerRohacellThickness + KaptonThickness*/,ChamberInnerRadius + ChamberThickness + InnerRohacellThickness + KaptonThickness+OuterRohacellThickness,ChamberLength/2.,0,360.); //size
- logicOuterRohacell = new G4LogicalVolume(solidOuterRohacell, //its solid
- OuterRohacellMaterial, //its material
- "OuterRohacell"); //its name
- logicOuterRohacell->SetVisAttributes(m_VisOuterRohacell);
- }
- return logicOuterRohacell;
-}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-//I redefine the Rohacell for bigger cylinder
-G4LogicalVolume* Minos::BuildOuterOuterRohacell(){
- if(logicOuterRohacell){
- //
- // Outer Rohacell
- //
- solidOuterRohacell = new G4Tubs("OuterRohacell", //its name
- TPCRadiusExt-OuterRohacellThickness-KaptonThickness, TPCRadiusExt ,ChamberLength/2.,0,360.); //size
- logicOuterRohacell = new G4LogicalVolume(solidOuterRohacell, //its solid
- OuterRohacellMaterial, //its material
+
+G4LogicalVolume* Minos::BuildOuterRohacell(){
+ if(logicRohacell){
+ solidRohacell = new G4Tubs("OuterRohacell", //its name
+ TPCRadiusExt-RohacellThickness, TPCRadiusExt ,ChamberLength/2.,0,360.); //size
+ logicRohacell = new G4LogicalVolume(solidRohacell, //its solid
+ RohacellMaterial, //its material
"OuterRohacell"); //its name
- logicOuterRohacell->SetVisAttributes(m_VisOuterOuterRohacell);
+ logicRohacell->SetVisAttributes(m_VisRohacell);
}
- return logicOuterRohacell;
+ return logicRohacell;
}
-
-
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4LogicalVolume* Minos::BuildKapton(){
if(!logicKapton){
- //
- // Kapton
- //
solidKapton = new G4Tubs("Kapton", //its name
- ChamberInnerRadius+ ChamberThickness +InnerRohacellThickness ,ChamberInnerRadius + ChamberThickness+InnerRohacellThickness+KaptonThickness,ChamberLength/2.,0,360.); //size
+ ChamberInnerRadius+RohacellThickness ,ChamberInnerRadius +RohacellThickness+KaptonThickness,ChamberLength/2.,0,360.); //size
logicKapton = new G4LogicalVolume(solidKapton, //its solid
KaptonMaterial, //its material
"Kapton"); //its name
logicKapton->SetVisAttributes(m_VisKapton);
}
-
return logicKapton;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4LogicalVolume* Minos::BuildOuterKapton(){
if(logicKapton){
- //
- // Kapton
- //
solidKapton = new G4Tubs("Kapton", //its name
- TPCRadiusExt-OuterRohacellThickness-KaptonThickness, TPCRadiusExt-OuterRohacellThickness,ChamberLength/2.,0,360.); //size
-
+ TPCRadiusExt-RohacellThickness-KaptonThickness, TPCRadiusExt-RohacellThickness,ChamberLength/2.,0,360.); //size
logicKapton = new G4LogicalVolume(solidKapton, //its solid
KaptonMaterial, //its material
"Kapton"); //its name
- logicKapton->SetVisAttributes(m_VisOuterKapton);
+ logicKapton->SetVisAttributes(m_VisOuterKapton);
}
-
-
return logicKapton;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-G4LogicalVolume* Minos::BuildTPC(){
- if(!logicTPC){
- //
- // TPC
- //
- solidTPC = new G4Tubs("TPC",
- ChamberInnerRadius /*+ ChamberThickness + InnerRohacellThickness + KaptonThickness + OuterRohacellThickness*/,TPCRadiusExt,ChamberLength/2.,0,360.);
- logicTPC = new G4LogicalVolume(solidTPC, //its solid
- TPCMaterial, //its material
- "TPC"); //name
- logicTPC->SetVisAttributes(m_VisTPC);
- }
- return logicTPC;
-}
+
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -503,7 +275,6 @@ G4LogicalVolume* Minos::BuildWindow0(){
if(!logicWindow0){
solidWindow0 = new G4Tubs("WindowTube", //its name
0.,TargetRadius+WindowThickness,TargetLength/2.+WindowThickness,0,360.);
-
logicWindow0 = new G4LogicalVolume(solidWindow0, //its solid
WindowMaterial, //its material
"WindowTube"); //name
@@ -511,36 +282,6 @@ G4LogicalVolume* Minos::BuildWindow0(){
}
return logicWindow0;
}
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-
-G4LogicalVolume* Minos::BuildWindow1(){
- if(!logicWindow1){
- solidWindow1 = new G4Tubs("WindowEntrance", //its name
- 0.,TargetRadius+WindowThickness,WindowThickness,0,360.);
-
- logicWindow1 = new G4LogicalVolume(solidWindow1, //its solid
- WindowMaterial, //its material
- "WindowEntrance"); //name
- logicWindow1->SetVisAttributes(m_VisTargetCell);
- }
- return logicWindow1;
-}
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-G4LogicalVolume* Minos::BuildWindow2(){
- if(!logicWindow2){
- solidWindow2 = new G4Tubs("WindowOutcoming", //its name
- 0.,TargetRadius+WindowThickness,WindowThickness,0,360.);
-
- logicWindow2 = new G4LogicalVolume(solidWindow2, //its solid
- WindowMaterial, //its material
- "WindowOutcoming"); //name
- logicWindow2->SetVisAttributes(m_VisTargetCell);
- }
- return logicWindow2;
-}
-
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -559,15 +300,15 @@ void Minos::ReadConfiguration(NPL::InputParser parser){
if(blocks[i]->HasTokenList(cart)){
if(NPOptionManager::getInstance()->GetVerboseLevel())
cout << endl << "//// Minos " << i+1 << endl;
-
- G4ThreeVector Pos = NPS::ConvertVector(blocks[i]->GetTVector3("POS","mm"));
- TargetLength = blocks[i]->GetDouble("TargetLength","mm");
- G4String TargetMaterialname = blocks[i]->GetString("TargetMaterial");
- G4String CellMaterial = blocks[i]->GetString("CellMaterial");
- TPCOnly = blocks[i]->GetInt("TPCOnly");
- AddDetector(Pos,TargetLength,TargetMaterialname, CellMaterial, TPCOnly);
- /* AddDetector(Pos,TargetLength, TPCOnly); */
-
+
+ G4ThreeVector Pos = NPS::ConvertVector(blocks[i]->GetTVector3("POS","mm"));
+ TargetLength = blocks[i]->GetDouble("TargetLength","mm");
+ G4String TargetMaterialname = blocks[i]->GetString("TargetMaterial");
+ G4String CellMaterial = blocks[i]->GetString("CellMaterial");
+ TPCOnly = blocks[i]->GetInt("TPCOnly");
+ AddDetector(Pos,TargetLength,TargetMaterialname, CellMaterial, TPCOnly);
+ /* AddDetector(Pos,TargetLength, TPCOnly); */
+
}
else{
cout << "ERROR: check your input file formatting " << endl;
@@ -585,22 +326,40 @@ void Minos::ConstructDetector(G4LogicalVolume* world){
TargetLength = m_TargetLength[i];
TPCOnly = m_TPCOnly[i];
- DefineMaterials();
-
- SetTargetMaterial(m_TargetMaterial[i]);
- SetChamberMaterial("Inox");
- SetTPCMaterial("mix");
- SetWindowMaterial(m_CellMaterial[i]);
- SetKaptonMaterial("Kapton");
- SetInnerRohacellMaterial("Rohacell");
- SetOuterRohacellMaterial("Rohacell");
+ /* TargetMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary("Vacuum"); */
+ TargetMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary("LH2");
+ WindowMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary("Mylar");
+ ChamberMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary("Al");
+ KaptonMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary("Kapton");
+ RohacellMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary("Rohacell");
+ TPCMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary("mixMINOS");
+
+ /* TargetMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary("Vacuum"); */
+ /* WindowMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary("Vacuum"); */
+ /* ChamberMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary("Vacuum"); */
+ /* KaptonMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary("Vacuum"); */
+ /* RohacellMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary("Vacuum"); */
+ /* TPCMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary("mixMINOS"); */
+
+ ///////// Drift properties
+ G4MaterialPropertiesTable* MPT = new G4MaterialPropertiesTable();
+ MPT->AddConstProperty("DE_PAIRENERGY",30*eV);
+ MPT->AddConstProperty("DE_ABSLENGTH",10*pc);
+ MPT->AddConstProperty("DE_DRIFTSPEED",3.475*cm/microsecond);
+ MPT->AddConstProperty("DE_TRANSVERSALSPREAD",14e-5*mm2/ns);
+ MPT->AddConstProperty("DE_LONGITUDINALSPREAD",14e-5*mm2/ns);
+
+ /* MPT->AddConstProperty("DE_TRANSVERSALSPREAD",0*mm2/ns); */
+ /* MPT->AddConstProperty("DE_LONGITUDINALSPREAD",0*mm2/ns); */
+
+ TPCMaterial->SetMaterialPropertiesTable(MPT);
G4ThreeVector Det_pos = m_POS[i] ;
-
+
double MinosX = Det_pos.x();
double MinosY= Det_pos.y();
double MinosZ = Det_pos.z() + m_TargetLength[i]/2. ;
-
+
new G4PVPlacement(0,//its name
G4ThreeVector(0,0,+ChamberLength/2 ),
/* G4ThreeVector(wX,wY, wZ + ChamberLength - m_TargetLength[i]-WindowThickness*2. - 5*mm ), // Z positioning putting TPC beginn and Target beginning w/ difference of 5mm */
@@ -610,7 +369,7 @@ void Minos::ConstructDetector(G4LogicalVolume* world){
false, //no boolean operation
0); //copy number
-//////// ELECTRIC FIELD
+ //////// ELECTRIC FIELD
G4ElectricField* field = new G4UniformElectricField(G4ThreeVector(0.0,0.0,200*volt/cm));
// Create an equation of motion for this field
@@ -630,95 +389,87 @@ void Minos::ConstructDetector(G4LogicalVolume* world){
FieldManager->SetChordFinder( ChordFinder );
G4Material* Cu
- = MaterialManager::getInstance()->GetMaterialFromLibrary("Cu");
-
- ///////// CONSTRUCT PADS using parametrized volumes or replica (uncomment the correct section)
-
+ = MaterialManager::getInstance()->GetMaterialFromLibrary("Cu");
+
+ ///////// CONSTRUCT PADS using parametrized volumes or replica (uncomment the correct section)
+
//// Using REPLICA
// Construct Pads ring by ring
for (int RingNbr = 0; RingNbr < 18; RingNbr++){
- int PadsPerRing[18]={144,152,156,164,172,176,184,192,196,204,212,216,224,228,236,244,248,256};
- /* G4double InnerRadius = (45.2+RingNbr*2.1+0.02)*mm;// 0.02 mm between each ring */
- G4double InnerRadius = (45.75+RingNbr*2.1)*mm;
- G4double OuterRadius = (47.75+RingNbr*2.1)*mm;
-
- //Volume where are placed replica pads
- G4VSolid* AnodeRing = new G4Tubs("ring",InnerRadius,OuterRadius,
- 0.01*mm,0.,360*deg);
- G4LogicalVolume * AnodeRing_log = new G4LogicalVolume(AnodeRing,
- TPCMaterial, "ringL", 0, 0, 0);
-
- {G4VisAttributes* atb= new G4VisAttributes(G4Colour(0.8, 0.4, 0.,0.4));
+ int PadsPerRing[18]={144,152,156,164,172,176,184,192,196,204,212,216,224,228,236,244,248,256};
+ /* G4double InnerRadius = (45.2+RingNbr*2.1+0.02)*mm;// 0.02 mm between each ring */
+ G4double InnerRadius = (45.75+RingNbr*2.1)*mm;
+ G4double OuterRadius = (47.75+RingNbr*2.1)*mm;
+
+ //Volume where are placed replica pads
+ G4VSolid* AnodeRing = new G4Tubs("ring",InnerRadius,OuterRadius,
+ 0.01*mm,0.,360*deg);
+ G4LogicalVolume * AnodeRing_log = new G4LogicalVolume(AnodeRing,
+ /* mix, "ringL", 0, 0, 0); */
+ TPCMaterial, "ringL", 0, 0, 0);
+
+ {G4VisAttributes* atb= new G4VisAttributes(G4Colour(0.8, 0.4, 0.,0.4));
AnodeRing_log->SetVisAttributes(atb);}
- new G4PVPlacement(0,G4ThreeVector(0,0,-ChamberLength/2+0.01*mm),
- AnodeRing_log,"ring", logicTPC, false, RingNbr);
+ new G4PVPlacement(0,G4ThreeVector(0,0,-ChamberLength/2+0.01*mm),
+ AnodeRing_log,"ring", logicTPC, false, RingNbr);
+
+ G4double Pad_dPhi = (360./(PadsPerRing[RingNbr]+1))*deg; // longitudinal component of Pad
+ G4double Pad_shift = (360./PadsPerRing[RingNbr])*deg; // dPhi between each Pads
- G4double Pad_dPhi = (360./(PadsPerRing[RingNbr]+1))*deg; // longitudinal component of Pad
- G4double Pad_shift = (360./PadsPerRing[RingNbr])*deg; // dPhi between each Pads
+ G4VSolid* Pad = new G4Tubs("div_ring", InnerRadius, OuterRadius,
+ 0.01*mm,0, Pad_dPhi);
- G4VSolid* Pad = new G4Tubs("div_ring", InnerRadius, OuterRadius,
- 0.01*mm,0, Pad_dPhi);
-
- G4LogicalVolume* Pad_log = new G4LogicalVolume(Pad,
- Cu,"div_ringL",0,0,0);
-
- {G4VisAttributes* atb= new G4VisAttributes(G4Colour(0.8, 0.4, 0.,0.4));
+ G4LogicalVolume* Pad_log = new G4LogicalVolume(Pad,
+ Cu,"div_ringL",0,0,0);
+
+ {G4VisAttributes* atb= new G4VisAttributes(G4Colour(0.8, 0.4, 0.,0.4));
Pad_log->SetVisAttributes(atb);}
- Pad_log->SetSensitiveDetector(m_MinosPadScorer);
-
- new G4PVReplica("div_ring_phys", Pad_log,
- AnodeRing_log, kPhi, PadsPerRing[RingNbr],Pad_shift ,0);
-
+ Pad_log->SetSensitiveDetector(m_MinosPadScorer);
+
+ new G4PVReplica("div_ring_phys", Pad_log,
+ AnodeRing_log, kPhi, PadsPerRing[RingNbr],Pad_shift ,0);
}
-
-//////////////////////////////////////
-
+
+ //////////////////////////////////////
+
new G4PVPlacement(0, //its name
G4ThreeVector(0,0,0), //at (0,0,0)
- BuildOuterRohacell(), //its logical volume
- "Rohacell", //its name
+ BuildInnerRohacell(), //its logical volume
+ "InnerRohacell", //its name
logicTPC, //its mother volume
false, //no boolean operation
0); //copy number
- //check the order and positioning of kapton , chamber and rohacell from pag 16 of Thesis Clementine
-
- new G4PVPlacement(0, G4ThreeVector(0,0,0), //at (0,0,0)
- BuildChamber(), //its logical volume
- "Chamber", //its name
- logicOuterRohacell, //its mother volume
+ new G4PVPlacement(0, //its name
+ G4ThreeVector(0,0,0), //at (0,0,0)
+ BuildOuterRohacell(), //its logical volume
+ "OuterRohacell", //its name
+ logicTPC, //its mother volume
+ /* logicTPC, //its mother volume */
false, //no boolean operation
0); //copy number
new G4PVPlacement(0, G4ThreeVector(0,0,0), //at (0,0,0)
BuildKapton(), //its logical volume
- "Kapton", //its name
- logicOuterRohacell, //its mother volume
- false, //no boolean operation
- 0); //copy number
-
- new G4PVPlacement(0, //its name
- G4ThreeVector(0,0,0), //at (0,0,0)
- BuildOuterOuterRohacell(), //its logical volume
- "Rohacell"/*"OuterRohacell"*/, //its name
- logicTPC/*world*/, //its mother volume
+ "InnerKapton", //its name
+ logicTPC, //its mother volume
false, //no boolean operation
0); //copy number
new G4PVPlacement(0, //its name
G4ThreeVector(0,0,0), //at (0,0,0)
BuildOuterKapton(), //its logical volume
- "Kapton", //its name
- logicOuterRohacell, //its mother volume
+ "OuterKapton", //its name
+ logicTPC, //its mother volume
false, //no boolean operation
0); //copy number
if(TPCOnly!=1){
new G4PVPlacement(0, //its name
- G4ThreeVector(MinosX,MinosY, MinosZ), //at (0,0,0)
+ G4ThreeVector(MinosX,MinosY,MinosZ), //at (0,0,0)
BuildWindow0(), //its logical volume
"WindowTube", //its name
world, //its mother volume
@@ -738,7 +489,7 @@ void Minos::ConstructDetector(G4LogicalVolume* world){
m_ReactionRegion -> AddRootLogicalVolume(logicTarget);
m_ReactionRegion->SetUserLimits(new G4UserLimits(1.2*mm)); //???
}
-
+
G4FastSimulationManager* mng = m_ReactionRegion->GetFastSimulationManager();
unsigned int size = m_ReactionModel.size();
for(unsigned int o = 0 ; o < size ; o++){
@@ -753,7 +504,7 @@ void Minos::ConstructDetector(G4LogicalVolume* world){
m_ReactionModel.push_back(fsm);
}
}
-
+
// Visualization attributes
world->SetVisAttributes (G4VisAttributes::Invisible);
@@ -777,21 +528,23 @@ void Minos::InitializeRootOutput(){
void Minos::ReadSensitive(const G4Event* ){
m_Event->Clear();
-///////////
-// DriftElectron scorer
+ ///////////
+ // DriftElectron scorer
CylinderTPCScorers::PS_TPCAnode* Scorer2= (CylinderTPCScorers::PS_TPCAnode*) m_MinosPadScorer->GetPrimitive(0);
- unsigned int size2 = Scorer2->GetMult();
+ unsigned int size2 = Scorer2->GetMult();
for(unsigned int i = 0 ; i < size2 ; i++){
-
+
int Pad = Scorer2->GetPad(i);
double X = Scorer2->GetX(i);
double Y = Scorer2->GetY(i);
-
+
m_Event->SetCharge(Pad,X,Y);
- /* m_Event->AddChargePoint(Scorer2->GetQ(i), Scorer2->GetT(i)); */
-
+ //m_Event->AddChargePoint(Scorer2->GetQ(i), Scorer2->GetT(i));
+
SimulateGainAndDigitizer(Scorer2->GetQ(i), Scorer2->GetT(i));
+
+ /* m_Event->AddChargePoint(Scorer2->GetQ(i),Scorer2->GetT(i)); */
}
}
@@ -800,81 +553,87 @@ void Minos::SimulateGainAndDigitizer(vector<double> rawQ, vector<double> rawT){
Charge2.clear();
Time.clear();
-
+
Raw_Signal->Reset();
Elec_Signal->Reset();
-
- // Reallocate electrons from each pack
-
+
// Add white noise;
- for( int bin = 0 ; bin < Elec_Signal->GetNbinsX() ; bin++)
- Elec_Signal->SetBinContent(bin,20+(10-gRandom->Uniform(20)*20));
-
- for (unsigned int j = 1; j < rawQ.size()-1; j++){
- time = (rawT[j+1]-rawT[j-1])/rawQ[0];
-
- for ( int k = -rawQ[0]/2; k < (rawQ[0]/2); k++){
- Raw_Signal->Fill((rawT[j]+k*time)/30.);
+ /* for( int bin = 0 ; bin < Elec_Signal->GetNbinsX() ; bin++) */
+ /* Elec_Signal->SetBinContent(bin,20+(10-gRandom->Uniform(20)*20)); */
+
+ // Reallocate electrons from each pack
+ /* if (rawQ.sze()==1) { */
+ /* for ( int j = 0; j < rawQ[0]; j++){ */
+ /* Raw_Signal->Fill(rawT[0]/30.); */
+ /* } */
+ /* } */
+ /* else{ */
+ /* for ( int j = 1; j < rawQ.size()-1; j++){ */
+ /* time = (rawT[j+1]-rawT[j-1])/rawQ[j]; */
+ /* for ( int k = -rawQ[0]/2; k < (rawQ[0]/2); k++){ */
+ /* Raw_Signal->Fill((rawT[j]+k*time)/30.); */
+ /* } */
+ /* } */
+ /* } */
+
+ for ( unsigned int j = 0; j < rawQ.size(); j++){
+ for ( int i = 0 ; i < rawQ[j]; i++ ){
+ Raw_Signal->Fill((rawT[j])/30.);
}
}
-
- // Application of the electronic reponse function
-
+ // Application of the electronic reponse function
for ( int x=0; x < Raw_Signal->GetNbinsX(); x++){
+
if(Raw_Signal->GetBinContent(x) < 0.5)
continue;
- else {
- start = Raw_Signal->GetBinCenter(x);
- end = Raw_Signal->GetBinCenter(x)+1200/30.;
- /* DriftTime = Raw_Signal->GetBinCenter(x); */
-
- fa1->SetRange(start, end);
- fa1->SetParameter(0,1500);
- fa1->SetParameter(1,start);
- fa1->SetParameter(2,450/30.);
- /* Elec_Signal->SetBinContent(bin,20+(10-gRandom->Uniform(20)*20)) */
- /* for (int p=0; p< 0; p++) */
- for (int p=0; p< Raw_Signal->GetBinContent(x)*1500; p++)
- Elec_Signal->Fill(fa1->GetRandom());
-
- }
- }
-
+ else{
+ start = Raw_Signal->GetBinCenter(x);
+ end = Raw_Signal->GetBinCenter(x)+1200/30.;
+ // DriftTime = Raw_Signal->GetBinCenter(x);
+ fa1->SetRange(start, end);
+ fa1->SetParameter(0,1500);
+ fa1->SetParameter(1,start);
+ fa1->SetParameter(2,450/30.);
+ // Elec_Signal->SetBinContent(bin,20+(10-gRandom->Uniform(20)*20))
+ // for (int p=0; p< 0; p++)
+ for (int p=0; p< Raw_Signal->GetBinContent(x)*1500; p++)
+ Elec_Signal->Fill(fa1->GetRandom());
+ }
+ }
for ( int bin=0; bin < Elec_Signal->GetNbinsX(); bin++){
Charge2.push_back(Elec_Signal->GetBinContent(bin));
Time.push_back(Elec_Signal->GetBinCenter(bin));
-
}
-
+
m_Event->AddChargePoint(Charge2,Time);
}
- //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
- ////////////////////////////////////////////////////////////////
- void Minos::InitializeScorers() {
- // This check is necessary in case the geometry is reloaded
- bool already_exist = false;
- bool already_exist2 = false;
- bool already_exist3 = false;
-
- m_MinosPadScorer = CheckScorer("MinosPadScorer",already_exist3) ;
-
- if(already_exist && already_exist2 && already_exist3 )
- return ;
-
- // Otherwise the scorer is initialised
- vector<int> level; level.push_back(0);
-
- G4VPrimitiveScorer* DriftElectronMinosTPCScorer= new CylinderTPCScorers::PS_TPCAnode("DriftElectronsScore",level, 0) ;
-
- //and register it to the multifunctionnal detector
- m_MinosPadScorer->RegisterPrimitive(DriftElectronMinosTPCScorer);
-
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+////////////////////////////////////////////////////////////////
+void Minos::InitializeScorers() {
+ // This check is necessary in case the geometry is reloaded
+ bool already_exist = false;
+ bool already_exist2 = false;
+ bool already_exist3 = false;
+
+ m_MinosPadScorer = CheckScorer("MinosPadScorer",already_exist3) ;
+
+ if(already_exist && already_exist2 && already_exist3 )
+ return ;
+
+ // Otherwise the scorer is initialised
+ vector<int> level; level.push_back(0);
+
+ G4VPrimitiveScorer* DriftElectronMinosTPCScorer= new CylinderTPCScorers::PS_TPCAnode("DriftElectronsScore",level, 0) ;
+
+ //and register it to the multifunctionnal detector
+ m_MinosPadScorer->RegisterPrimitive(DriftElectronMinosTPCScorer);
+
G4VPrimitiveScorer* PadScorer= new CylinderTPCScorers::PS_TPCAnode("MinosTPCAnode",level, 0);
m_MinosPadScorer->RegisterPrimitive(PadScorer);
-
+
G4SDManager::GetSDMpointer()->AddNewDetector(m_MinosPadScorer) ;
}
diff --git a/NPSimulation/Detectors/Minos/Minos.hh b/NPSimulation/Detectors/Minos/Minos.hh
index 5c5db0349..71ca4078b 100644
--- a/NPSimulation/Detectors/Minos/Minos.hh
+++ b/NPSimulation/Detectors/Minos/Minos.hh
@@ -65,21 +65,12 @@ class Minos : public NPS::VDetector{
void DefineMaterials();
public:
- void SetTargetMaterial(G4String materialChoice);
- void SetChamberMaterial(G4String materialChoice);
- void SetTPCMaterial(G4String materialChoice);
- void SetWindowMaterial(G4String materialChoice);
- void SetInnerRohacellMaterial(G4String materialChoice);
- void SetOuterRohacellMaterial(G4String materialChoice);
- void SetKaptonMaterial(G4String materialChoice);
-
G4LogicalVolume* BuildSquareDetector();
G4LogicalVolume* BuildCylindricalDetector();
G4LogicalVolume* BuildTarget();
G4LogicalVolume* BuildChamber();
G4LogicalVolume* BuildInnerRohacell();
G4LogicalVolume* BuildOuterRohacell();
- G4LogicalVolume* BuildOuterOuterRohacell();
G4LogicalVolume* BuildKapton();
G4LogicalVolume* BuildOuterKapton();
G4LogicalVolume* BuildTPC();
@@ -97,12 +88,13 @@ class Minos : public NPS::VDetector{
G4Material* TargetMaterial;
G4Material* WindowMaterial;
G4Material* ChamberMaterial;
- G4Material* InnerRohacellMaterial;
- G4Material* OuterRohacellMaterial;
+ G4Material* RohacellMaterial;
G4Material* KaptonMaterial;
G4Material* TPCMaterial;
G4Material* defaultMaterial;
-
+ G4Material* Mylar;
+
+
G4LogicalVolume* m_SquareDetector;
G4LogicalVolume* m_CylindricalDetector;
@@ -123,21 +115,13 @@ class Minos : public NPS::VDetector{
G4LogicalVolume* logicWindow0;
// G4VPhysicalVolume* physiWindow0;
- G4Tubs* solidWindow1;
- G4LogicalVolume* logicWindow1;
- // G4VPhysicalVolume* physiWindow1;
-
- G4Tubs* solidWindow2;
- G4LogicalVolume* logicWindow2;
- // G4VPhysicalVolume* physiWindow2;
-
- G4Tubs* solidInnerRohacell;
- G4LogicalVolume* logicInnerRohacell;
+ G4Tubs* solidRohacell;
+ G4LogicalVolume* logicRohacell;
// G4VPhysicalVolume* physiInnerRohacell;
-
- G4Tubs* solidOuterRohacell;
- G4LogicalVolume* logicOuterRohacell;
- // G4VPhysicalVolume* physiOuterRohacell;
+
+ /* G4Tubs* solidOuterRohacell; */
+ /* G4LogicalVolume* logicOuterRohacell; */
+ /* // G4VPhysicalVolume* physiOuterRohacell; */
G4Tubs* solidKapton;
G4LogicalVolume* logicKapton;
@@ -203,9 +187,7 @@ class Minos : public NPS::VDetector{
G4VisAttributes* m_VisTarget;
G4VisAttributes* m_VissimpleBox;
G4VisAttributes* m_VisTPC;
- G4VisAttributes* m_VisInnerRohacell;
- G4VisAttributes* m_VisOuterRohacell;
- G4VisAttributes* m_VisOuterOuterRohacell;
+ G4VisAttributes* m_VisRohacell;
G4VisAttributes* m_VisKapton;
G4VisAttributes* m_VisTargetCell;
G4VisAttributes* m_VisOuterKapton;
diff --git a/Projects/Dali/Dali.detector b/Projects/Dali/Dali.detector
index 36f679550..b865ae207 100644
--- a/Projects/Dali/Dali.detector
+++ b/Projects/Dali/Dali.detector
@@ -1,97 +1,78 @@
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-Target
- THICKNESS= 20 mm
- RADIUS= 20 mm
- MATERIAL= CD2
- ANGLE= 0 deg
- X= 0 mm
- Y= 0 mm
- Z= 220 mm
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Dali
R = 212.4 mm
Alpha = 240 deg
- Zeta = -90.03 mm
- Shape = Square
- Material = NaI_Tl
+ Zeta = 298,28 mm
+ Ring = 1
+ Material = NaI(Tl)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Dali
R = 212.4 mm
Alpha = 180 deg
- Zeta = -90.03 mm
- Shape = Square
- Material = NaI_Tl
+ Zeta = 298,28 mm
+ Ring = 1
+ Material = NaI(Tl)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Dali
- R = 212.4 mm
+ R= 212.4 mm
Alpha = 120 deg
- Zeta = -90.03 mm
- Shape = Square
- Material = NaI_Tl
+ Zeta = 298,28 mm
+ Ring = 1
+ Material = NaI(Tl)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Dali
R = 212.4 mm
Alpha = 60 deg
- Zeta = -90.03 mm
- Shape = Square
- Material = NaI_Tl
+ Zeta = 298,28 mm
+ Ring = 1
+ Material = NaI(Tl)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Dali
R = 212.4 mm
Alpha = 0 deg
- Zeta = -90.03 mm
- Shape = Square
- Material = NaI_Tl
+ Zeta = 298,28 mm
+ Ring = 1
+ Material = NaI(Tl)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Dali
R = 212.4 mm
Alpha = -60 deg
- Zeta = -90.03 mm
- Shape = Square
- Material = NaI_Tl
-
-
-
+ Zeta = 298,28 mm
+ Ring = 1
+ Material = NaI(Tl)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Dali
R = 212.4 mm
Alpha = 240 deg
- Zeta = 90.03 mm
- Shape = Square
- Material = NaI_Tl
+ Zeta = 469.28 mm
+ Material = NaI(Tl)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Dali
R = 212.4 mm
Alpha = 180 deg
- Zeta = 90.03 mm
- Shape = Square
- Material = NaI_Tl
+ Zeta = 469.28 mm
+ Material = NaI(Tl)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Dali
R = 212.4 mm
Alpha = 120 deg
- Zeta = 90.03 mm
- Shape = Square
- Material = NaI_Tl
+ Zeta = 469.28 mm
+ Material = NaI(Tl)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Dali
R = 212.4 mm
Alpha = 60 deg
- Zeta = 90.03 mm
- Shape = Square
- Material = NaI_Tl
+ Zeta = 469.28 mm
+ Material = NaI(Tl)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Dali
R = 212.4 mm
Alpha = 0 deg
- Zeta = 90.03 mm
- Shape = Square
- Material = NaI_Tl
+ Zeta = 469.28 mm
+ Material = NaI(Tl)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Dali
R = 212.4 mm
Alpha = -60 deg
- Zeta = 90.03 mm
- Shape = Square
- Material = NaI_Tl
+ Zeta = 469.28 mm
+ Material = NaI(Tl)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/Projects/Dali/DaliMinos.detector b/Projects/Dali/DaliMinos.detector
old mode 100644
new mode 100755
--
GitLab