diff --git a/NPLib/Detectors/Dali/TDaliData.cxx b/NPLib/Detectors/Dali/TDaliData.cxx
index 542de349e870945908a19791576e294c234fefa6..70bc43ecbad125333353ae25e26f97e845fb84c2 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 43c26549f4c3d8fd256d839ddc66fe280d79b2d6..6c58a2fd80bb51adb1f546c557edf22edfbbdf4a 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 11df06300100125dee5587f1212aabe3933dfd70..72f870f57b62e15db01fae952424ee062c4f36ce 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 60d4a5cade4844fd193e7d1a579b3a909ffb6574..39a60b4e3a555d749dd99517edbbd59bdf81717b 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/CMakeLists.txt b/NPLib/Detectors/Minos/CMakeLists.txt
index 4e7173cabc758116b79916d94233c3fa5266f128..49ea047f815daef3723e7b599d040d6d910e6f3c 100644
--- a/NPLib/Detectors/Minos/CMakeLists.txt
+++ b/NPLib/Detectors/Minos/CMakeLists.txt
@@ -5,7 +5,7 @@ add_custom_command(OUTPUT TMinosDataDict.cxx COMMAND ../../scripts/build_dict.sh
add_custom_command(OUTPUT TMinosClustDict.cxx COMMAND ../../scripts/build_dict.sh TMinosClust.h TMinosClustDict.cxx TMinosClust.rootmap libNPMinos.dylib DEPENDS TMinosClust.h)
add_custom_command(OUTPUT TMinosResultDict.cxx COMMAND ../../scripts/build_dict.sh TMinosResult.h TMinosResultDict.cxx TMinosResult.rootmap libNPMinos.dylib DEPENDS TMinosResult.h)
-## Check if MINUIT2 is installed along with ROOT
+## Check if MINUIT is installed along with ROOT
find_library(libMinuit_FOUND NAMES Minuit HINTS "$ENV{ROOTSYS}/lib")
if(libMinuit_FOUND)
message(STATUS "Minuit support enabled for Minos.")
@@ -14,15 +14,18 @@ else()
message(STATUS "Minuit support disabled for Minos. Not found in $ENV{ROOTSYS}/lib")
endif()
-add_library(NPMinos SHARED TMinosClust.cxx TMinosData.cxx TMinosPhysics.cxx TMinosResult TMinosDataDict.cxx TMinosPhysicsDict.cxx TMinosClustDict.cxx TMinosResultDict.cxx)
if(libMinuit_FOUND)
-## Link to Minuit2 library
- target_link_libraries(NPMinos ${ROOT_LIBRARIES} NPCore NPTrackReconstruction Minuit)
+## Link to Minuit library
+ add_library(NPMinos SHARED TMinosClust.cxx TMinosData.cxx TMinosPhysics.cxx TMinosResult TMinosDataDict.cxx TMinosPhysicsDict.cxx TMinosClustDict.cxx TMinosResultDict.cxx)
+ target_link_libraries(NPMinos ${ROOT_LIBRARIES} Minuit NPCore NPTrackReconstruction )
+ install(FILES TMinosClust.h TMinosData.h TMinosPhysics.h TMinosResult.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
else()
-## No Minuit2 library, skip linking
- target_link_libraries(NPMinos ${ROOT_LIBRARIES} NPCore NPTrackReconstruction)
+## No Minuit library, skip linking
+ add_library(NPMinos SHARED TMinosData.cxx TMinosDataDict.cxx)
+ target_link_libraries(NPMinos ${ROOT_LIBRARIES} NPCore )
+ install(FILES TMinosData.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
+
endif()
-install(FILES TMinosClust.h TMinosData.h TMinosPhysics.h TMinosResult.h DESTINATION ${CMAKE_INCLUDE_OUTPUT_DIRECTORY})
diff --git a/NPLib/Detectors/Minos/TMinosPhysics.cxx b/NPLib/Detectors/Minos/TMinosPhysics.cxx
index 6435da703fedac6e4dc1207c9292b43feaeef8e6..ceda3da527066e0c569c0e9098d1c427e97c4427 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 23d9d57e5c6ed5739a27b8b959cf6628daed0325..2ac1d0cb9963adc00d1e9b74e5c63a6bff7f3acc 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/Physics/NPQFS.cxx b/NPLib/Physics/NPQFS.cxx
index 391ffa6c0571bf27e68b88a1156b365ab328e086..1b20bfbde3c40df7b0cb83bac4de05360e7e40d0 100644
--- a/NPLib/Physics/NPQFS.cxx
+++ b/NPLib/Physics/NPQFS.cxx
@@ -7,29 +7,32 @@
/*****************************************************************************
* *
- * Original Author : F. Flavigny contact address: flavigny@ipno.in2p3.fr *
+ * Original Author : F. Flavigny contact address: flavigny@lpccaen.in2p3.fr *
* *
* Creation Date : April 2019 *
- * Last update : April 2019 *
+ * Last update : Nov 2019 *
*---------------------------------------------------------------------------*
* Decription: *
* This class deal with Quasi Free Scattering Reaction in which a cluster *
* or a nucleon is removed from a projectile by interaction with a target *
* nucleon (proton target in general) *
* *
- * Labeling is: *
- * A --> i ==> B + (i -> a) => B + 1 + 2 *
+ * First step (dissociation): A -> B + c *
+ * Second step (scattering) : c + T -> 1 + 2 *
+ * Labeling is: *
+ * *
+ * A --> T ==> B + (c -> T) => B + 1 + 2 *
* *
* where: *
* + A is the beam nucleus *
- * + i is the target nucleon (proton) *
+ * + T is the target nucleon (proton) *
* *
* + B is the residual fragment (beam-like) *
- * + 1 is the scattered target nucleon (former i) *
- * + 2 is the knocked-out cluster/nucleon (noted a) in the intermediate *
+ * + 1 is the scattered target nucleon (former T) *
+ * + 2 is the knocked-out cluster/nucleon (noted c) in the intermediate *
*---------------------------------------------------------------------------*
* Comment: *
- * + *
+ * + Adapted from original event generator from V. Panin (R3B collab) *
* *
*****************************************************************************/
@@ -55,16 +58,625 @@ using namespace NPUNITS;
// ROOT
#include"TF1.h"
+#include"TRandom3.h"
ClassImp(QFS)
-/*
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+////////////////////////////////////////////////////////////////////////////////
+
QFS::QFS(){
- fVerboseLevel=0;
+
+ //------------- Default Constructor -------------
+ fVerboseLevel = NPOptionManager::getInstance()->GetVerboseLevel();
+ fBeamEnergy = 0;
+ fThetaCM = 0;
+ fPhiCM = 0;
+
+ fExcitationA = 0;
+ fExcitationB = 0;
+ fMomentumSigma = 0;
+ fshootB=false;
+ fshoot1=true;
+ fshoot2=true;
+ fisotropic = true;
+
+ fTheta2VsTheta1 = 0;
+ fPhi2VsPhi1 = 0;
}
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+////////////////////////////////////////////////////////////////////////////////
+
QFS::~QFS(){
+
+}
+
+/////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+void QFS::ReadConfigurationFile(string Path){
+ ifstream ReactionFile;
+ string GlobalPath = getenv("NPTOOL");
+ string StandardPath = GlobalPath + "/Inputs/EventGenerator/" + Path;
+ ReactionFile.open(Path.c_str());
+ if (!ReactionFile.is_open()) {
+ ReactionFile.open(StandardPath.c_str());
+ if(ReactionFile.is_open()) {
+ Path = StandardPath;
+ }
+ else {cout << "QFS File " << Path << " not found" << endl;exit(1);}
+ }
+ NPL::InputParser parser(Path);
+ ReadConfigurationFile(parser);
+}
+
+////////////////////////////////////////////////////////////////////////////////
+
+void QFS::ReadConfigurationFile(NPL::InputParser parser){
+
+ cout << " In QFS ReadConfiguration " << endl;
+ vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("QFSReaction");
+ if(blocks.size()>0 && NPOptionManager::getInstance()->GetVerboseLevel())
+ cout << endl << "\033[1;35m//// QFS reaction found " << endl;
+
+ vector<string> token1 = {"Beam","Target","Scattered","KnockedOut","Heavy"};
+ for(unsigned int i = 0 ; i < blocks.size() ; i++){
+ if(blocks[i]->HasTokenList(token1)){
+ int v = NPOptionManager::getInstance()->GetVerboseLevel();
+ NPOptionManager::getInstance()->SetVerboseLevel(0);
+ fNucleiA.ReadConfigurationFile(parser);
+ NPOptionManager::getInstance()->SetVerboseLevel(v);
+
+ fBeamEnergy= fNucleiA.GetEnergy();
+ GetNucleus(blocks[i]->GetString("Beam"),parser);
+ fNucleiT = GetNucleus(blocks[i]->GetString("Target"),parser);
+ fNucleiB = GetNucleus(blocks[i]->GetString("Heavy"),parser);
+ fNuclei1 = GetNucleus(blocks[i]->GetString("Scattered"),parser);
+ fNuclei2 = GetNucleus(blocks[i]->GetString("KnockedOut"),parser);
+ }
+ else{
+ cout << "ERROR: check your input file formatting \033[0m" << endl;
+ exit(1);
+ }
+ if(blocks[i]->HasToken("ExcitationEnergyBeam")){
+ fExcitationA = blocks[i]->GetDouble("ExcitationEnergyBeam","MeV");
+ }
+ if(blocks[i]->HasToken("ExcitationEnergyHeavy")){
+ fExcitationB = blocks[i]->GetDouble("ExcitationEnergyHeavy","MeV");
+ }
+ if(blocks[i]->HasToken("MomentumSigma")){
+ fMomentumSigma = blocks[i]->GetDouble("MomentumSigma","MeV");
+ }
+ if(blocks[i]->HasToken("ShootHeavy")){
+ fshootB = blocks[i]->GetInt("ShootHeavy");
+ }
+ if(blocks[i]->HasToken("ShootLight")){
+ fshoot1 = blocks[i]->GetInt("ShootLight");
+ fshoot2 = blocks[i]->GetInt("ShootLight");
+ }
+ }
+
+ cout << "\033[0m" ;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+Nucleus QFS::GetNucleus(string name, NPL::InputParser parser){
+
+ vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithTokenAndValue("DefineNucleus",name);
+ unsigned int size = blocks.size();
+ if(size==0)
+ return NPL::Nucleus(name);
+ else if(size==1){
+ cout << " -- User defined nucleus " << name << " -- " << endl;
+ vector<string> token = {"SubPart","BindingEnergy"};
+ if(blocks[0]->HasTokenList(token)){
+ NPL::Nucleus N(name,blocks[0]->GetVectorString("SubPart"),blocks[0]->GetDouble("BindingEnergy","MeV"));
+ if(blocks[0]->HasToken("ExcitationEnergy"))
+ N.SetExcitationEnergy(blocks[0]->GetDouble("ExcitationEnergy","MeV"));
+ if(blocks[0]->HasToken("SpinParity"))
+ N.SetSpinParity(blocks[0]->GetString("SpinParity").c_str());
+ if(blocks[0]->HasToken("Spin"))
+ N.SetSpin(blocks[0]->GetDouble("Spin",""));
+ if(blocks[0]->HasToken("Parity"))
+ N.SetParity(blocks[0]->GetString("Parity").c_str());
+ if(blocks[0]->HasToken("LifeTime"))
+ N.SetLifeTime(blocks[0]->GetDouble("LifeTime","s"));
+
+ cout << " -- -- -- -- -- -- -- -- -- -- --" << endl;
+ return N;
+ }
+ }
+ else{
+ NPL::SendErrorAndExit("NPL::QFS","Too many nuclei define with the same name");
+ }
+}
+
+
+////////////////////////////////////////////////////////////////////////////////////////////
+void QFS::CalculateVariables(){
+
+ if(fBeamEnergy < 0)
+ fBeamEnergy = 0 ;
+
+ //cout<<"---- COMPUTE ------"<<endl;
+ // cout<<"--CM--"<<endl;
+
+ mA = fNucleiA.Mass(); // Beam mass in MeV
+ mT = fNucleiT.Mass(); // Target mass in MeV
+ mB = fNucleiB.Mass(); // Heavy residual mass in MeV
+ m1 = mT; // scattered target nucleon (same mass);
+ m2 = fNuclei2.Mass(); // knocked out cluster mass in MeV
+ ma = m2; // intermediate cluster mass in MeV (same);
+
+ double TA = fBeamEnergy; // Beam kinetic energy
+ double PA = sqrt(TA*(TA+2*mA)); // Beam momentum (norm)
+ double EA = sqrt(mA*mA + PA*PA); // Beam total energy
+ fEnergyImpulsionLab_A = TLorentzVector(0.,0.,PA,EA);
+
+ //Internal momentum of removed cluster/nucleon
+ static TRandom3 r;
+ //r.SetSeed(0);
+ //double mom_sigma = 0; // MeV/c
+ Pa.SetX(r.Gaus(0.,fMomentumSigma));
+ Pa.SetY(r.Gaus(0.,fMomentumSigma));
+ Pa.SetZ(r.Gaus(0.,fMomentumSigma));
+
+ //Internal momentum of heavy recoil after removal
+ PB.SetXYZ( (-Pa.X()) , (-Pa.Y()) , (-Pa.Z()) );
+
+ //cout<<"Pa_cm=\t("<<Pa.Px()<<","<<Pa.Py()<<","<<Pa.Pz()<<")"<<endl;
+ //cout<<"||Pa||^2=\t("<<Pa.Mag2()<<endl;
+ //cout<<"mA^2=\t"<<mA*mA<<endl;
+ //cout<<"mB^2=\t"<<mB*mB<<endl;
+
+
+ // Off-shell mass of the bound nucleon from E conservation
+ // in virtual dissociation of A -> B + a
+ double buffer = mA*mA + mB*mB - 2*mA*sqrt(mB*mB+Pa.Mag2()) ;
+ if(buffer<=0) { cout<<"ERROR off shell mass ma_off=\t"<<buffer<<endl; return;}
+ ma_off = sqrt(buffer);
+
+ //deduced total energies of "a" and "B" in restframe of A
+ double Ea = sqrt(ma_off*ma_off + Pa.Mag2());
+ double EB = sqrt(mB*mB + PB.Mag2());
+
+ //cout<<"ma_off^2=\t"<<buffer<<endl;
+ //cout<<"Ea=\t"<<Ea<<endl;
+ //cout<<"EB=\t"<<EB<<endl;
+
+ fEnergyImpulsionLab_a = TLorentzVector(Pa,Ea);
+ fEnergyImpulsionLab_B = TLorentzVector(PB,EB);
+ fEnergyImpulsionLab_a.Boost(0,0,fEnergyImpulsionLab_A.Beta());
+ fEnergyImpulsionLab_B.Boost(0,0,fEnergyImpulsionLab_A.Beta());
+ Ea_lab = fEnergyImpulsionLab_a.E();
+ EB_lab = fEnergyImpulsionLab_B.E();
+ Pa = fEnergyImpulsionLab_a.Vect();
+ PB = fEnergyImpulsionLab_B.Vect();
+
+ // Scattering part (2-body kinematics)
+ // virtual cluster of mass "ma_off" scattering on target T
+ // to give scattered cluster with real mass (ma=m2)
+ // and scattered target (mT=m1)
+
+ fQValue =ma_off+mT-m1-m2;
+
+ s = ma_off*ma_off + mT*mT + 2*mT*Ea_lab ;
+ fTotalEnergyImpulsionCM = TLorentzVector(0,0,0,sqrt(s));
+ fEcm = sqrt(s) - m1 -m2;
+ if(fEcm<=0) { cout<<"ERROR Ecm negative =\t"<<fEcm<<endl; return;}
+
+ ECM_a = (s + ma_off*ma_off - mT*mT)/(2*sqrt(s));
+ ECM_T = (s + mT*mT - ma_off*ma_off)/(2*sqrt(s));
+ ECM_1 = (s + m1*m1 - m2*m2)/(2*sqrt(s));
+ ECM_2 = (s + m2*m2 - m1*m1)/(2*sqrt(s));
+
+ pCM_a = sqrt(ECM_a*ECM_a - ma_off*ma_off);
+ pCM_T = sqrt(ECM_T*ECM_T - mT*mT);
+ pCM_1 = sqrt(ECM_1*ECM_1 - m1*m1);
+ pCM_2 = sqrt(ECM_2*ECM_2 - m2*m2);
+
+ BetaCM = Pa.Mag() / (Ea_lab + mT);
+
+ //if(ECM_a<0 || ECM_T<0 || ECM_1<0 || ECM_2<0 ) {
+ /*
+ if(pCM_a<0 || pCM_T<0 || pCM_1<0 || pCM_2<0 ) {
+ cout<<"--LAB after Boost--"<<endl;
+ cout<<"Pa_lab=\t("<<Pa.Px()<<","<<Pa.Py()<<","<<Pa.Pz()<<")"<<endl;
+ cout<<"PB_lab=\t("<<PB.Px()<<","<<PB.Py()<<","<<PB.Pz()<<")"<<endl;
+ cout<<"Ea_lab=\t"<<Ea_lab<<endl;
+ cout<<"EB_lab=\t"<<EB_lab<<endl;
+ cout<<"--Back to CM--"<<endl;
+ cout<<"fQValue=\t"<<fQValue<<endl;
+ cout<<"s=\t"<<s<<endl;
+ cout<<"Ecm=\t"<<fEcm<<endl;
+ cout<<"ea*=\t"<<ECM_a<<endl;
+ cout<<"eT*=\t"<<ECM_T<<endl;
+ cout<<"e1*=\t"<<ECM_1<<endl;
+ cout<<"p1*=\t"<<pCM_1<<endl;
+ cout<<"e2*=\t"<<ECM_2<<endl;
+ cout<<"p2*=\t"<<pCM_2<<endl;
+ cout<<"beta_cm=\t"<<BetaCM<<endl;
+ }
+ */
+}
+
+////////////////////////////////////////////////////////////////////////////////////////////
+
+void QFS::KineRelativistic(double& ThetaLab1, double& PhiLab1, double& KineticEnergyLab1, double& ThetaLab2, double& PhiLab2, double& KineticEnergyLab2){
+
+ CalculateVariables();
+
+ //double thetaCM_1 = fThetaCM*TMath::Pi()/180.;
+ double thetaCM_1 = fThetaCM;
+ double thetaCM_2 = M_PI - thetaCM_1;
+ //double phiCM_1 = fPhiCM*TMath::Pi()/180.;
+ double phiCM_1 = fPhiCM;
+ double phiCM_2 = 2*M_PI - phiCM_1;
+
+ TVector3 z_axis(0.,0.,1.);
+
+ fEnergyImpulsionCM_2 = TLorentzVector(
+ pCM_2*sin(thetaCM_2)*cos(phiCM_2),
+ pCM_2*cos(thetaCM_2)*sin(phiCM_2),
+ pCM_2*cos(thetaCM_2),
+ ECM_2);
+ fEnergyImpulsionCM_1 = fTotalEnergyImpulsionCM - fEnergyImpulsionCM_2;
+
+ fEnergyImpulsionLab_1 = fEnergyImpulsionCM_1;
+ fEnergyImpulsionLab_1.Boost(0,0,BetaCM);
+ fEnergyImpulsionLab_2 = fEnergyImpulsionCM_2;
+ fEnergyImpulsionLab_2.Boost(0,0,BetaCM);
+
+ // Angle in the lab frame
+ //ThetaLab1 = fEnergyImpulsionLab_1.Angle(fEnergyImpulsionLab_A.Vect());
+ ThetaLab1 = fEnergyImpulsionLab_1.Angle(z_axis);
+ if (ThetaLab1 < 0) ThetaLab1 += M_PI;
+ //ThetaLab2 = fEnergyImpulsionLab_4.Angle(fEnergyImpulsionLab_A.Vect());
+ ThetaLab2 = fEnergyImpulsionLab_2.Angle(z_axis);
+ if (fabs(ThetaLab1) < 1e-6) ThetaLab1 = 0;
+ ThetaLab2 = fabs(ThetaLab2);
+ if (fabs(ThetaLab2) < 1e-6) ThetaLab2 = 0;
+
+ PhiLab1 = M_PI + fEnergyImpulsionLab_1.Vect().Phi();
+ if (fabs(PhiLab1) < 1e-6) PhiLab1 = 0;
+ PhiLab2 = M_PI + fEnergyImpulsionLab_2.Vect().Phi();
+ if (fabs(PhiLab2) < 1e-6) PhiLab2 = 0;
+
+ // Kinetic Energy in the lab frame
+ KineticEnergyLab1 = fEnergyImpulsionLab_1.E() - m1;
+ KineticEnergyLab2 = fEnergyImpulsionLab_2.E() - m2;
+
+ // test for total energy conversion
+ //if (fabs(fTotalEnergyImpulsionLab.E() - (fEnergyImpulsionLab_1.E()+fEnergyImpulsionLab_2.E())) > 1e-6)
+ // cout << "Problem for energy conservation" << endl;
+
+/*
+ cout<<"--KINE RELATIVISTIC--"<<endl;
+ cout<<"theta_cm:"<<fThetaCM*180./TMath::Pi()<<endl;
+ cout<<"phi_cm:"<<fPhiCM*180./TMath::Pi()<<endl;
+ cout<<"P1_CM=\t("<<fEnergyImpulsionCM_1.Px()<<","<<fEnergyImpulsionCM_1.Py()<<","<<fEnergyImpulsionCM_1.Pz()<<")"<<endl;
+ cout<<"P2_CM=\t("<<fEnergyImpulsionCM_2.Px()<<","<<fEnergyImpulsionCM_2.Py()<<","<<fEnergyImpulsionCM_2.Pz()<<")"<<endl;
+ cout<<"P1_lab=\t("<<fEnergyImpulsionLab_1.Px()<<","<<fEnergyImpulsionLab_1.Py()<<","<<fEnergyImpulsionLab_1.Pz()<<")"<<endl;
+ cout<<"P2_lab=\t("<<fEnergyImpulsionLab_2.Px()<<","<<fEnergyImpulsionLab_2.Py()<<","<<fEnergyImpulsionLab_2.Pz()<<")"<<endl;
+ cout<<"Theta1:\t"<<fEnergyImpulsionLab_1.Vect().Theta()*180./TMath::Pi()<<endl;
+ cout<<"Theta2:\t"<<fEnergyImpulsionLab_2.Vect().Theta()*180./TMath::Pi()<<endl;
+ cout<<"Phi1:\t"<<M_PI+fEnergyImpulsionLab_1.Vect().Phi()*180./TMath::Pi()<<endl;
+ cout<<"Phi2:\t"<<M_PI+fEnergyImpulsionLab_2.Vect().Phi()*180./TMath::Pi()<<endl;
+*/
+
}
+
+////////////////////////////////////////////////////////////////////////////////////////////
+
+double QFS::ShootRandomThetaCM(){
+
+ double theta; // CM
+/*
+ if(fDoubleDifferentialCrossSectionHist){
+ // Take a slice in energy
+ TAxis* Y = fDoubleDifferentialCrossSectionHist->GetYaxis();
+ int binY;
+
+ // Those test are there for the tail event of the energy distribution
+ // In case the energy is outside the range of the 2D histo we take the
+ // closest available CS
+ if(Y->FindBin(fBeamEnergy) > Y->GetLast())
+ binY = Y->GetLast()-1;
+
+ else if(Y->FindBin(fBeamEnergy) < Y->GetFirst())
+ binY = Y->GetFirst()+1;
+
+ else
+ binY = Y->FindBin(fBeamEnergy);
+
+ TH1D* Proj = fDoubleDifferentialCrossSectionHist->ProjectionX("proj",binY,binY);
+ SetThetaCM( theta=Proj->GetRandom()*deg );
+ }
+ else if (fLabCrossSection){
+ double thetalab=-1;
+ double energylab=-1;
+ while(energylab<0){
+ thetalab=fCrossSectionHist->GetRandom()*deg; //shoot in lab
+ energylab=EnergyLabFromThetaLab(thetalab); //get corresponding energy
+ }
+ theta = EnergyLabToThetaCM(energylab, thetalab); //transform to theta CM
+ SetThetaCM( theta );
+ }
+ else{
+ // When root perform a Spline interpolation to shoot random number out of
+ // the distribution, it can over shoot and output a number larger that 180
+ // this lead to an additional signal at 0-4 deg Lab, especially when using a
+ // flat distribution.
+ // This fix it.
+ theta=181;
+ if(theta/deg>180)
+ theta=fCrossSectionHist->GetRandom();
+ //cout << " Shooting Random ThetaCM " << theta << endl;
+ SetThetaCM( theta*deg );
+ }
*/
+ return theta;
+}
+
+////////////////////////////////////////////////////////////////////////////////////////////
+
+TGraph* QFS::GetTheta2VsTheta1(double AngleStep_CM){
+
+ vector<double> vx;
+ vector<double> vy;
+ double theta1,phi1,E1,theta2,phi2,E2;
+
+ for (double angle=0 ; angle <= 180 ; angle+=AngleStep_CM){
+ SetThetaCM(angle*TMath::Pi()/180.);
+ KineRelativistic(theta1, phi1, E1, theta2, phi2, E2);
+
+ vx.push_back(theta1*180./M_PI);
+ vy.push_back(theta2*180./M_PI);
+ }
+ fTheta2VsTheta1 = new TGraph(vx.size(),&vx[0],&vy[0]);
+
+ return(fTheta2VsTheta1);
+}
+
+////////////////////////////////////////////////////////////////////////////////////////////
+
+TGraph* QFS::GetPhi2VsPhi1(double AngleStep_CM){
+
+ vector<double> vx;
+ vector<double> vy;
+ double theta1,phi1,E1,theta2,phi2,E2;
+
+ for (double theta=0 ; theta <= 180 ; theta+=AngleStep_CM){
+ for (double angle=0 ; angle <= 180 ; angle+=AngleStep_CM){
+ SetThetaCM(theta*TMath::Pi()/180.);
+ SetPhiCM(angle*TMath::Pi()/180.);
+ KineRelativistic(theta1, phi1, E1, theta2, phi2, E2);
+ vx.push_back(phi1*180./M_PI);
+ vy.push_back(phi2*180./M_PI);
+ }
+ }
+ fPhi2VsPhi1 = new TGraph(vx.size(),&vx[0],&vy[0]);
+
+ return(fPhi2VsPhi1);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// Check whenever the reaction is allowed at a given energy
+bool QFS::IsAllowed(){//double Energy){
+ //double AvailableEnergy = Energy + fNuclei1.Mass() + fNuclei2.Mass();
+ //double RequiredEnergy = fNuclei3.Mass() + fNuclei4.Mass();
+
+ //if(AvailableEnergy>RequiredEnergy)
+ return true;
+ //else
+ // return false;
+}
+
+////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////
+///////////////// Old R3B method not using TLorentz Vector ////////////////////////////////
+/////////////////// (used as a reference)///////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////
+/*
+void QFS::CalculateVariablesOld(){
+
+ initializePrecomputeVariable();
+
+ vector<double> theta1;
+ vector<double> theta2;
+ vector<double> phi1;
+ vector<double> phi2;
+
+ //for(int i=0; i<=180; i++){
+ int i = 29;
+ KineR3B(s, ma_off, mT, ma, (double)i);
+ if(!good) { cout<<"ERROR CM calculations!!!"<<endl; return;}
+
+ TVector3 P1_cm(0.,0.,1.), P2_cm(0.,0.,1.);
+ P2_cm.SetMag(p_clust);
+ P2_cm.SetTheta(theta_clust);
+ TRandom3 ra;
+ ra.SetSeed(0);
+ //P2_cm.SetPhi(ra.Uniform(-1*TMath::Pi(),+1*TMath::Pi()));
+ P2_cm.SetPhi(0.);
+ P1_cm.SetX(-P2_cm.X());
+ P1_cm.SetY(-P2_cm.Y());
+ P1_cm.SetZ(-P2_cm.Z());
+ cout<<"FFFFFFFFFFFFFFFFF"<<endl;
+ cout<<"P1_CM=\t("<<P1_cm.X()<<","<<P1_cm.Y()<<","<<P1_cm.Z()<<")"<<endl;
+ cout<<"P2_CM=\t("<<P2_cm.X()<<","<<P2_cm.Y()<<","<<P2_cm.Z()<<")"<<endl;
+
+
+ // Calculate relative to direction of quasi-particle (cluster)
+
+ double beta_cm = -Pa.Mag() / (Ea_lab + mT);
+ double gamma_cm = 1/sqrt(1-beta_cm*beta_cm);
+
+ pair<double,double> lor_a1 = Lorentz(gamma_cm,beta_cm,e_scat,P1_cm.Z());
+ pair<double,double> lor_a2 = Lorentz(gamma_cm,beta_cm,e_clust,P2_cm.Z());
+
+ P1_cm.SetZ(lor_a1.second);
+ P2_cm.SetZ(lor_a2.second);
+
+ //Rotating back to beam direction
+ //TVector3 P1_L = Rotations(P1_cm, Pa);
+ //TVector3 P2_L = Rotations(P2_cm, Pa);
+
+ TVector3 P1_L = P1_cm;
+ TVector3 P2_L = P2_cm;
+ TVector3 direction = Pa.Unit();
+ //P1_L.RotateUz(direction);
+ //P1_L.RotateUz(direction);
+
+ cout<<"----Calculate variables output------"<<endl;
+ cout<<"--CM--"<<endl;
+ cout<<"theta1*=\t"<<theta_scat*180/TMath::Pi()<<endl;
+ cout<<"theta2*=\t"<<theta_clust*180/TMath::Pi()<<endl;
+ cout<<"e1*=\t"<<e_scat<<endl;
+ cout<<"p1*=\t"<<p_scat<<endl;
+ cout<<"e2*=\t"<<e_clust<<endl;
+ cout<<"p2*=\t"<<p_clust<<endl;
+ cout<<"T=\t"<<T<<endl;
+ cout<<"beta_cm=\t"<<beta_cm<<endl;
+ cout<<"gamma_cm=\t"<<gamma_cm<<endl;
+
+ cout<<"--LAB (cluster dir)--"<<endl;
+ cout<<"P1_lab=\t("<<P1_cm.X()<<","<<P1_cm.Y()<<","<<P1_cm.Z()<<")"<<endl;
+ cout<<"P2_lab=\t("<<P2_cm.X()<<","<<P2_cm.Y()<<","<<P2_cm.Z()<<")"<<endl;
+ cout<<"Theta1:\t"<<P1_cm.Theta()*180./TMath::Pi()<<endl;
+ cout<<"Theta2:\t"<<P2_cm.Theta()*180./TMath::Pi()<<endl;
+
+ cout<<"--LAB--"<<endl;
+ cout<<"Theta1L:\t"<<P1_L.Theta()*180./TMath::Pi()<<endl;
+ cout<<"Theta2L:\t"<<P2_L.Theta()*180./TMath::Pi()<<endl;
+ cout<<"Phi1L:\t"<<P1_L.Phi()*180./TMath::Pi()<<endl;
+ cout<<"Phi2L:\t"<<P2_L.Phi()*180./TMath::Pi()<<endl;
+
+ //cout<<P1_cm.Theta()*180./TMath::Pi()<<"\t"<<P2_cm.Theta()*180./TMath::Pi()<<endl;
+ //cout<<P1_L.Phi()*180./TMath::Pi()<<"\t"<<P2_L.Phi()*180./TMath::Pi()<<endl;
+
+ theta1.push_back(P1_L.Theta()*180./TMath::Pi());
+ theta2.push_back(P2_L.Theta()*180./TMath::Pi());
+
+ double temp_phi1 = P1_L.Phi();
+ double temp_phi2 = P2_L.Phi();
+ phi1.push_back(180. + P1_L.Phi()*180./TMath::Pi());
+ phi2.push_back(180. + P2_L.Phi()*180./TMath::Pi());
+
+ //}
+ TGraph* fTheta2VsTheta1 = new TGraph(theta1.size(),&theta1[0],&theta2[0]);
+ TGraph* fPhi2VsPhi1 = new TGraph(phi1.size(),&phi1[0],&phi2[0]);
+ fTheta2VsTheta1->SetName("Theta2VsTheta1");
+ fPhi2VsPhi1->SetName("Phi2VsPhi1");
+ TFile* f = new TFile("graphs.root","RECREATE");
+ fTheta2VsTheta1->Write();
+ fPhi2VsPhi1->Write();
+ f->Close();
+
+
+}
+
+// Calculate elastic scattering kinematics in CM-system (1-target proton, 2-cluster)
+void QFS::KineR3B(double s,double m2off,double m1,double m2,double thetacm)
+{
+ if(thetacm>180 || thetacm<0){
+ cout << "\nERROR! ThetaCM (in deg) should be between 0 and 180"<<endl;
+ return;
+ }
+
+ e_clust = 0;
+ p_clust = 0;
+ theta_clust = 0;
+ e_scat = 0;
+ p_scat = 0;
+ theta_scat = 0;
+ T = 0;
+ good = false;
+
+
+ double X = s;
+ double Y = m2off*m2off;
+ double Z = m1*m1;
+ double sqrt_s = sqrt(s);
+
+ // Kinematics before the scattering process
+ // (with one off-shell mass)
+ double p2_off = sqrt(function(X,Y,Z))/2/sqrt_s;
+ double p1_off = p2_off;
+ // CM energies
+ double e1_off = (s+Z-Y)/2/sqrt_s;
+ double e2_off = (s+Y-Z)/2/sqrt_s;
+
+ // Now take the real masses (after scattering)
+ Y = m2*m2; Z = m1*m1;
+ //And check whether the kinematical function is ok
+ //for this specific kinematical case
+ double ERROR_CI = function(X,Y,Z);
+ if(ERROR_CI <= 0.){
+ cout << "\nERROR!!! Kinematical function is negative!";
+ return;
+ }
+
+ // Kinematics after the scattering process
+ // (with all real masses)
+ double p2 = sqrt(function(X,Y,Z))/2/sqrt_s;
+ double p1 = p2;
+ double e1 = (s+Z-Y)/2/sqrt_s;
+ double e2 = (s+Y-Z)/2/sqrt_s;
+
+ // Let's consider momentum transfer <t> from the
+ // target particle 1 to the cluster 2
+ double t = 2*(m1*m1 - e1_off*e1 + p1_off*p1*cos(thetacm*TMath::Pi()/180.));
+
+ //CM scattering angles
+ double theta1 = thetacm*TMath::Pi()/180.;
+ double theta2 = TMath::Pi() - theta1;
+
+ e_clust = e2;
+ p_clust = p2;
+ theta_clust = theta2;
+
+ e_scat = e1;
+ p_scat = p1;
+ theta_scat = theta1;
+
+ T = t;
+ good = true;
+
+
+}
+
+
+
+double QFS::function(double x,double y,double z)
+{
+ double lambda = x*x + y*y + z*z - 2*x*y - 2*x*z - 2*y*z;
+ return lambda;
+}
+
+//---- Two consecutive rotations
+//first around Z on <phi>, then around new X' on <theta> (1=Pcm, 2=Pa in lab)
+TVector3 QFS::Rotations(TVector3 v1,TVector3 v2)
+{
+ double CT = v2.Z()/v2.Mag(); // cos(theta) of v2 wrt. Z-axis
+ double ST = sqrt(1-CT*CT); // sin(theta)
+ double CF = v2.X()/v2.Mag()/ST;
+ double SF = v2.Y()/v2.Mag()/ST;
+
+ TVector3 v3;
+ double _v3x = v1.X()*CT*CF - v1.Y()*SF + v1.Z()*ST*CF;
+ double _v3y = v1.X()*CT*SF + v1.Y()*CF + v1.Z()*ST*SF;
+ double _v3z = -v1.X()*ST + v1.Z()*CT;
+ v3.SetXYZ(_v3x,_v3y,_v3z);
+ return v3;
+}
+
+
+
+pair<double, double> QFS::Lorentz(double gamma,double beta,double e,double p)
+{
+ double eL = gamma*e - gamma*beta*p;
+ double pL = gamma*p - gamma*beta*e;
+ return make_pair(eL, pL);
+}
+*/
diff --git a/NPLib/Physics/NPQFS.h b/NPLib/Physics/NPQFS.h
index b0107b6c85c3def61911987017ee613fa918ed36..0d7f63c753157247acd905f10497bcd11f7eb4a7 100644
--- a/NPLib/Physics/NPQFS.h
+++ b/NPLib/Physics/NPQFS.h
@@ -1,6 +1,5 @@
#ifndef NPQFS_h
#define NPQFS_h
-/*****************************************************************************/
/*****************************************************************************
* Copyright (C) 2009-2019 this file is part of the NPTool Project *
* *
@@ -10,29 +9,32 @@
/*****************************************************************************
* *
- * Original Author : F. Flavigny contact address: flavigny@ipno.in2p3.fr *
+ * Original Author : F. Flavigny contact address: flavigny@lpccaen.in2p3.fr *
* *
* Creation Date : April 2019 *
- * Last update : April 2019 *
+ * Last update : Nov 2019 *
*---------------------------------------------------------------------------*
* Decription: *
* This class deal with Quasi Free Scattering Reaction in which a cluster *
* or a nucleon is removed from a projectile by interaction with a target *
* nucleon (proton target in general) *
* *
- * Labeling is: *
- * A --> i ==> B + (i -> a) => B + 1 + 2 *
+ * First step (dissociation): A -> B + c *
+ * Second step (scattering) : c + T -> 1 + 2 *
+ * Labeling is: *
+ * *
+ * A --> T ==> B + (c -> T) => B + 1 + 2 *
* *
* where: *
* + A is the beam nucleus *
- * + i is the target nucleon (proton) *
+ * + T is the target nucleon (proton) *
* *
* + B is the residual fragment (beam-like) *
- * + 1 is the scattered target nucleon (former i) *
- * + 2 is the knocked-out cluster/nucleon (noted a) in the intermediate *
+ * + 1 is the scattered target nucleon (former T) *
+ * + 2 is the knocked-out cluster/nucleon (noted c) in the intermediate *
*---------------------------------------------------------------------------*
* Comment: *
- * + *
+ * + Adapted from original event generator from V. Panin (R3B collab) *
* *
*****************************************************************************/
@@ -41,6 +43,7 @@
// NPL
#include "NPNucleus.h"
+#include "NPBeam.h"
#include "NPInputParser.h"
using namespace NPL;
@@ -60,12 +63,138 @@ namespace NPL{
class QFS{
public: // Constructors and Destructors
- QFS(){};
- ~QFS(){};
+ QFS();
+ ~QFS();
private:
int fVerboseLevel;
+ Beam fNucleiA; // Beam (A)
+ Nucleus fNucleiT; // Target (T)
+ Nucleus fNucleiB; // Beam-like ejectile (B)
+ Nucleus fNuclei1; // Target-like ejectile (1)
+ Nucleus fNuclei2; // Knocked-out nucleon/cluster (2)
+ double fQValue; // Q-value in MeV
+ double fEcm; // Ecm in MeV
+ double fThetaCM; // Center-of-mass theta angle in radian
+ double fPhiCM; // Center-of-mass Phi angle in radian
+ double fBeamEnergy; // Beam energy in MeV
+ double fExcitationA; // Excitation energy in MeV of the beam, useful for isomers
+ double fExcitationB; // Excitation energy in MeV of beam-like ejectile
+ double fMomentumSigma; // Width of the momentum distribution of the removed cluster (sigma)
+ double fisotropic;
+
+ TGraph* fTheta2VsTheta1;
+ TGraph* fPhi2VsPhi1;
+
+ // used for MC simulations
+ bool fshootB; // shoot beam-like ejectile
+ bool fshoot1; // shoot light ejectile &
+ bool fshoot2; // shoot light ejectile 2
+
+ public:
+ Nucleus GetNucleus(string name, NPL::InputParser parser);
+ void ReadConfigurationFile(string Path);
+ void ReadConfigurationFile(NPL::InputParser);
+ void CalculateVariables();
+ void KineRelativistic(double &ThetaLab1, double &PhiLab1, double &KineticEnergyLab1, double &ThetaLab2, double &PhiLab2, double &KineticEnergyLab2);
+ double ShootRandomThetaCM();
+ bool IsAllowed();
+
+ private: // intern precompute variable
+ double mA;
+ double mB;
+ double ma_off;
+ double ma;
+ double mT;
+ double m1;
+ double m2;
+
+ TVector3 Pa, PB;
+ double Ea_lab,EB_lab;
+
+ // Lorentz Vector
+ TLorentzVector fEnergyImpulsionLab_A;
+ TLorentzVector fEnergyImpulsionLab_B;
+ TLorentzVector fEnergyImpulsionLab_a;
+ TLorentzVector fEnergyImpulsionLab_T;
+ TLorentzVector fEnergyImpulsionLab_1;
+ TLorentzVector fEnergyImpulsionLab_2;
+ TLorentzVector fTotalEnergyImpulsionLab;
+
+ TLorentzVector fEnergyImpulsionCM_a;
+ TLorentzVector fEnergyImpulsionCM_T;
+ TLorentzVector fEnergyImpulsionCM_1;
+ TLorentzVector fEnergyImpulsionCM_2;
+ TLorentzVector fTotalEnergyImpulsionCM;
+
+ // Impulsion Vector3
+ TVector3 fImpulsionLab_a;
+ TVector3 fImpulsionLab_T;
+ TVector3 fImpulsionLab_1;
+ TVector3 fImpulsionLab_2;
+
+ TVector3 fImpulsionCM_a;
+ TVector3 fImpulsionCM_T;
+ TVector3 fImpulsionCM_1;
+ TVector3 fImpulsionCM_2;
+
+ // CM Energy composante & CM impulsion norme
+ Double_t ECM_a;
+ Double_t ECM_T;
+ Double_t ECM_1;
+ Double_t ECM_2;
+ Double_t pCM_a;
+ Double_t pCM_T;
+ Double_t pCM_1;
+ Double_t pCM_2;
+
+ // Mandelstam variable
+ Double_t s;
+
+ // Center of Mass Kinematic
+ Double_t BetaCM;
+
+ public:
+ //SETTERS
+ void SetBeamEnergy(const double& eBeam) {fBeamEnergy = eBeam;}
+ void SetThetaCM(const double& angle) {fThetaCM = angle;}
+ void SetPhiCM(const double& angle) {fPhiCM = angle;}
+
+ //GETTERS
+ Nucleus* GetNucleusA() {return &fNucleiA;}
+ Nucleus* GetNucleusT() {return &fNucleiT;}
+ Nucleus* GetNucleusB() {return &fNucleiB;}
+ Nucleus* GetNucleus1() {return &fNuclei1;}
+ Nucleus* GetNucleus2() {return &fNuclei2;}
+ bool GetShoot1() const {return fshoot1;}
+ bool GetShoot2() const {return fshoot2;}
+ bool GetShootB() const {return fshootB;}
+ TLorentzVector* GetEnergyImpulsionLab_B() {return &fEnergyImpulsionLab_B;}
+ TGraph* GetTheta2VsTheta1(double AngleStep_CM=1);
+ TGraph* GetPhi2VsPhi1(double AngleStep_CM=1);
+
+
+
+ /*
+ //TO REMOVE AT SOME POINT WHEN CLASS IS ROBUSTLY TESTED
+ private:
+ // R3B Methods and Variables used as a starting point for this class (useful for checks)
+ void CalculateVariablesOld();
+ pair<double,double> Lorentz(double, double, double, double);
+ double function(double, double, double);
+ void KineR3B(double, double, double, double, double);
+ TVector3 Rotations(TVector3,TVector3);
+ double e_clust;
+ double p_clust;
+ double theta_clust;
+ double e_scat;
+ double p_scat;
+ double theta_scat;
+ bool good;
+ double T;
+ */
+
ClassDef(QFS,0)
};
}
diff --git a/NPLib/TrackReconstruction/Tracking.cxx b/NPLib/TrackReconstruction/Tracking.cxx
index 02bda352a39d2814fe2ec96193cfe7af00f61eae..3edd675d5d302d1fd0f55f0fed83fecf46bb7763 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 ee158a5bb2061cc60c6dd5f29512cba9c36ebd08..a198da8a90339dc31c3426fea2771a29cbf19048 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 80c4ae3610492bf49b8a4cf3f7be98e53e182b93..9f80915593327e6e003b4362ad476fc31223f780 100644
--- a/NPSimulation/Core/MaterialManager.cc
+++ b/NPSimulation/Core/MaterialManager.cc
@@ -148,11 +148,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;
}
@@ -358,6 +359,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)
@@ -747,6 +799,15 @@ 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;
@@ -754,6 +815,7 @@ G4Material* MaterialManager::GetMaterialFromLibrary(string Name,
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;
}
@@ -990,6 +1052,15 @@ G4Material* MaterialManager::GetGasFromLibrary(string Name, double Pressure, dou
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);
diff --git a/NPSimulation/Detectors/Dali/Dali.cc b/NPSimulation/Detectors/Dali/Dali.cc
index 84c1eb247c175828562e59f22c37809277ba84ab..c8848e87108abf3733198ef73ff4883ed0a417ad 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 02fb8c6dd394c31c5b3180526c1ca79f131e3276..9b265abcecce41596b228bf9ebd46d77374f0deb 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 121b390ae7d0a0403ab6412eb943044a3d119f0a..a9cc583fb81c778c83dfa7fce8e979390f03bb2c 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 5c5db0349e5335049833b53a41236d7a513a6233..71ca4078b5e99e9a2fd74708da22646b7250eb62 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/NPSimulation/Process/BeamReaction.cc b/NPSimulation/Process/BeamReaction.cc
index 5a94c3020a0e63b004ac278a9c5469f4e5780139..07f42fcf5a43f69cea71093467eedfa6b5a1896a 100644
--- a/NPSimulation/Process/BeamReaction.cc
+++ b/NPSimulation/Process/BeamReaction.cc
@@ -58,19 +58,37 @@ void NPS::BeamReaction::AttachReactionConditions() {
////////////////////////////////////////////////////////////////////////////////
void NPS::BeamReaction::ReadConfiguration() {
NPL::InputParser input(NPOptionManager::getInstance()->GetReactionFile());
- m_Reaction.ReadConfigurationFile(input);
- m_BeamName = NPL::ChangeNameToG4Standard(m_Reaction.GetNucleus1()->GetName());
-
- if (m_Reaction.GetNucleus3()->GetName() != "") {
- m_active = true;
- m_ReactionConditions = new TReactionConditions();
- AttachReactionConditions();
- if (!RootOutput::getInstance()->GetTree()->FindBranch("ReactionConditions"))
- RootOutput::getInstance()->GetTree()->Branch(
- "ReactionConditions", "TReactionConditions", &m_ReactionConditions);
-
- } else {
- m_active = false;
+
+ //vector<NPL::InputBlock*> blocks;
+ //blocks = input.GetAllBlocksWithToken("TwoBodyReaction");
+ //if(blocks.size()>0) m_ReactionType ="TwoBodyReaction";
+ //
+ //blocks = input.GetAllBlocksWithToken("QFSReaction");
+ //if(blocks.size()>0) m_ReactionType ="QFSReaction";
+
+ if(input.GetAllBlocksWithToken("TwoBodyReaction").size()>0) m_ReactionType ="TwoBodyReaction";
+ if(input.GetAllBlocksWithToken("QFSReaction").size()>0) m_ReactionType ="QFSReaction";
+
+
+ if (m_ReactionType=="TwoBodyReaction" ) {
+ m_Reaction.ReadConfigurationFile(input);
+ m_BeamName = NPL::ChangeNameToG4Standard(m_Reaction.GetNucleus1()->GetName());
+ if(m_Reaction.GetNucleus3()->GetName() != ""){
+ m_active = true;
+ m_ReactionConditions = new TReactionConditions();
+ AttachReactionConditions();
+ if (!RootOutput::getInstance()->GetTree()->FindBranch("ReactionConditions"))
+ RootOutput::getInstance()->GetTree()->Branch(
+ "ReactionConditions", "TReactionConditions", &m_ReactionConditions);
+ }
+ } else if (m_ReactionType=="QFSReaction") {
+ m_QFS.ReadConfigurationFile(input);
+ m_BeamName = NPL::ChangeNameToG4Standard(m_QFS.GetNucleusA()->GetName());
+ m_active = true;
+ m_ReactionConditions = new TReactionConditions();
+ }
+ else {
+ m_active = false;
}
}
@@ -90,6 +108,8 @@ NPS::BeamReaction::IsApplicable(const G4ParticleDefinition& particleType) {
////////////////////////////////////////////////////////////////////////////////
G4bool NPS::BeamReaction::ModelTrigger(const G4FastTrack& fastTrack) {
+
+ //cout<< "MODEL TRIG"<<endl;
static bool shoot = false;
static double rand = 0;
const G4Track* PrimaryTrack = fastTrack.GetPrimaryTrack();
@@ -103,9 +123,13 @@ G4bool NPS::BeamReaction::ModelTrigger(const G4FastTrack& fastTrack) {
double out = solid->DistanceToOut(P, -V);
double ratio = in / (out + in);
+ //cout<< "in:"<<in<<std::scientific<<endl;
+ //cout<< "ou:"<<out<<std::scientific<<endl;
+ //cout<< "ratio:"<<ratio<<std::scientific<<endl;
// m_StepSize = m_StepSize/100000.;
if (out == 0) { // first step of current event
+ //cout<< "FIRST"<<endl;
rand = G4RandFlat::shoot();
m_PreviousLength = m_StepSize;
m_PreviousEnergy = PrimaryTrack->GetKineticEnergy();
@@ -113,21 +137,35 @@ G4bool NPS::BeamReaction::ModelTrigger(const G4FastTrack& fastTrack) {
m_ReactionConditions->Clear();
shoot = true;
}
-
- else if (in == 0) { // last step
+ else if ((in-m_StepSize) <= 1E-9) { // last step
+ //cout<< "LAST"<<endl;
return true;
}
+ //cout.precision(17);
+ //cout<< "rand:"<<rand<<std::scientific<<endl;
+
// If the condition is met, the event is generated
if (ratio < rand) {
+
// Reset the static for next event
// shoot = false;
- if (shoot && m_Reaction.IsAllowed(PrimaryTrack->GetKineticEnergy())) {
- shoot = false;
- return true;
-
- } else {
- return false;
+ if(m_ReactionType=="QFSReaction"){
+ if ( shoot && m_QFS.IsAllowed() ) {
+ shoot = false;
+ return true;
+ } else {
+ return false;
+ }
+ }else if(m_ReactionType=="TwoBodyReaction"){
+ if ( shoot && m_Reaction.IsAllowed(PrimaryTrack->GetKineticEnergy()) ) {
+ shoot = false;
+ return true;
+ } else {
+ return false;
+ }
+ }else{
+ return false;
}
}
@@ -135,6 +173,7 @@ G4bool NPS::BeamReaction::ModelTrigger(const G4FastTrack& fastTrack) {
// so it can be used in the next one
if (!PrimaryTrack->GetStep()->IsLastStepInVolume()) {
m_PreviousLength = PrimaryTrack->GetStep()->GetStepLength();
+ //cout<< "PreviousLength="<<m_PreviousLength<<endl;
m_PreviousEnergy = PrimaryTrack->GetKineticEnergy();
}
@@ -143,207 +182,337 @@ G4bool NPS::BeamReaction::ModelTrigger(const G4FastTrack& fastTrack) {
////////////////////////////////////////////////////////////////////////////////
void NPS::BeamReaction::DoIt(const G4FastTrack& fastTrack,
- G4FastStep& fastStep) {
-
- // Get the track info
- const G4Track* PrimaryTrack = fastTrack.GetPrimaryTrack();
- G4ThreeVector pdirection = PrimaryTrack->GetMomentum().unit();
- G4ThreeVector localdir = fastTrack.GetPrimaryTrackLocalDirection();
-
- G4ThreeVector worldPosition = PrimaryTrack->GetPosition();
- G4ThreeVector localPosition = fastTrack.GetPrimaryTrackLocalPosition();
-
- double energy = PrimaryTrack->GetKineticEnergy();
- double time = PrimaryTrack->GetGlobalTime();
-
- // Randomize within the step
- // Assume energy loss is linear within the step
- // Assume no scattering
- double rand = G4RandFlat::shoot();
- double length = rand * (m_PreviousLength);
- energy -= (1 - rand) * (m_PreviousEnergy - energy);
- G4ThreeVector ldir = pdirection;
- ldir *= length;
- localPosition = localPosition - ldir;
- // Set the end of the step conditions
- fastStep.SetPrimaryTrackFinalKineticEnergyAndDirection(0, pdirection);
- fastStep.SetPrimaryTrackFinalPosition(worldPosition);
- fastStep.SetTotalEnergyDeposited(0);
- fastStep.SetPrimaryTrackFinalTime(time);
- fastStep.KillPrimaryTrack();
- fastStep.SetPrimaryTrackPathLength(0.0);
-
- /////////////////////////////////////////////////
- //////Define the kind of particle to shoot////////
- //////////////////////////////////////////////////
-
- // Nucleus 3
- int LightZ = m_Reaction.GetNucleus3()->GetZ();
- int LightA = m_Reaction.GetNucleus3()->GetA();
- static G4IonTable* IonTable
- = G4ParticleTable::GetParticleTable()->GetIonTable();
-
- G4ParticleDefinition* LightName;
-
- if (LightZ == 0 && LightA == 1) // neutron is special case
- {
- LightName = G4Neutron::Definition();
- } else {
- if (m_Reaction.GetUseExInGeant4())
- LightName
- = IonTable->GetIon(LightZ, LightA, m_Reaction.GetExcitation3() * MeV);
- else
- LightName = IonTable->GetIon(LightZ, LightA);
- }
+ G4FastStep& fastStep) {
+
+ // Get the track info
+ const G4Track* PrimaryTrack = fastTrack.GetPrimaryTrack();
+ G4ThreeVector pdirection = PrimaryTrack->GetMomentum().unit();
+ G4ThreeVector localdir = fastTrack.GetPrimaryTrackLocalDirection();
+
+ G4ThreeVector worldPosition = PrimaryTrack->GetPosition();
+ G4ThreeVector localPosition = fastTrack.GetPrimaryTrackLocalPosition();
+
+ double energy = PrimaryTrack->GetKineticEnergy();
+ double time = PrimaryTrack->GetGlobalTime();
+
+ // Randomize within the step
+ // Assume energy loss is linear within the step
+ // Assume no scattering
+ double rand = G4RandFlat::shoot();
+ double length = rand * (m_PreviousLength);
+ energy -= (1 - rand) * (m_PreviousEnergy - energy);
+ G4ThreeVector ldir = pdirection;
+ ldir *= length;
+ localPosition = localPosition - ldir;
+ // Set the end of the step conditions
+ fastStep.SetPrimaryTrackFinalKineticEnergyAndDirection(0, pdirection);
+ fastStep.SetPrimaryTrackFinalPosition(worldPosition);
+ fastStep.SetTotalEnergyDeposited(0);
+ fastStep.SetPrimaryTrackFinalTime(time);
+ fastStep.KillPrimaryTrack();
+ fastStep.SetPrimaryTrackPathLength(0.0);
+
+
+ if (m_ReactionType=="TwoBodyReaction" ) {
+
+ ///////////////////////////////
+ // Two-Body Reaction Case /////
+ ///////////////////////////////
+
+ //////Define the kind of particle to shoot////////
+ // Nucleus 3
+ int LightZ = m_Reaction.GetNucleus3()->GetZ();
+ int LightA = m_Reaction.GetNucleus3()->GetA();
+ static G4IonTable* IonTable
+ = G4ParticleTable::GetParticleTable()->GetIonTable();
+
+ G4ParticleDefinition* LightName;
+
+ if (LightZ == 0 && LightA == 1) // neutron is special case
+ {
+ LightName = G4Neutron::Definition();
+ } else {
+ if (m_Reaction.GetUseExInGeant4())
+ LightName
+ = IonTable->GetIon(LightZ, LightA, m_Reaction.GetExcitation3() * MeV);
+ else
+ LightName = IonTable->GetIon(LightZ, LightA);
+ }
+
+ // Nucleus 4
+ G4int HeavyZ = m_Reaction.GetNucleus4()->GetZ();
+ G4int HeavyA = m_Reaction.GetNucleus4()->GetA();
+
+ // Generate the excitation energy if a distribution is given
+ m_Reaction.ShootRandomExcitationEnergy();
+
+ // Use to clean up the IonTable in case of the Ex changing at every event
+ G4ParticleDefinition* HeavyName;
+
+ if (m_Reaction.GetUseExInGeant4())
+ HeavyName
+ = IonTable->GetIon(HeavyZ, HeavyA, m_Reaction.GetExcitation4() * MeV);
+ else
+ HeavyName = IonTable->GetIon(HeavyZ, HeavyA);
+
+ // Set the Energy of the reaction
+ m_Reaction.SetBeamEnergy(energy);
+
+ double Beam_theta = pdirection.theta();
+ double Beam_phi = pdirection.phi();
+
+ ///////////////////////////
+ ///// Beam Parameters /////
+ ///////////////////////////
+ m_ReactionConditions->SetBeamParticleName(
+ PrimaryTrack->GetParticleDefinition()->GetParticleName());
+
+ m_ReactionConditions->SetBeamReactionEnergy(energy);
+ m_ReactionConditions->SetVertexPositionX(localPosition.x());
+ m_ReactionConditions->SetVertexPositionY(localPosition.y());
+ m_ReactionConditions->SetVertexPositionZ(localPosition.z());
+
+ G4ThreeVector U(1, 0, 0);
+ G4ThreeVector V(0, 1, 0);
+ G4ThreeVector ZZ(0, 0, 1);
+ m_ReactionConditions->SetBeamEmittanceTheta(
+ PrimaryTrack->GetMomentumDirection().theta() / deg);
+ m_ReactionConditions->SetBeamEmittancePhi(
+ PrimaryTrack->GetMomentumDirection().phi() / deg);
+ m_ReactionConditions->SetBeamEmittanceThetaX(
+ PrimaryTrack->GetMomentumDirection().angle(U) / deg);
+ m_ReactionConditions->SetBeamEmittancePhiY(
+ PrimaryTrack->GetMomentumDirection().angle(V) / deg);
+
+ //////////////////////////////////////////////////////////
+ ///// Build rotation matrix to go from the incident //////
+ ///// beam frame to the "world" frame //////
+ //////////////////////////////////////////////////////////
+
+
+ // G4ThreeVector col1(cos(Beam_theta) * cos(Beam_phi),
+ // cos(Beam_theta) * sin(Beam_phi),
+ // -sin(Beam_theta));
+ // G4ThreeVector col2(-sin(Beam_phi),
+ // cos(Beam_phi),
+ // 0);
+ // G4ThreeVector col3(sin(Beam_theta) * cos(Beam_phi),
+ // sin(Beam_theta) * sin(Beam_phi),
+ // cos(Beam_theta));
+ // G4RotationMatrix BeamToWorld(col1, col2, col3);
+
+ /////////////////////////////////////////////////////////////////
+ ///// Angles for emitted particles following Cross Section //////
+ ///// Angles are in the beam frame //////
+ /////////////////////////////////////////////////////////////////
+
+ // Angles
+ // Shoot and Set a Random ThetaCM
+ m_Reaction.ShootRandomThetaCM();
+ double phi = RandFlat::shoot() * 2. * pi;
+
+ //////////////////////////////////////////////////
+ ///// Momentum and angles from kinematics /////
+ ///// Angles are in the beam frame /////
+ //////////////////////////////////////////////////
+ // Variable where to store results
+ double Theta3, Energy3, Theta4, Energy4;
+
+ // Compute Kinematic using previously defined ThetaCM
+ m_Reaction.KineRelativistic(Theta3, Energy3, Theta4, Energy4);
+ // Momentum in beam frame for light particle
+ G4ThreeVector momentum_kine3_beam(sin(Theta3) * cos(phi),
+ sin(Theta3) * sin(phi), cos(Theta3));
+ // Momentum in World frame //to go from the incident beam frame to the "world"
+ // frame
+ G4ThreeVector momentum_kine3_world = momentum_kine3_beam;
+ momentum_kine3_world.rotate(Beam_theta,
+ V); // rotation of Beam_theta on Y axis
+ momentum_kine3_world.rotate(Beam_phi, ZZ); // rotation of Beam_phi on Z axis
+
+ // Momentum in beam frame for heavy particle
+ G4ThreeVector momentum_kine4_beam(sin(Theta4) * cos(phi + pi),
+ sin(Theta4) * sin(phi + pi), cos(Theta4));
+ // Momentum in World frame
+ G4ThreeVector momentum_kine4_world = momentum_kine4_beam;
+ momentum_kine4_world.rotate(Beam_theta,
+ V); // rotation of Beam_theta on Y axis
+ momentum_kine4_world.rotate(Beam_phi, ZZ); // rotation of Beam_phi on Z axis
+
+ // Emitt secondary
+ if (m_Reaction.GetShoot3()) {
+ G4DynamicParticle particle3(LightName, momentum_kine3_world, Energy3);
+ fastStep.CreateSecondaryTrack(particle3, localPosition, time);
+ }
+
+ if (m_Reaction.GetShoot4()) {
+ G4DynamicParticle particle4(HeavyName, momentum_kine4_world, Energy4);
+ fastStep.CreateSecondaryTrack(particle4, localPosition, time);
+ }
+
+ // Reinit for next event
+ m_PreviousEnergy = 0;
+ m_PreviousLength = 0;
+
+ ///////////////////////////////////////
+ ///// Emitted particle Parameters /////
+ ///////////////////////////////////////
+ // Names 3 and 4//
+ m_ReactionConditions->SetParticleName(LightName->GetParticleName());
+ m_ReactionConditions->SetParticleName(HeavyName->GetParticleName());
+ // Angle 3 and 4 //
+ m_ReactionConditions->SetTheta(Theta3 / deg);
+ m_ReactionConditions->SetTheta(Theta4 / deg);
+
+ m_ReactionConditions->SetPhi(phi / deg);
+ if ((phi + pi) / deg > 360)
+ m_ReactionConditions->SetPhi((phi - pi) / deg);
+ else
+ m_ReactionConditions->SetPhi((phi + pi) / deg);
+
+ // Energy 3 and 4 //
+ m_ReactionConditions->SetKineticEnergy(Energy3);
+ m_ReactionConditions->SetKineticEnergy(Energy4);
+ // ThetaCM and Ex//
+ m_ReactionConditions->SetThetaCM(m_Reaction.GetThetaCM() / deg);
+ m_ReactionConditions->SetExcitationEnergy3(m_Reaction.GetExcitation3());
+ m_ReactionConditions->SetExcitationEnergy4(m_Reaction.GetExcitation4());
+ // Momuntum X 3 and 4 //
+ m_ReactionConditions->SetMomentumDirectionX(momentum_kine3_world.x());
+ m_ReactionConditions->SetMomentumDirectionX(momentum_kine4_world.x());
+ // Momuntum Y 3 and 4 //
+ m_ReactionConditions->SetMomentumDirectionY(momentum_kine3_world.y());
+ m_ReactionConditions->SetMomentumDirectionY(momentum_kine4_world.y());
+ // Momuntum Z 3 and 4 //
+ m_ReactionConditions->SetMomentumDirectionZ(momentum_kine3_world.z());
+ m_ReactionConditions->SetMomentumDirectionZ(momentum_kine4_world.z());
+
+ }// end if TwoBodyReaction
+
+ else if (m_ReactionType=="QFSReaction" ) {
+
+ //////Define the kind of particle to shoot////////
+ // A --> T ==> B + (c -> T) => B + 1 + 2
+
+ int Light1_Z = m_QFS.GetNucleus1()->GetZ();
+ int Light1_A = m_QFS.GetNucleus1()->GetA();
+ int Light2_Z = m_QFS.GetNucleus2()->GetZ();
+ int Light2_A = m_QFS.GetNucleus2()->GetA();
+
+ static G4IonTable* IonTable
+ = G4ParticleTable::GetParticleTable()->GetIonTable();
+
+ G4ParticleDefinition* Light1Name;
+ G4ParticleDefinition* Light2Name;
+
+ if (Light1_Z == 0 && Light1_A == 1) // neutron is special case
+ {
+ Light1Name = G4Neutron::Definition();
+ } else {
+ Light1Name = IonTable->GetIon(Light1_Z, Light1_A);
+ }
+
+ if (Light2_Z == 0 && Light2_A == 1) // neutron is special case
+ {
+ Light2Name = G4Neutron::Definition();
+ } else {
+ Light2Name = IonTable->GetIon(Light2_Z, Light2_A);
+ }
+
+ // Nucleus B
+ G4int Heavy_Z = m_QFS.GetNucleusB()->GetZ();
+ G4int Heavy_A = m_QFS.GetNucleusB()->GetA();
+
+ G4ParticleDefinition* HeavyName;
+ HeavyName = IonTable->GetIon(Heavy_Z, Heavy_A);
+
+ // Set the Energy of the reaction
+ m_QFS.SetBeamEnergy(energy);
+
+ double Beam_theta = pdirection.theta();
+ double Beam_phi = pdirection.phi();
+
+
+ /////////////////////////////////////////////////////////////////
+ ///// Angles for emitted particles following Cross Section //////
+ ///// Angles are in the beam frame //////
+ /////////////////////////////////////////////////////////////////
+
+ // Angles
+ // Shoot and Set a Random ThetaCM
+ //m_QFS.ShootRandomThetaCM();
+ //m_QFS.ShootRandomPhiCM();
+ double theta = RandFlat::shoot() * pi;
+ double phi = RandFlat::shoot() * 2. * pi;
+
+ m_QFS.SetThetaCM(theta);
+ m_QFS.SetPhiCM(phi);
+
+ //////////////////////////////////////////////////
+ ///// Momentum and angles from kinematics /////
+ //////////////////////////////////////////////////
+ // Variable where to store results
+ double Theta1, Phi1, Energy1, Theta2, Phi2, Energy2;
+
+ // Compute Kinematic using previously defined ThetaCM
+ m_QFS.KineRelativistic(Theta1, Phi1, Energy1, Theta2, Phi2, Energy2);
+
+ G4ThreeVector U(1, 0, 0);
+ G4ThreeVector V(0, 1, 0);
+ G4ThreeVector ZZ(0, 0, 1);
+
+ // Momentum in beam and world frame for light particle 1
+ G4ThreeVector momentum_kine1_beam(sin(Theta1) * cos(Phi1),
+ sin(Theta1) * sin(Phi1), cos(Theta1));
+ G4ThreeVector momentum_kine1_world = momentum_kine1_beam;
+ momentum_kine1_world.rotate(Beam_theta, V); // rotation of Beam_theta on Y axis
+ momentum_kine1_world.rotate(Beam_phi, ZZ); // rotation of Beam_phi on Z axis
+
+ // Momentum in beam and world frame for light particle 2
+ G4ThreeVector momentum_kine2_beam(sin(Theta2) * cos(Phi2),
+ sin(Theta2) * sin(Phi2), cos(Theta2));
+ G4ThreeVector momentum_kine2_world = momentum_kine2_beam;
+ momentum_kine2_world.rotate(Beam_theta, V); // rotation of Beam_theta on Y axis
+ momentum_kine2_world.rotate(Beam_phi, ZZ); // rotation of Beam_phi on Z axis
+
+ // Momentum in beam and world frame for heavy residual
+ //
+ //G4ThreeVector momentum_kineB_beam(sin(ThetaB) * cos(PhiB + pi),
+ // sin(ThetaB) * sin(PhiB + pi), cos(ThetaB));
+ //G4ThreeVector momentum_kineB_world = momentum_kineB_beam;
+ //momentum_kineB_world.rotate(Beam_theta, V); // rotation of Beam_theta on Y axis
+ //momentum_kineB_world.rotate(Beam_phi, ZZ); // rotation of Beam_phi on Z axis
+
+ TLorentzVector* P_B = m_QFS.GetEnergyImpulsionLab_B();
+
+ G4ThreeVector momentum_kineB_beam( P_B->Px(), P_B->Py(), P_B->Pz() );
+ momentum_kineB_beam = momentum_kineB_beam.unit();
+ double EnergyB = m_QFS.GetEnergyImpulsionLab_B()->Energy();
+ G4ThreeVector momentum_kineB_world = momentum_kineB_beam;
+ momentum_kineB_world.rotate(Beam_theta, V); // rotation of Beam_theta on Y axis
+ momentum_kineB_world.rotate(Beam_phi, ZZ); // rotation of Beam_phi on Z axis
- // Nucleus 4
- G4int HeavyZ = m_Reaction.GetNucleus4()->GetZ();
- G4int HeavyA = m_Reaction.GetNucleus4()->GetA();
-
- // Generate the excitation energy if a distribution is given
- m_Reaction.ShootRandomExcitationEnergy();
-
- // Use to clean up the IonTable in case of the Ex changing at every event
- G4ParticleDefinition* HeavyName;
-
- if (m_Reaction.GetUseExInGeant4())
- HeavyName
- = IonTable->GetIon(HeavyZ, HeavyA, m_Reaction.GetExcitation4() * MeV);
- else
- HeavyName = IonTable->GetIon(HeavyZ, HeavyA);
-
- // Set the Energy of the reaction
- m_Reaction.SetBeamEnergy(energy);
-
- double Beam_theta = pdirection.theta();
- double Beam_phi = pdirection.phi();
-
- ///////////////////////////
- ///// Beam Parameters /////
- ///////////////////////////
- m_ReactionConditions->SetBeamParticleName(
- PrimaryTrack->GetParticleDefinition()->GetParticleName());
-
- m_ReactionConditions->SetBeamReactionEnergy(energy);
- m_ReactionConditions->SetVertexPositionX(localPosition.x());
- m_ReactionConditions->SetVertexPositionY(localPosition.y());
- m_ReactionConditions->SetVertexPositionZ(localPosition.z());
-
- G4ThreeVector U(1, 0, 0);
- G4ThreeVector V(0, 1, 0);
- G4ThreeVector ZZ(0, 0, 1);
- m_ReactionConditions->SetBeamEmittanceTheta(
- PrimaryTrack->GetMomentumDirection().theta() / deg);
- m_ReactionConditions->SetBeamEmittancePhi(
- PrimaryTrack->GetMomentumDirection().phi() / deg);
- m_ReactionConditions->SetBeamEmittanceThetaX(
- PrimaryTrack->GetMomentumDirection().angle(U) / deg);
- m_ReactionConditions->SetBeamEmittancePhiY(
- PrimaryTrack->GetMomentumDirection().angle(V) / deg);
-
- //////////////////////////////////////////////////////////
- ///// Build rotation matrix to go from the incident //////
- ///// beam frame to the "world" frame //////
- //////////////////////////////////////////////////////////
-
- /*
-
- G4ThreeVector col1(cos(Beam_theta) * cos(Beam_phi),
- cos(Beam_theta) * sin(Beam_phi),
- -sin(Beam_theta));
- G4ThreeVector col2(-sin(Beam_phi),
- cos(Beam_phi),
- 0);
- G4ThreeVector col3(sin(Beam_theta) * cos(Beam_phi),
- sin(Beam_theta) * sin(Beam_phi),
- cos(Beam_theta));
- G4RotationMatrix BeamToWorld(col1, col2, col3);
-
- */
-
- /////////////////////////////////////////////////////////////////
- ///// Angles for emitted particles following Cross Section //////
- ///// Angles are in the beam frame //////
- /////////////////////////////////////////////////////////////////
-
- // Angles
- // Shoot and Set a Random ThetaCM
- m_Reaction.ShootRandomThetaCM();
- double phi = RandFlat::shoot() * 2. * pi;
-
- //////////////////////////////////////////////////
- ///// Momentum and angles from kinematics /////
- ///// Angles are in the beam frame /////
- //////////////////////////////////////////////////
- // Variable where to store results
- double Theta3, Energy3, Theta4, Energy4;
-
- // Compute Kinematic using previously defined ThetaCM
- m_Reaction.KineRelativistic(Theta3, Energy3, Theta4, Energy4);
- // Momentum in beam frame for light particle
- G4ThreeVector momentum_kine3_beam(sin(Theta3) * cos(phi),
- sin(Theta3) * sin(phi), cos(Theta3));
- // Momentum in World frame //to go from the incident beam frame to the "world"
- // frame
- G4ThreeVector momentum_kine3_world = /*BeamToWorld */ momentum_kine3_beam;
- momentum_kine3_world.rotate(Beam_theta,
- V); // rotation of Beam_theta on Y axis
- momentum_kine3_world.rotate(Beam_phi, ZZ); // rotation of Beam_phi on Z axis
-
- // Momentum in beam frame for heavy particle
- G4ThreeVector momentum_kine4_beam(sin(Theta4) * cos(phi + pi),
- sin(Theta4) * sin(phi + pi), cos(Theta4));
- // Momentum in World frame
- G4ThreeVector momentum_kine4_world = /*BeamToWorld */ momentum_kine4_beam;
- momentum_kine4_world.rotate(Beam_theta,
- V); // rotation of Beam_theta on Y axis
- momentum_kine4_world.rotate(Beam_phi, ZZ); // rotation of Beam_phi on Z axis
-
- // Emitt secondary
- if (m_Reaction.GetShoot3()) {
- G4DynamicParticle particle3(LightName, momentum_kine3_world, Energy3);
- fastStep.CreateSecondaryTrack(particle3, localPosition, time);
- }
+ // Emitt secondary
+ if (m_QFS.GetShoot1()) {
+ G4DynamicParticle particle1(Light1Name, momentum_kine1_world, Energy1);
+ fastStep.CreateSecondaryTrack(particle1, localPosition, time);
+ }
- if (m_Reaction.GetShoot4()) {
- G4DynamicParticle particle4(HeavyName, momentum_kine4_world, Energy4);
- fastStep.CreateSecondaryTrack(particle4, localPosition, time);
- }
+ if (m_QFS.GetShoot2()) {
+ G4DynamicParticle particle2(Light2Name, momentum_kine2_world, Energy2);
+ fastStep.CreateSecondaryTrack(particle2, localPosition, time);
+ }
+ if (m_QFS.GetShootB()) {
+ G4DynamicParticle particleB(HeavyName, momentum_kineB_world, EnergyB);
+ fastStep.CreateSecondaryTrack(particleB, localPosition, time);
+ }
- // Reinit for next event
- m_PreviousEnergy = 0;
- m_PreviousLength = 0;
- ///////////////////////////////////////
- ///// Emitted particle Parameters /////
- ///////////////////////////////////////
- // Names 3 and 4//
- m_ReactionConditions->SetParticleName(LightName->GetParticleName());
- m_ReactionConditions->SetParticleName(HeavyName->GetParticleName());
- // Angle 3 and 4 //
- m_ReactionConditions->SetTheta(Theta3 / deg);
- m_ReactionConditions->SetTheta(Theta4 / deg);
-
- m_ReactionConditions->SetPhi(phi / deg);
- if ((phi + pi) / deg > 360)
- m_ReactionConditions->SetPhi((phi - pi) / deg);
- else
- m_ReactionConditions->SetPhi((phi + pi) / deg);
-
- // Energy 3 and 4 //
- m_ReactionConditions->SetKineticEnergy(Energy3);
- m_ReactionConditions->SetKineticEnergy(Energy4);
- // ThetaCM and Ex//
- m_ReactionConditions->SetThetaCM(m_Reaction.GetThetaCM() / deg);
- m_ReactionConditions->SetExcitationEnergy3(m_Reaction.GetExcitation3());
- m_ReactionConditions->SetExcitationEnergy4(m_Reaction.GetExcitation4());
- // Momuntum X 3 and 4 //
- m_ReactionConditions->SetMomentumDirectionX(momentum_kine3_world.x());
- m_ReactionConditions->SetMomentumDirectionX(momentum_kine4_world.x());
- // Momuntum Y 3 and 4 //
- m_ReactionConditions->SetMomentumDirectionY(momentum_kine3_world.y());
- m_ReactionConditions->SetMomentumDirectionY(momentum_kine4_world.y());
- // Momuntum Z 3 and 4 //
- m_ReactionConditions->SetMomentumDirectionZ(momentum_kine3_world.z());
- m_ReactionConditions->SetMomentumDirectionZ(momentum_kine4_world.z());
+ // Reinit for next event
+ m_PreviousEnergy = 0;
+ m_PreviousLength = 0;
+
+
+
+ }
}
diff --git a/NPSimulation/Process/BeamReaction.hh b/NPSimulation/Process/BeamReaction.hh
index cb81602f18eb28cde9f4fcbd6af5d146f307e288..b77c424d8885112d5152c7b1f2bf32b37999e16d 100644
--- a/NPSimulation/Process/BeamReaction.hh
+++ b/NPSimulation/Process/BeamReaction.hh
@@ -27,6 +27,7 @@
#include "G4VFastSimulationModel.hh"
#include "PhysicsList.hh"
#include "NPReaction.h"
+#include "NPQFS.h"
#include "TReactionConditions.h"
class G4VPhysicalVolume;
namespace NPS{
@@ -44,7 +45,9 @@ namespace NPS{
private:
NPL::Reaction m_Reaction;
+ NPL::QFS m_QFS;
string m_BeamName;
+ string m_ReactionType;
double m_PreviousEnergy;
double m_PreviousLength;
bool m_active;// is the process active
diff --git a/NPSimulation/Simulation.cc b/NPSimulation/Simulation.cc
index 9c4f542b79178f72af0684da75eed06aa5ff56f3..61dcaca7639d3a387d4f7e71b1e49ab7d4d1b3c8 100644
--- a/NPSimulation/Simulation.cc
+++ b/NPSimulation/Simulation.cc
@@ -115,29 +115,29 @@ int main(int argc, char** argv){
G4VisManager* visManager=NULL;
- if(!OptionManager->GetG4BatchMode()){
+
#ifdef G4UI_USE
+ if(!OptionManager->GetG4BatchMode()){
string Path_Macro = getenv("NPTOOL");
Path_Macro+="/NPSimulation/ressources/macro/";
UImanager->ApplyCommand("/control/execute " +Path_Macro+"verbose.mac");
-#ifdef G4VIS_USE
+ #ifdef G4VIS_USE
UImanager->ApplyCommand("/control/execute " +Path_Macro+"aliases.mac");
visManager = new G4VisExecutive("Quiet");
visManager->Initialize();
UImanager->ApplyCommand("/control/execute " +Path_Macro+"vis.mac");
-#endif
+ #endif
if (ui->IsGUI()){
UImanager->ApplyCommand("/control/execute " +Path_Macro+"gui.mac");
}
-#ifdef __APPLE__
+ #ifdef __APPLE__
string command= "osascript ";
command+= getenv("NPTOOL");
command+="/NPSimulation/ressources/scripts/bringtofront.osa & ";
int res =system(command.c_str());
res =0;
-
-#endif
+ #endif
}
else{// if batch mode do not accumulate any track
UImanager->ApplyCommand("/vis/scene/endOfEventAction accumulate 0");
diff --git a/Projects/Dali/Dali.detector b/Projects/Dali/Dali.detector
index 36f67955030657bfdfbb0f76cf478a7d4f499868..b865ae2073e3e6fbf4cf0a92fc07e1d6e3282655 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