From 37b75c260d09968e2b3091a8600a0581e0fabcb9 Mon Sep 17 00:00:00 2001
From: adrien matta <matta@lpccaen.in2p3.fr>
Date: Tue, 10 Jan 2017 13:41:17 +0100
Subject: [PATCH] * Checking behaviour of DE transport with DE scorer -
The Time of arrival of electron was too long - fixed the issue by
skipping DETransport step going above geometrical safety
---
NPLib/Detectors/Chio/TChio_anData.cxx | 33 ++----
NPLib/Detectors/Chio/TChio_anData.h | 24 ++--
NPLib/Detectors/Chio/TChio_anPhysics.cxx | 4 +-
NPLib/Detectors/Chio/TChio_digData.cxx | 30 +++--
NPLib/Detectors/Chio/TChio_digData.h | 23 ++--
NPLib/Detectors/Chio/TChio_digPhysics.cxx | 11 +-
NPSimulation/Detectors/Chio/Chio.cc | 115 ++++++++++++++++---
NPSimulation/Detectors/Chio/Chio.hh | 8 +-
NPSimulation/Process/G4DETransport.cc | 26 +++--
NPSimulation/Process/G4IonizationWithDE.cc | 23 ++--
NPSimulation/Scorers/DriftElectronScorers.cc | 104 +++++++++++++++--
NPSimulation/Scorers/DriftElectronScorers.hh | 37 +++++-
Projects/BeLise/9Li.reaction | 12 +-
Projects/BeLise/chio.detector | 2 +-
14 files changed, 322 insertions(+), 130 deletions(-)
diff --git a/NPLib/Detectors/Chio/TChio_anData.cxx b/NPLib/Detectors/Chio/TChio_anData.cxx
index 12d3fbb9d..7c38f0d88 100644
--- a/NPLib/Detectors/Chio/TChio_anData.cxx
+++ b/NPLib/Detectors/Chio/TChio_anData.cxx
@@ -27,28 +27,19 @@ using namespace std;
ClassImp(TChio_anData)
-TChio_anData::TChio_anData()
-{
- // Default constructor
-
- // (E)
- fChio_an_Energy.clear();
- fChio_an_Energy_pileup.clear();
-}
+TChio_anData::TChio_anData(){
+ }
-TChio_anData::~TChio_anData()
-{
+TChio_anData::~TChio_anData(){
}
-
-void TChio_anData::Clear()
-{
- // (E)
+void TChio_anData::Clear(){
fChio_an_Energy.clear();
- fChio_an_Energy_pileup.clear();
+ fChio_an_Time.clear();
+ fChio_an_DetectorNbr.clear();
}
@@ -56,14 +47,12 @@ void TChio_anData::Clear()
void TChio_anData::Dump() const
{
cout << "XXXXXXXXXXXXXXXXXXXXXXXX New Event XXXXXXXXXXXXXXXXX" << endl;
-
// Chio_an
// (E)
cout << "Chio_an_MultE = " << fChio_an_Energy.size() << endl;
- for (UShort_t i = 0; i < fChio_an_Energy.size(); i++)
- cout << " Energy: " << fChio_an_Energy[i] << endl;
-
- cout << "Chio_an_MultE pileup = " << fChio_an_Energy_pileup.size() << endl;
- for (UShort_t i = 0; i < fChio_an_Energy_pileup.size(); i++)
- cout << " Energy: " << fChio_an_Energy_pileup[i] << endl;
+ for (UShort_t i = 0; i < fChio_an_Energy.size(); i++){
+ cout << " Detector: " << fChio_an_DetectorNbr[i] ;
+ cout << " Energy: " << fChio_an_Energy[i];
+ cout << " Time: " << fChio_an_Time[i] << endl;
+ }
}
diff --git a/NPLib/Detectors/Chio/TChio_anData.h b/NPLib/Detectors/Chio/TChio_anData.h
index 808e66477..ce0160bba 100644
--- a/NPLib/Detectors/Chio/TChio_anData.h
+++ b/NPLib/Detectors/Chio/TChio_anData.h
@@ -28,30 +28,34 @@ using namespace std;
class TChio_anData : public TObject {
private:
// ADC
- vector<UShort_t> fChio_an_Energy;
- vector<UShort_t> fChio_an_Energy_pileup;
+ vector<double> fChio_an_Energy;
+ vector<double> fChio_an_Time;
+ vector<double> fChio_an_DetectorNbr;
public:
TChio_anData();
virtual ~TChio_anData();
void Clear();
- void Clear(const Option_t*) {};
+ void Clear(const Option_t*) {};
void Dump() const;
///////////////////// GETTERS ////////////////////////
// (E)
- UShort_t GetMultE() {return fChio_an_Energy.size();}
- UShort_t GetEnergy(Int_t i) {return fChio_an_Energy.at(i);}
- UShort_t GetMultE_pileup() {return fChio_an_Energy_pileup.size();}
- UShort_t GetEnergy_pileup(Int_t i) {return fChio_an_Energy_pileup.at(i);}
+ UShort_t GetMult() {return fChio_an_Energy.size();}
+ UShort_t GetEnergy(unsigned int i) {return fChio_an_Energy[i];}
+ UShort_t GetTime(unsigned int i) {return fChio_an_Time[i];}
+ UShort_t GetDetectorNbr(unsigned i) {return fChio_an_DetectorNbr[i];}
///////////////////// SETTERS ////////////////////////
// (E)
- void SetEnergy(UShort_t E) {fChio_an_Energy.push_back(E);}
- void SetEnergy_pileup(UShort_t E) {fChio_an_Energy_pileup.push_back(E);}
+ inline void SetEnergyAndTime(double E, double T, unsigned int Det) {
+ fChio_an_Energy.push_back(E);
+ fChio_an_Time.push_back(T);
+ fChio_an_DetectorNbr.push_back(Det);
+ }
- ClassDef(TChio_anData,1) // Chio_anData structure
+ ClassDef(TChio_anData,2) // Chio_anData structure
};
#endif
diff --git a/NPLib/Detectors/Chio/TChio_anPhysics.cxx b/NPLib/Detectors/Chio/TChio_anPhysics.cxx
index 09edc31e7..4d060a2a7 100644
--- a/NPLib/Detectors/Chio/TChio_anPhysics.cxx
+++ b/NPLib/Detectors/Chio/TChio_anPhysics.cxx
@@ -101,9 +101,9 @@ void TChio_anPhysics::AddParameterToCalibrationManager()
void TChio_anPhysics::InitializeRootInputRaw()
{
TChain* inputChain = RootInput::getInstance()->GetChain() ;
- inputChain->SetBranchStatus("CHIO_AN" , true) ;
+ inputChain->SetBranchStatus("ChioAn" , true) ;
inputChain->SetBranchStatus("fChio_an_*" , true) ;
- inputChain->SetBranchAddress("CHIO_AN" , &EventData) ;
+ inputChain->SetBranchAddress("ChioAn" , &EventData) ;
}
// Activated associated Branches and link it to the private member DetectorPhysics address
diff --git a/NPLib/Detectors/Chio/TChio_digData.cxx b/NPLib/Detectors/Chio/TChio_digData.cxx
index 5bc21ed61..6d6fa2d36 100644
--- a/NPLib/Detectors/Chio/TChio_digData.cxx
+++ b/NPLib/Detectors/Chio/TChio_digData.cxx
@@ -26,43 +26,39 @@ using namespace std;
ClassImp(TChio_digData)
-
+////////////////////////////////////////////////////////////////////////////////
TChio_digData::TChio_digData()
{
- // Default constructor
-
- // (E)
- fChio_dig_Energy.clear();
- //sum
- fChio_dig_Sum = 0;
}
-
+////////////////////////////////////////////////////////////////////////////////
TChio_digData::~TChio_digData()
{
}
-
+////////////////////////////////////////////////////////////////////////////////
void TChio_digData::Clear()
{
- // (E)
fChio_dig_Energy.clear();
- //sum
- fChio_dig_Sum = 0;
-}
+ fChio_dig_Time.clear();
+}
+////////////////////////////////////////////////////////////////////////////////
void TChio_digData::Dump() const
{
cout << "XXXXXXXXXXXXXXXXXXXXXXXX New Event XXXXXXXXXXXXXXXXX" << endl;
// (E)
cout << "Chio_dig_Esize = " << fChio_dig_Energy.size() << endl;
- for (UShort_t i = 0; i < fChio_dig_Energy.size(); i++)
- cout << " Energy: " << fChio_dig_Energy[i] << endl;
- //sum
- cout<<"Chio_dig_sum" << fChio_dig_Sum << endl;
+}
+////////////////////////////////////////////////////////////////////////////////
+TGraph* TChio_digData::GetEnergyAsGraph(){
+ TGraph* res = new TGraph (fChio_dig_Energy.size(),&fChio_dig_Time[0],&fChio_dig_Energy[0]);
+ res->Sort();
+ return res;
+
}
diff --git a/NPLib/Detectors/Chio/TChio_digData.h b/NPLib/Detectors/Chio/TChio_digData.h
index 3c1069acc..fe1d56585 100644
--- a/NPLib/Detectors/Chio/TChio_digData.h
+++ b/NPLib/Detectors/Chio/TChio_digData.h
@@ -3,16 +3,16 @@
#include <vector>
#include "TObject.h"
+#include "TGraph.h"
using namespace std;
class TChio_digData : public TObject {
private:
// ADC
- vector<UShort_t> fChio_dig_Energy;
+ vector<double> fChio_dig_Energy;
+ vector<double> fChio_dig_Time;
- //sum
- Int_t fChio_dig_Sum;
public:
TChio_digData();
@@ -24,19 +24,18 @@ class TChio_digData : public TObject {
///////////////////// GETTERS ////////////////////////
// (E)
- UShort_t GetEsize() {return fChio_dig_Energy.size();}
- UShort_t GetEnergy(Int_t i) {return fChio_dig_Energy.at(i);}
-
- //sum
- Int_t GetSum() {return fChio_dig_Sum;}
+ inline unsigned int GetEsize() {return fChio_dig_Energy.size();}
+ inline vector<double> GetEnergy() {return fChio_dig_Energy;}
+ TGraph* GetEnergyAsGraph();
///////////////////// SETTERS ////////////////////////
- //sum
- void SetSum(UShort_t E) {fChio_dig_Sum += (Int_t)E;}
// (E)
- void SetEnergy(UShort_t E) {fChio_dig_Energy.push_back(E);}
+ void SetEnergy(vector<double>& E) {fChio_dig_Energy=E;}
+ void AddEnergyPoint(double& E,double& T){
+ fChio_dig_Energy.push_back(E);
+ fChio_dig_Time.push_back(T);}
- ClassDef(TChio_digData,1) // Chio_digData structure
+ ClassDef(TChio_digData,2) // Chio_digData structure
};
#endif
diff --git a/NPLib/Detectors/Chio/TChio_digPhysics.cxx b/NPLib/Detectors/Chio/TChio_digPhysics.cxx
index 471d3507a..dc23e116a 100644
--- a/NPLib/Detectors/Chio/TChio_digPhysics.cxx
+++ b/NPLib/Detectors/Chio/TChio_digPhysics.cxx
@@ -99,9 +99,9 @@ void TChio_digPhysics::AddParameterToCalibrationManager()
void TChio_digPhysics::InitializeRootInputRaw()
{
TChain* inputChain = RootInput::getInstance()->GetChain() ;
- inputChain->SetBranchStatus("CHIO_DIG" , true) ;
+ inputChain->SetBranchStatus("ChioDig" , true) ;
inputChain->SetBranchStatus("fChio_dig_*" , true) ;
- inputChain->SetBranchAddress("CHIO_DIG" , &EventData) ;
+ inputChain->SetBranchAddress("ChioDig" , &EventData) ;
}
// Activated associated Branches and link it to the private member DetectorPhysics address
@@ -237,12 +237,7 @@ void TChio_digPhysics::BuildSimplePhysicalEvent()
// this value is defined by RC cuircuit in pre-amplifier.
// ask technitian for R and C values
- // read digitizer data into "Energy"
- for(int ch=0;ch<Nch;ch++){
- Energy.push_back((double)EventData->GetEnergy(ch));
- // Energy[ch] = (int)EventData->GetEnergy(ch);
- // cout << EventData->GetEnergy(ch) << endl;
- }
+ Energy = EventData->GetEnergy();
if(Nch!=0){
int ch = 0;
diff --git a/NPSimulation/Detectors/Chio/Chio.cc b/NPSimulation/Detectors/Chio/Chio.cc
index 3c3de7223..458b56a92 100644
--- a/NPSimulation/Detectors/Chio/Chio.cc
+++ b/NPSimulation/Detectors/Chio/Chio.cc
@@ -1,5 +1,5 @@
/*****************************************************************************
- * Copyright (C) 2009-2106 this file is part of the NPTool Project *
+ * Copyright (C) 2009-2016 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
@@ -61,6 +61,10 @@
// CLHEP header
#include "CLHEP/Random/RandGauss.h"
+// ROOT
+#include "TH1D.h"
+#include "TF1.h"
+
using namespace std;
using namespace CLHEP;
@@ -81,7 +85,9 @@ namespace Chio_NS{
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Chio Specific Method
Chio::Chio(){
- m_Event = new TChio_anData() ;
+ m_Event_an = new TChio_anData() ;
+ m_Event_dig = new TChio_digData() ;
+
m_ChioScorer = 0;
m_SquareDetector = 0;
m_CylindricalDetector = 0;
@@ -147,10 +153,10 @@ G4LogicalVolume* Chio::BuildDetector(){
G4MaterialPropertiesTable* MPT = new G4MaterialPropertiesTable();
MPT->AddConstProperty("DE_PAIRENERGY",30*eV);
- MPT->AddConstProperty("DE_YIELD",1e-1);
+ MPT->AddConstProperty("DE_YIELD",1e-2);
// MPT->AddConstProperty("DE_AMPLIFICATION",1e4);
MPT->AddConstProperty("DE_ABSLENGTH",1*pc);
- MPT->AddConstProperty("DE_DRIFTSPEED",8e-3*mm/ns);
+ MPT->AddConstProperty("DE_DRIFTSPEED",13*cm/us);
MPT->AddConstProperty("DE_TRANSVERSALSPREAD",6e-5*mm2/ns);
MPT->AddConstProperty("DE_LONGITUDINALSPREAD",4e-5*mm2/ns);
@@ -175,11 +181,11 @@ G4LogicalVolume* Chio::BuildDetector(){
logicGas,
"ChioGas",m_SquareDetector,false,0);
- new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(2.5*cm,0,0)),
+ new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(-2.5*cm,0,0)),
logicGrid,
"ChioGrid",logicGas,false,0);
- new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(3*cm,0,0)),
+ new G4PVPlacement(G4Transform3D(*Rot,G4ThreeVector(3*cm-0.5*1*um,0,0)),
logicCathode,
"ChioCathode",logicGas,false,0);
@@ -298,41 +304,108 @@ void Chio::ConstructDetector(G4LogicalVolume* world){
void Chio::InitializeRootOutput(){
RootOutput *pAnalysis = RootOutput::getInstance();
TTree *pTree = pAnalysis->GetTree();
- if(!pTree->FindBranch("Chio")){
- pTree->Branch("Chio", "TChioData", &m_Event) ;
+ if(!pTree->FindBranch("ChioAn")){
+ pTree->Branch("ChioAn", "TChio_anData", &m_Event_an) ;
+ }
+ pTree->SetBranchAddress("ChioAn", &m_Event_an) ;
+
+ ///////////////
+ if(!pTree->FindBranch("ChioDig")){
+ pTree->Branch("ChioDig", "TChio_digData", &m_Event_dig) ;
}
- pTree->SetBranchAddress("Chio", &m_Event) ;
+ pTree->SetBranchAddress("ChioDig", &m_Event_dig) ;
+
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Read sensitive part and fill the Root tree.
// Called at in the EventAction::EndOfEventAvtion
void Chio::ReadSensitive(const G4Event* event){
- m_Event->Clear();
+ m_Event_an->Clear();
///////////
- // Cathoderimeter scorer
+ // Cathode analogic scorer
NPS::HitsMap<G4double*>* CathodeHitMap;
std::map<G4int, G4double**>::iterator Cathode_itr;
- G4int CathodeCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("ChioScorer/Cathode");
+ G4int CathodeCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("ChioScorer/Cathode_an");
CathodeHitMap = (NPS::HitsMap<G4double*>*)(event->GetHCofThisEvent()->GetHC(CathodeCollectionID));
// Loop on the Cathode map
for (Cathode_itr = CathodeHitMap->GetMap()->begin() ; Cathode_itr != CathodeHitMap->GetMap()->end() ; Cathode_itr++){
-
G4double* Info = *(Cathode_itr->second);
double Count= Info[0];
if(Count>Chio_NS::ChargeThreshold-1){
double Time = RandGauss::shoot(Info[1],Chio_NS::ResoTime);
- int DetectorNbr = (int) Info[2];
- //cout << Count << " " << Time/ns << endl;
- //m_Event->SetEnergy(DetectorNbr,Energy);
- //m_Event->SetTime(DetectorNbr,Time);
+ m_Event_an->SetEnergyAndTime(Count,Time,Info[2]);
}
}
// clear map for next event
CathodeHitMap->clear();
+
+ m_Event_dig->Clear();
+
+ ///////////
+ CathodeCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("ChioScorer/Cathode_dig");
+ CathodeHitMap = (NPS::HitsMap<G4double*>*)(event->GetHCofThisEvent()->GetHC(CathodeCollectionID));
+
+ // Loop on the Cathode map
+ TH1D* h = new TH1D("h","h",25000,0,25000);
+ for (Cathode_itr = CathodeHitMap->GetMap()->begin() ; Cathode_itr != CathodeHitMap->GetMap()->end() ; Cathode_itr++){
+ G4double* Info = *(Cathode_itr->second);
+ if(Info[0]){
+ h->Fill(Info[1],Info[0]);
+ }
+ }
+
+ vector<double> E,T;
+ for(int i = 0 ; i < h->GetNbinsX() ; i++){
+ double count = h->GetBinContent(i);
+ double time = h->GetBinCenter(i);
+ if(count)
+ // m_Event_dig->AddEnergyPoint(count,time);
+ E.push_back(count);
+ T.push_back(time+5000);
+ }
+
+ SimulateDigitizer(E,T,275*ns,1.40*us,0,12500,25);
+
+ delete h;
+ // clear map for next event
+ CathodeHitMap->clear();
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void Chio::SimulateDigitizer(vector<double> E, vector<double> T,double riseTime, double fallTime,double start,double stop, double step){
+
+ string formula = "";
+ string Es,Ts,var,cond1,cond2;
+ string fall=std::to_string(fallTime);
+ string rise=std::to_string(riseTime);
+
+ for(unsigned int i = 0 ; i < E.size() ; i++){
+ if(E[i]!=0 && T[i]!=0){
+ Es = std::to_string(E[i]);
+ Ts = std::to_string(T[i]);
+ cond1 = ")*(x<"+Ts+"&& x>"+Ts+"-"+rise+ ")+";
+ cond2 = ")*(x>"+Ts+")+";
+ var = "(x-"+Ts+")";
+ formula += Es+"*-1*exp("+var+"/"+rise+cond1
+ +Es+"*-1*exp(-"+var+"/"+fall+cond2;
+ }
+ }
+ formula+="0";
+ TF1* f = new TF1("f",formula.c_str(),start,stop);
+ unsigned int size = (stop-start)/step;
+ for(unsigned int i = 0 ; i < size ; i++){
+ double time = start+i*step;
+ double energy = f->Eval(time);
+
+ m_Event_dig->AddEnergyPoint(energy,time);
+ }
+
+ delete f;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -346,9 +419,13 @@ void Chio::InitializeScorers() {
return ;
// Otherwise the scorer is initialised
- G4VPrimitiveScorer* Cathode= new DRIFTELECTRONSCORERS::PS_DriftElectron("Cathode",0) ;
+ G4VPrimitiveScorer* Cathode_an= new DRIFTELECTRONSCORERS::PS_DECathode("Cathode_an",0) ;
+ G4VPrimitiveScorer* Cathode_dig= new DRIFTELECTRONSCORERS::PS_DEDigitizer("Cathode_dig",0) ;
+
//and register it to the multifunctionnal detector
- m_ChioScorer->RegisterPrimitive(Cathode);
+ m_ChioScorer->RegisterPrimitive(Cathode_an);
+ m_ChioScorer->RegisterPrimitive(Cathode_dig);
+
G4SDManager::GetSDMpointer()->AddNewDetector(m_ChioScorer) ;
}
diff --git a/NPSimulation/Detectors/Chio/Chio.hh b/NPSimulation/Detectors/Chio/Chio.hh
index 65ab3a9cf..76f82cc41 100644
--- a/NPSimulation/Detectors/Chio/Chio.hh
+++ b/NPSimulation/Detectors/Chio/Chio.hh
@@ -1,7 +1,7 @@
#ifndef Chio_h
#define Chio_h 1
/*****************************************************************************
- * Copyright (C) 2009-XYEARX this file is part of the NPTool Project *
+ * Copyright (C) 2009-2017 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
@@ -35,6 +35,8 @@ using namespace std;
// NPTool header
#include "NPSVDetector.hh"
#include "TChio_anData.h"
+#include "TChio_digData.h"
+
#include "NPInputParser.h"
class Chio : public NPS::VDetector{
@@ -84,6 +86,7 @@ class Chio : public NPS::VDetector{
public: // Scorer
// Initialize all Scorer used by the MUST2Array
void InitializeScorers() ;
+ void SimulateDigitizer(vector<double> E, vector<double> T,double riseTime, double fallTime,double start,double stop,double step);
// Associated Scorer
G4MultiFunctionalDetector* m_ChioScorer ;
@@ -91,7 +94,8 @@ class Chio : public NPS::VDetector{
///////////Event class to store Data////////////////
////////////////////////////////////////////////////
private:
- TChio_anData* m_Event ;
+ TChio_anData* m_Event_an ;
+ TChio_digData* m_Event_dig ;
////////////////////////////////////////////////////
///////////////Private intern Data//////////////////
diff --git a/NPSimulation/Process/G4DETransport.cc b/NPSimulation/Process/G4DETransport.cc
index f8d5db264..57cc69460 100644
--- a/NPSimulation/Process/G4DETransport.cc
+++ b/NPSimulation/Process/G4DETransport.cc
@@ -117,7 +117,7 @@ G4DETransport::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
if(step_length<100*micrometer){
return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
}
-
+
// Get the material table
const G4Material* aMaterial = aTrack.GetMaterial();
@@ -150,35 +150,41 @@ G4DETransport::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
// Normalised the drift direction
driftDir = driftDir.unit();
G4double v_drift = aMaterialPropertiesTable->GetConstProperty("DE_DRIFTSPEED");
- // Should be equal to delta length
- G4double step = step_length/v_drift;
-
+
G4double sigmaTrans = sqrt(2*step_length*aMaterialPropertiesTable->GetConstProperty("DE_TRANSVERSALSPREAD")/v_drift);
G4double sigmaLong = sqrt(2*step_length*aMaterialPropertiesTable->GetConstProperty("DE_LONGITUDINALSPREAD")/v_drift);
G4double d_long = G4RandGauss::shoot(0,sigmaLong);
G4double d_trans = G4RandGauss::shoot(0,sigmaTrans);
-
+ G4double d_drift = step_length+d_long;
+
G4ThreeVector trans(G4RandGauss::shoot(0,d_trans),0,0);
trans.rotateY(twopi*G4UniformRand());
- G4ThreeVector d_Pos = (step_length+d_long)*driftDir+trans;
+ G4ThreeVector d_Pos = (d_drift)*driftDir+trans;
- // New particle Position with matching time
+ // Should be equal to delta length
+ G4double step = (d_drift)/v_drift;
+
+ // New particle Position with matching time
G4ThreeVector pos = x0 + d_Pos;
G4double time = t0 + step;
+
// Garanty that the electron does not jump outside the current volume
G4double safety = m_SafetyHelper->ComputeSafety(x0,d_Pos.mag());
+ // If the distance travelled if above safety, the step is not taken
if(d_Pos.mag()>safety){
- pos = x0+safety*d_Pos.unit();
+ return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
}
aParticleChange.ProposeMomentumDirection(d_Pos.unit());
aParticleChange.ProposeEnergy(Energy);
aParticleChange.ProposePosition(pos);
aParticleChange.ProposeGlobalTime(time);
- aParticleChange.ProposeVelocity(v_drift/c_light);
+ aParticleChange.ProposeLocalTime(time);
+ aParticleChange.ProposeVelocity(v_drift/c_light);
+ m_SafetyHelper->ReLocateWithinVolume(pos);
return &aParticleChange;
}
@@ -193,5 +199,5 @@ G4double G4DETransport::GetMeanFreePath(const G4Track& ,
// Typical distance after which the electron position should be reevaluated
// to take into account the diffusivity of the charge cloud inside the
// medium
- return 1*mm;
+ return 5*mm;
}
diff --git a/NPSimulation/Process/G4IonizationWithDE.cc b/NPSimulation/Process/G4IonizationWithDE.cc
index b2aa522a4..b4a2abe6c 100644
--- a/NPSimulation/Process/G4IonizationWithDE.cc
+++ b/NPSimulation/Process/G4IonizationWithDE.cc
@@ -49,7 +49,8 @@
#include "G4SystemOfUnits.hh"
#include "G4ParticleTypes.hh"
#include "G4EmProcessSubType.hh"
-
+#include "G4ElectroMagneticField.hh"
+#include "G4FieldManager.hh"
#include "G4IonizationWithDE.hh"
#include <iostream>
using namespace std;
@@ -183,12 +184,20 @@ G4IonizationWithDE::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
// Create the secondary tracks
for(G4int i = 0 ; i < number_electron ; i++){
- // Random direction at creation
- G4double cost = 1-2*G4UniformRand();
- G4double theta = acos(cost);
- G4double phi = twopi*G4UniformRand();
- G4ThreeVector p;
- p.setRThetaPhi(1,theta,phi);
+ // Electron follow the field direction
+ // The field direction is taken from the field manager
+ G4double* fieldArr = new G4double[6];
+ G4double Point[4]={x0.x(),x0.y(),x0.z(),t0};
+ G4FieldManager* fMng = pPreStepPoint->GetTouchableHandle()->GetVolume()->GetLogicalVolume()->
+ GetFieldManager();
+
+ G4ElectroMagneticField* field = (G4ElectroMagneticField*)fMng->GetDetectorField();
+ field->GetFieldValue(Point,fieldArr) ;
+
+ // Electron move opposite to the field direction, hance the minus sign
+ G4ThreeVector p(-fieldArr[3],-fieldArr[4],-fieldArr[5]);
+ // Normalised the drift direction
+ p = p.unit();
// Random Position along the step with matching time
G4double rand = G4UniformRand();
diff --git a/NPSimulation/Scorers/DriftElectronScorers.cc b/NPSimulation/Scorers/DriftElectronScorers.cc
index ef4847db2..4c0548f42 100644
--- a/NPSimulation/Scorers/DriftElectronScorers.cc
+++ b/NPSimulation/Scorers/DriftElectronScorers.cc
@@ -24,19 +24,18 @@
using namespace DRIFTELECTRONSCORERS ;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-PS_DriftElectron::PS_DriftElectron(G4String name,G4int depth)
+PS_DECathode::PS_DECathode(G4String name,G4int depth)
:G4VPrimitiveScorer(name, depth),HCID(-1){
m_Index = -1 ;
m_Level = depth;
- m_NumberOfDriftElectron = 0;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-PS_DriftElectron::~PS_DriftElectron(){
+PS_DECathode::~PS_DECathode(){
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-G4bool PS_DriftElectron::ProcessHits(G4Step* aStep, G4TouchableHistory*){
+G4bool PS_DECathode::ProcessHits(G4Step* aStep, G4TouchableHistory*){
// contain Energy Time, DetNbr, StripFront and StripBack
G4double* Infos = new G4double[9];
@@ -67,6 +66,7 @@ G4bool PS_DriftElectron::ProcessHits(G4Step* aStep, G4TouchableHistory*){
if(it!=EvtMap->GetMap()->end()){
G4double* dummy = *(it->second);
Infos[0]+=dummy[0];
+ Infos[1]=dummy[1];
}
EvtMap->set(m_Index, Infos);
@@ -74,7 +74,7 @@ G4bool PS_DriftElectron::ProcessHits(G4Step* aStep, G4TouchableHistory*){
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void PS_DriftElectron::Initialize(G4HCofThisEvent* HCE){
+void PS_DECathode::Initialize(G4HCofThisEvent* HCE){
EvtMap = new NPS::HitsMap<G4double*>(GetMultiFunctionalDetector()->GetName(), GetName());
if (HCID < 0) {
HCID = GetCollectionID(0);
@@ -83,11 +83,11 @@ void PS_DriftElectron::Initialize(G4HCofThisEvent* HCE){
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void PS_DriftElectron::EndOfEvent(G4HCofThisEvent*){
+void PS_DECathode::EndOfEvent(G4HCofThisEvent*){
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void PS_DriftElectron::clear(){
+void PS_DECathode::clear(){
std::map<G4int, G4double**>::iterator MapIterator;
for (MapIterator = EvtMap->GetMap()->begin() ; MapIterator != EvtMap->GetMap()->end() ; MapIterator++){
delete *(MapIterator->second);
@@ -97,12 +97,98 @@ void PS_DriftElectron::clear(){
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void PS_DriftElectron::DrawAll(){
+void PS_DECathode::DrawAll(){
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void PS_DriftElectron::PrintAll(){
+void PS_DECathode::PrintAll(){
+ G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl ;
+ G4cout << " PrimitiveScorer " << GetName() << G4endl ;
+ G4cout << " Number of entries " << EvtMap->entries() << G4endl ;
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+PS_DEDigitizer::PS_DEDigitizer(G4String name,G4int depth)
+:G4VPrimitiveScorer(name, depth),HCID(-1){
+ m_Index = -1 ;
+ m_Level = depth;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+PS_DEDigitizer::~PS_DEDigitizer(){
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+G4bool PS_DEDigitizer::ProcessHits(G4Step* aStep, G4TouchableHistory*){
+
+ // contain Energy Time, DetNbr, StripFront and StripBack
+ G4double* Infos = new G4double[9];
+ Infos[0] = 0;
+ Infos[1] = aStep->GetPreStepPoint()->GetGlobalTime();
+
+ m_DetectorNumber = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(m_Level);
+ m_Position = aStep->GetPreStepPoint()->GetPosition();
+
+ // Interaction coordinates (used to fill the InteractionCoordinates branch)
+ Infos[2] = m_Position.x();
+ Infos[3] = m_Position.y();
+ Infos[4] = m_Position.z();
+ Infos[5] = m_Position.theta();
+ Infos[6] = m_Position.phi();
+ Infos[7] = m_DetectorNumber;
+
+ m_Index = aStep->GetTrack()->GetTrackID() + m_DetectorNumber*1e6 ;
+ G4String PID = aStep->GetTrack()->GetDefinition()->GetParticleName();
+
+ if(PID=="driftelectron"){
+ Infos[0] = 1;
+ }
+
+ // Check if the particle has interact before, if yes, add up the number of electron.
+ map<G4int, G4double**>::iterator it;
+ it= EvtMap->GetMap()->find(m_Index);
+ if(it!=EvtMap->GetMap()->end()){
+ G4double* dummy = *(it->second);
+ Infos[0]+=dummy[0];
+ Infos[1]=dummy[1];
+ }
+
+ EvtMap->set(m_Index, Infos);
+ return TRUE;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PS_DEDigitizer::Initialize(G4HCofThisEvent* HCE){
+ EvtMap = new NPS::HitsMap<G4double*>(GetMultiFunctionalDetector()->GetName(), GetName());
+ if (HCID < 0) {
+ HCID = GetCollectionID(0);
+ }
+ HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PS_DEDigitizer::EndOfEvent(G4HCofThisEvent*){
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PS_DEDigitizer::clear(){
+ std::map<G4int, G4double**>::iterator MapIterator;
+ for (MapIterator = EvtMap->GetMap()->begin() ; MapIterator != EvtMap->GetMap()->end() ; MapIterator++){
+ delete *(MapIterator->second);
+ }
+
+ EvtMap->clear();
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PS_DEDigitizer::DrawAll(){
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void PS_DEDigitizer::PrintAll(){
G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl ;
G4cout << " PrimitiveScorer " << GetName() << G4endl ;
G4cout << " Number of entries " << EvtMap->entries() << G4endl ;
diff --git a/NPSimulation/Scorers/DriftElectronScorers.hh b/NPSimulation/Scorers/DriftElectronScorers.hh
index 50f58d88c..6b49c8a67 100644
--- a/NPSimulation/Scorers/DriftElectronScorers.hh
+++ b/NPSimulation/Scorers/DriftElectronScorers.hh
@@ -1,4 +1,4 @@
-#ifndef DriftElectronScorers_h
+#ifndef DECathodeScorers_h
#define SiliconScorers_h 1
/*****************************************************************************
* Copyright (C) 2009-2016 this file is part of the NPTool Project *
@@ -33,11 +33,39 @@ using namespace CLHEP;
namespace DRIFTELECTRONSCORERS {
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
- class PS_DriftElectron : public G4VPrimitiveScorer{
+ class PS_DECathode : public G4VPrimitiveScorer{
public: // with description
- PS_DriftElectron(G4String name,G4int depth=0);
- ~PS_DriftElectron();
+ PS_DECathode(G4String name,G4int depth=0);
+ ~PS_DECathode();
+
+ protected: // with description
+ G4bool ProcessHits(G4Step*, G4TouchableHistory*);
+
+ public:
+ void Initialize(G4HCofThisEvent*);
+ void EndOfEvent(G4HCofThisEvent*);
+ void clear();
+ void DrawAll();
+ void PrintAll();
+
+private: // inherited from G4VPrimitiveScorer
+ G4int HCID;
+ NPS::HitsMap<G4double*>* EvtMap;
+
+ private: // Needed for intermediate calculation (avoid multiple instantiation in Processing Hit)
+ G4int m_Index;
+ G4int m_Level;
+ G4int m_DetectorNumber;
+ G4ThreeVector m_Position;
+
+ };
+////////////////////////////////////////////////////////////////////////////////
+ class PS_DEDigitizer : public G4VPrimitiveScorer{
+
+ public: // with description
+ PS_DEDigitizer(G4String name,G4int depth=0);
+ ~PS_DEDigitizer();
protected: // with description
G4bool ProcessHits(G4Step*, G4TouchableHistory*);
@@ -58,7 +86,6 @@ private: // inherited from G4VPrimitiveScorer
G4int m_Level;
G4int m_DetectorNumber;
G4ThreeVector m_Position;
- G4double m_NumberOfDriftElectron;
};
diff --git a/Projects/BeLise/9Li.reaction b/Projects/BeLise/9Li.reaction
index 02db44e32..8e29608a7 100644
--- a/Projects/BeLise/9Li.reaction
+++ b/Projects/BeLise/9Li.reaction
@@ -7,12 +7,12 @@ Beam
SigmaEnergy= 20 MeV
SigmaThetaX= 0.6138 deg
SigmaPhiY= 0.3812 deg
- SigmaX= 5.1 mm
- SigmaY= 8.5 mm
+ SigmaX= 6.216 mm
+ SigmaY= 6.109 mm
MeanThetaX= 0 deg
MeanPhiY= 0 deg
- MeanX= 2.5 mm
- MeanY= -0.5 mm
+ MeanX= 0 mm
+ MeanY= 0 mm
%EnergyProfilePath=
%XThetaXProfilePath=
%YPhiYProfilePath=
@@ -24,8 +24,8 @@ TwoBodyReaction
Light= 3He
Heavy= 9Li
ExcitationEnergyLight= 0.0 MeV
- ExcitationEnergyHeavy= 0 MeV
- CrossSectionPath= 11Li(d,3He)10He.txt CS10He
+ ExcitationEnergyHeavy= 0.0 MeV
+ CrossSectionPath= 11Li(d,3He)10He.txt CS2
ShootLight= 1
ShootHeavy= 1
diff --git a/Projects/BeLise/chio.detector b/Projects/BeLise/chio.detector
index c7cfe4ce9..c38fa688c 100644
--- a/Projects/BeLise/chio.detector
+++ b/Projects/BeLise/chio.detector
@@ -9,7 +9,7 @@ Target
Z= 0 mm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Chio
- POS= 0 0 350 mm
+ POS= 0 0 400 mm
Shape= Square
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
--
GitLab