diff --git a/VeikonKone.dox b/VeikonKone.dox
index 00957094131304ba0bd8087b1eb1e3ce3d4579a9..6af4a560a0dbddd1861f60537903b95487fd7018 100644
--- a/VeikonKone.dox
+++ b/VeikonKone.dox
@@ -753,7 +753,7 @@ WARN_LOGFILE =
# spaces.
# Note: If this tag is empty the current directory is searched.
-INPUT = /home/oliskir/Desktop/f20/sim/VeikonKone/include/
+INPUT = /home/oliskir/Desktop/F20/sim/VeikonKone/include/
# This tag can be used to specify the character encoding of the source files
# that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses
diff --git a/analysis/VKAdd b/analysis/VKAdd
index 0ed82ea52e580552e2df6a8a0edc3567e9601bb2..90a79678fd6dfa5cdf8c813e5e4fdc92e932a745 100755
Binary files a/analysis/VKAdd and b/analysis/VKAdd differ
diff --git a/analysis/VKAnalysis.dox b/analysis/VKAnalysis.dox
index 276ac49501b062d0a98692434cca60d9cebf3bae..9e71fd97a97a11cbe37e83f0e2e0458eb63c7a2c 100644
--- a/analysis/VKAnalysis.dox
+++ b/analysis/VKAnalysis.dox
@@ -753,7 +753,7 @@ WARN_LOGFILE =
# spaces.
# Note: If this tag is empty the current directory is searched.
-INPUT = /home/oliskir/Desktop/f20/sim/VeikonKone/analysis/include/
+INPUT = /home/oliskir/Desktop/F20/sim/VeikonKone/analysis/include/
# This tag can be used to specify the character encoding of the source files
# that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses
diff --git a/analysis/include/VKBaseData.hh b/analysis/include/VKBaseData.hh
index 187e8bdf2314730a92d817f19fc1244d80b0a188..5cfc95dd535162587d638661553548fa10839e60 100644
--- a/analysis/include/VKBaseData.hh
+++ b/analysis/include/VKBaseData.hh
@@ -14,6 +14,10 @@ struct VKEvent {
Float_t EV; // Energy Veto
UInt_t PS; // Photons Signal
UInt_t PV; // Photons Veto
+
+ Float_t EF2; // Energy Front v2.0
+ Float_t ES2; // Energy Signal v2.0
+ Float_t EV2; // Energy Veto v2.0
};
@@ -50,6 +54,9 @@ class VKBaseData {
Float_t _EVout;
UInt_t _PSout;
UInt_t _PVout;
+ Float_t _EF2out;
+ Float_t _ES2out;
+ Float_t _EV2out;
};
#endif
diff --git a/analysis/include/VKSimData.hh b/analysis/include/VKSimData.hh
index a809fc758e628e820f7891de5f6a551c39a7b524..fabe0d617d92d935f8e473dbe4b716db774d8c38 100644
--- a/analysis/include/VKSimData.hh
+++ b/analysis/include/VKSimData.hh
@@ -31,6 +31,9 @@ class VKSimData : public VKBaseData {
Float_t _EVin;
UInt_t _PSin;
UInt_t _PVin;
+ Float_t _EF2in;
+ Float_t _ES2in;
+ Float_t _EV2in;
};
#endif //VKSIMDATA_H
diff --git a/analysis/src/VKBaseData.cpp b/analysis/src/VKBaseData.cpp
index dfdc710f5bedf3d2d5fb56184ce09262b6412d92..b150c15eb17c0ea81f42a4afe46f2ec42f8c0961 100644
--- a/analysis/src/VKBaseData.cpp
+++ b/analysis/src/VKBaseData.cpp
@@ -56,10 +56,16 @@ void VKBaseData::CreateOutput(string filename)
_EVout = ev.EV;
_PSout = ev.PS;
_PVout = ev.PV;
+ _EF2out = ev.EF2;
+ _ES2out = ev.ES2;
+ _EV2out = ev.EV2;
// sample experimental resolution
SampleResolution(_ESout);
SampleResolution(_EVout);
+ SampleResolution(_EF2out);
+ SampleResolution(_ES2out);
+ SampleResolution(_EV2out);
// fill tree
fOutputTree -> Fill();
@@ -77,10 +83,13 @@ void VKBaseData::CreateOutput(string filename)
void VKBaseData::CreateOutputBranches(TTree * t)
{
- t -> Branch("EDepSignal", &_ESout, "_ESout/F");
- t -> Branch("EDepVeto" , &_EVout, "_EVout/F");
- t -> Branch("PhotonsSignal", &_PSout, "_PSout/i");
- t -> Branch("PhotonsVeto" , &_PVout , "_PVout/i");
+ t->Branch("EDepSignal", &_ESout, "_ESout/F");
+ t->Branch("EDepVeto" , &_EVout, "_EVout/F");
+ t->Branch("PhotonsSignal", &_PSout, "_PSout/i");
+ t->Branch("PhotonsVeto" , &_PVout , "_PVout/i");
+ t->Branch("EDepV2Front", &_EF2out, "_EF2out/F");
+ t->Branch("EDepV2Signal", &_ES2out, "_ES2out/F");
+ t->Branch("EDepV2Veto" , &_EV2out, "_EV2out/F");
}
@@ -88,11 +97,10 @@ void VKBaseData::SampleResolution(Float_t& e)
{
if (!fAddRes) return;
-
double fwhm = sqrt(pow(fP0,2) + pow(fP1,2) * e);
double sigma = fwhm / 2.35;
-
- e += fRNG -> Gaus(0, sigma);
+
+ if (e > 0) e += fRNG -> Gaus(0, sigma);
if (e<=0) e = 0;
}
diff --git a/analysis/src/VKPileupData.cpp b/analysis/src/VKPileupData.cpp
index 3d8dcb4497ea5ac8a51f628f7d64e6cfc44de9a5..6dab791ce01349f4c3593560a94af147074fdb34 100644
--- a/analysis/src/VKPileupData.cpp
+++ b/analysis/src/VKPileupData.cpp
@@ -35,5 +35,8 @@ const VKEvent& VKPileupData::GetEvent(size_t id)
fEvent.EV = ev1.EV + ev2.EV;
fEvent.PS = ev1.PS + ev2.PS;
fEvent.PV = ev1.PV + ev2.PV;
+ fEvent.EF2 = ev1.EF2 + ev2.EF2;
+ fEvent.ES2 = ev1.ES2 + ev2.ES2;
+ fEvent.EV2 = ev1.EV2 + ev2.EV2;
return fEvent;
}
diff --git a/analysis/src/VKSimData.cpp b/analysis/src/VKSimData.cpp
index 260fb153add0b5e9127644971ab5548ad3c2d940..6e1e3d26a14cfe1ded6e369d3b917e609a4ceee1 100644
--- a/analysis/src/VKSimData.cpp
+++ b/analysis/src/VKSimData.cpp
@@ -29,6 +29,14 @@ VKSimData::VKSimData(string filename)
// area
fArea =
((TParameter<double>*)(fInputTree->GetUserInfo()->FindObject("SOURCE_AREA/cm2")))->GetVal();
+
+ _ESin = 0;
+ _EVin = 0;
+ _PSin = 0;
+ _PVin = 0;
+ _EF2in = 0;
+ _ES2in = 0;
+ _EV2in = 0;
}
@@ -45,10 +53,26 @@ VKSimData::~VKSimData()
void VKSimData::SetInputBranches(TTree * t)
{
- t -> SetBranchAddress("EDepSignal", &_ESin);
- t -> SetBranchAddress("EDepVeto" , &_EVin);
- t -> SetBranchAddress("PhotonsSignal", &_PSin);
- t -> SetBranchAddress("PhotonsVeto" , &_PVin);
+ if (t->GetBranch("EDepSignal"))
+ t->SetBranchAddress("EDepSignal", &_ESin);
+
+ if (t->GetBranch("EDepVeto"))
+ t->SetBranchAddress("EDepVeto" , &_EVin);
+
+ if (t->GetBranch("PhotonsSignal"))
+ t->SetBranchAddress("PhotonsSignal", &_PSin);
+
+ if (t->GetBranch("PhotonsSignal"))
+ t->SetBranchAddress("PhotonsVeto" , &_PVin);
+
+ if (t->GetBranch("EDepV2Front"))
+ t->SetBranchAddress("EDepV2Front", &_EF2in);
+
+ if (t->GetBranch("EDepV2Signal"))
+ t->SetBranchAddress("EDepV2Signal", &_ES2in);
+
+ if (t->GetBranch("EDepV2Veto"))
+ t->SetBranchAddress("EDepV2Veto" , &_EV2in);
}
@@ -59,5 +83,8 @@ const VKEvent& VKSimData::GetEvent(size_t id)
fEvent.EV = _EVin;
fEvent.PS = _PSin;
fEvent.PV = _PVin;
+ fEvent.EF2 = _EF2in;
+ fEvent.ES2 = _ES2in;
+ fEvent.EV2 = _EV2in;
return fEvent;
}
diff --git a/analysis/src/VKSummingData.cpp b/analysis/src/VKSummingData.cpp
index a45b557ea55b0af7593a2ac13bc88fbfe1aeded5..a41dd04d62f01ee87f05c555a20daedfd2d26f9b 100644
--- a/analysis/src/VKSummingData.cpp
+++ b/analysis/src/VKSummingData.cpp
@@ -44,5 +44,8 @@ const VKEvent& VKSummingData::GetEvent(size_t id)
fEvent.EV = ev1.EV + ev2.EV;
fEvent.PS = ev1.PS + ev2.PS;
fEvent.PV = ev1.PV + ev2.PV;
+ fEvent.EF2 = ev1.EF2 + ev2.EF2;
+ fEvent.ES2 = ev1.ES2 + ev2.ES2;
+ fEvent.EV2 = ev1.EV2 + ev2.EV2;
return fEvent;
}
diff --git a/include/VKDetectorConstruction.hh b/include/VKDetectorConstruction.hh
index 176ae8e3c33619b329d4e1b87d1d8deed12a4b31..f21ccc6088047fa686c0acc85420e32dc831e177 100644
--- a/include/VKDetectorConstruction.hh
+++ b/include/VKDetectorConstruction.hh
@@ -42,6 +42,7 @@ class G4Box;
class G4Tubs;
class VKMagnet;
class VKBetaDetector;
+class VKBetaDetectorV2;
class VKMaterials;
class VKCatcherFoil;
class VKShield;
@@ -78,13 +79,14 @@ class VKDetectorConstruction : public G4VUserDetectorConstruction
* Add or remove the beta detector (v1.0) from setup
* @param b If true, detector is created. If false, detector is removed/not created
*/
- void SetBetaDetectorOn(G4bool b);
+ void SetBetaDetectorOn(G4bool b, bool defaultPos = true);
/*
* Add or remove the beta detector (v2.0) from setup
* @param b If true, detector is created. If false, detector is removed/not created
+ * @param defaultPos Return detector to default position and orientation
*/
- void SetBetaDetectorV2On(G4bool b);
+ void SetBetaDetectorV2On(G4bool b, bool defaultPos = true);
/*
* Add or remove the muon telescope from setup
@@ -164,10 +166,15 @@ class VKDetectorConstruction : public G4VUserDetectorConstruction
inline G4bool GetMagnetOn() const { return fMagnetOn; }
/*
- * Is the beta detector included in the setup?
+ * Is the beta detector v1.0 included in the setup?
*/
inline G4bool GetBetaDetectorOn() const { return fBetaDetectorOn; }
+ /*
+ * Is the beta detector v2.0 included in the setup?
+ */
+ inline G4bool GetBetaDetectorV2On() const { return fBetaDetectorV2On; }
+
/*
* Is the catcher foil included in the setup?
*/
@@ -244,10 +251,15 @@ class VKDetectorConstruction : public G4VUserDetectorConstruction
inline const VKVirtualDetector* GetMidwayPlane() const { return fMidwayPlane; }
/*
- * Get pointer to beta detector
+ * Get pointer to beta detector v1.0
*/
inline const VKBetaDetector* GetBetaDetector() const { return fBetaDetector; }
+ /*
+ * Get pointer to beta detector v2.0
+ */
+ inline const VKBetaDetectorV2* GetBetaDetectorV2() const { return fBetaDetectorV2; }
+
private:
@@ -283,11 +295,22 @@ class VKDetectorConstruction : public G4VUserDetectorConstruction
bool fMagnetOn;
VKMagnet *fMagnet;
+ G4double fRefl;
+
+ G4ThreeVector fBetaDetectorDefaultPosition;
+ G4RotationMatrix fBetaDetectorDefaultRotation;
+
// Beta detector
G4ThreeVector fBetaDetectorPosition;
+ G4RotationMatrix fBetaDetectorRotation;
bool fBetaDetectorOn;
VKBetaDetector *fBetaDetector;
- G4double fRefl;
+
+ // Beta detector V2
+ G4ThreeVector fBetaDetectorV2Position;
+ G4RotationMatrix fBetaDetectorV2Rotation;
+ bool fBetaDetectorV2On;
+ VKBetaDetectorV2 *fBetaDetectorV2;
// Catcher foil
G4ThreeVector fCatcherFoilPosition;
@@ -311,10 +334,15 @@ class VKDetectorConstruction : public G4VUserDetectorConstruction
bool fVirtualDetectorsOn;
VKVirtualDetector *fFocalPlane, *fMidwayPlane;
- // Sensitive Detectors
+ // Sensitive Detectors (BetaDetector)
G4bool fSetSDs;
G4Cache<VKScintillatorSensitiveDetector*> fBetaDetSignal_SD, fBetaDetVeto_SD;
G4Cache<VKSiPMSensitiveDetector*> fSiPM_SD;
+
+ // Sensitive Detectors (BetaDetectorV2)
+ G4bool fSetV2SDs;
+ G4Cache<VKScintillatorSensitiveDetector*> fBetaDetV2Front_SD, fBetaDetV2Signal_SD, fBetaDetV2Veto_SD;
+ G4Cache<VKSiPMSensitiveDetector*> fSiPMV2_SD;
};
#endif
diff --git a/include/VKEventAction.hh b/include/VKEventAction.hh
index ab58288e26c7d5a12281a6c8c5de6ca2813ff56d..d5b24a04456bbb4d3f72ff4503f6516d67471ed0 100644
--- a/include/VKEventAction.hh
+++ b/include/VKEventAction.hh
@@ -62,9 +62,16 @@ class VKEventAction : public G4UserEventAction
VKEventMessenger* fEventMessenger;
G4int fSaveThreshold;
+
G4int fScintCollID;
G4int fVetoCollID;
G4int fSiPMCollID;
+
+ G4int fFrontV2CollID;
+ G4int fScintV2CollID;
+ G4int fVetoV2CollID;
+ G4int fSiPMV2CollID;
+
G4int fVerbose;
G4int fSiPMThreshold;
G4bool fForcedrawphotons;
diff --git a/include/VKRecorderRoot.hh b/include/VKRecorderRoot.hh
index 25b3a5295c4bf870fc54c8adb6ca35dd8db9a2aa..77bb8aa88a6979f402c226b7c5b8de2a83a22c9a 100644
--- a/include/VKRecorderRoot.hh
+++ b/include/VKRecorderRoot.hh
@@ -117,10 +117,11 @@ class VKRecorderRoot : public VKRecorderBase {
G4String fRunMacroFileName;
G4String fFilename;
- G4bool fBetaOn, fVirtualOn;
+ G4bool fBetaOn, fBetaV2On, fVirtualOn;
const VKDetectorConstruction *fDetector;
const VKVirtualDetector *fFocalPlane, *fMidwayPlane;
const G4VPhysicalVolume *fSignalPV, *fVetoPV;
+ const G4VPhysicalVolume *fFrontV2PV, *fSignalV2PV, *fVetoV2PV;
G4PrimaryVertex *fPrimVertex;
G4PrimaryParticle *fPrimParticle;
@@ -133,6 +134,9 @@ class VKRecorderRoot : public VKRecorderBase {
G4double fEDepVeto;
G4double fPhotonsSignal;
G4double fPhotonsVeto;
+ G4double fEDepV2Front;
+ G4double fEDepV2Signal;
+ G4double fEDepV2Veto;
// ROOT output
TFile *fTFile;
@@ -144,6 +148,9 @@ class VKRecorderRoot : public VKRecorderBase {
Float_t _eDepVeto;
UInt_t _photonsSignal;
UInt_t _photonsVeto;
+ Float_t _eDepV2Front;
+ Float_t _eDepV2Signal;
+ Float_t _eDepV2Veto;
Float_t _Xfp, _Yfp, _Zfp, _Efp, _THfp, _PHIfp;
Float_t _Xmp, _Ymp, _Zmp, _Emp, _THmp, _PHImp;
diff --git a/mac/electron.mac b/mac/electron.mac
index 345713e1aab7f1bd959a3e458c7159ed4a0ab953..f48e8be592be230c79b76a281d8b2e3bc321b737 100644
--- a/mac/electron.mac
+++ b/mac/electron.mac
@@ -2,6 +2,7 @@
/control/execute mac/init.mac
/VK/eventVerbose 0
+
# output
/VK/output/zeroSuppression true
/VK/storeTrajectoryData true
@@ -18,13 +19,14 @@
/gps/ang/maxtheta 180 deg
#/gps/energy 481.697 keV
-/gps/energy 975.657 keV
-#/gps/energy 1682.232 keV
+#/gps/energy 975.657 keV
+/gps/energy 1682.232 keV
-/VK/field/maxbfield 0.12 tesla
+/VK/field/maxbfield 0.15 tesla
/VK/output/openFile output/test.root
/run/beamOn 10000
+
#/VK/field/maxbfield 0.0750 tesla
#/VK/output/openFile output/scan/e2/bmax0750.root
#/run/beamOn 10000
diff --git a/mac/init.mac b/mac/init.mac
index 8b1326855ee5d85a1fd39eab89bea558e9178d7b..91638e58a54e521252bd52c5f49f3df1bddf0ba6 100644
--- a/mac/init.mac
+++ b/mac/init.mac
@@ -10,7 +10,8 @@
/VK/physics/opticalPhysics false
# setup
-/VK/detector/volumes/betaDetector true
+/VK/detector/volumes/betaDetector false
+/VK/detector/volumes/betaDetectorV2 true
/VK/detector/volumes/magnet true
/VK/detector/volumes/catcherFoil true
/VK/detector/volumes/shield true
diff --git a/mac/muon.mac b/mac/muon.mac
index a9f3a2e90dead2f0a518577e48a4bc8330a1d04f..38e5f17a502528b68be533676cf2a6ee5b7b1e1f 100644
--- a/mac/muon.mac
+++ b/mac/muon.mac
@@ -3,7 +3,8 @@
/control/execute mac/init.mac
# simplify setup
-/VK/detector/volumes/betaDetector true
+/VK/detector/volumes/muonTelescope true
+###/VK/detector/volumes/betaDetector true
/VK/detector/volumes/magnet false
/VK/detector/volumes/catcherFoil false
/VK/detector/volumes/shield false
@@ -15,20 +16,22 @@
# source
/VK/source/cosmicRayAngDist true
/VK/source/cosmicRayMinZenithAngle 0 deg
-/VK/source/cosmicRayMaxZenithAngle 90 deg
+/VK/source/cosmicRayMaxZenithAngle 20 deg
/VK/source/diskShapedSource true
/VK/source/innerDiameter 0 cm
-/VK/source/outerDiameter 60 cm
+###/VK/source/outerDiameter 60 cm
+/VK/source/outerDiameter 5.0 cm
/VK/source/surfaceNormal 0 1 0
# radiation
/gps/particle mu-
/gps/energy 1 GeV
-/gps/position 0 3.1 32 cm
+###/gps/position 0 3.1 32 cm
+/gps/position 0 29.0 0 cm
# output
/VK/output/zeroSuppression true
/VK/output/openFile output/muon.root
# run simulations
-/run/beamOn 1000000
+/run/beamOn 100000
diff --git a/mac/visMuonTelescope.mac b/mac/visMuonTelescope.mac
new file mode 100644
index 0000000000000000000000000000000000000000..9a412fe7c5b9c1e490acb7cee87fc83ea9c6e544
--- /dev/null
+++ b/mac/visMuonTelescope.mac
@@ -0,0 +1,18 @@
+
+/vis/viewer/set/viewpointVector 0 1 0
+/vis/viewer/zoom 1.0
+
+/VK/detector/volumes/betaDetector false
+/VK/detector/volumes/betaDetectorV2 false
+/VK/detector/volumes/magnet true
+/VK/detector/volumes/catcherFoil false
+/VK/detector/volumes/shield false
+/VK/detector/volumes/blinkers false
+
+/VK/storeTrajectoryData false
+
+/VK/field/enable false
+
+/VK/detector/volumes/muonTelescope true
+
+#/run/initialize
diff --git a/src/VKBetaDetector.cpp b/src/VKBetaDetector.cpp
index 502af0d9c8a9a7266a2a55cdfab84cc98748d607..470a6afd57949765f2fd35a37e4c992b61df3556 100644
--- a/src/VKBetaDetector.cpp
+++ b/src/VKBetaDetector.cpp
@@ -251,7 +251,9 @@ VKBetaDetector::VKBetaDetector(G4RotationMatrix *pRot,
// Move the detector so that 'tlate' coincides with
// the front of the detector (including coating)
- pos = tlate + G4ThreeVector (0, 0, fContainerLength/2);
+ G4ThreeVector dx(0, 0, fContainerLength/2);
+ dx = (*pRot).inverse() * dx;
+ pos = tlate + dx;
this -> SetTranslation(pos);
// Check for overlaps
diff --git a/src/VKBetaDetectorV2.cpp b/src/VKBetaDetectorV2.cpp
index 3e149d5c0a1baeb4a4233434acf70b4baebac10f..d15873b923eafaaca9854ff657da81bc6ebbe9d5 100644
--- a/src/VKBetaDetectorV2.cpp
+++ b/src/VKBetaDetectorV2.cpp
@@ -79,7 +79,7 @@ VKBetaDetectorV2::VKBetaDetectorV2(G4RotationMatrix *pRot,
: G4PVPlacement(pRot,
tlate,
CreateTemporaryLogicalVolume(), // <-- see VKUtil.cpp
- "betaDetector",
+ "betaDetectorV2",
pMother_log,
pMany,
pCopyNo),
@@ -123,11 +123,11 @@ VKBetaDetectorV2::VKBetaDetectorV2(G4RotationMatrix *pRot,
// *** All other volumes should fit inside this volume !
// solid
- fContainer_sol = new G4Tubs("betaDetector_sol", 0., fContainerDiameter/2,
+ fContainer_sol = new G4Tubs("betaDetectorV2_sol", 0., fContainerDiameter/2,
fContainerLength/2, 0.*deg, 360.*deg);
// logical
- fContainer_log = new G4LogicalVolume(fContainer_sol, fLabVacuum, "betaDetector_log");
+ fContainer_log = new G4LogicalVolume(fContainer_sol, fLabVacuum, "betaDetectorV2_log");
// ----------------------------------------------- //
@@ -135,16 +135,16 @@ VKBetaDetectorV2::VKBetaDetectorV2(G4RotationMatrix *pRot,
// ----------------------------------------------- //
// solid
- fFront_sol = new G4Tubs("front_sol", 0., fFrontDiameter/2,
+ fFront_sol = new G4Tubs("frontV2_sol", 0., fFrontDiameter/2,
fFrontLength/2, 0.*deg, 360.*deg);
// logical
- fFront_log = new G4LogicalVolume(fFront_sol, fEJ200, "front_log");
+ fFront_log = new G4LogicalVolume(fFront_sol, fEJ200, "frontV2_log");
// physical
displ = -(fContainerLength - fFrontLength)/2 + fReflectorThickness;
pos = G4ThreeVector(0, 0, displ);
- new G4PVPlacement(0, pos, fFront_log, "front", fContainer_log, false, 0, fCheckOverlaps);
+ new G4PVPlacement(0, pos, fFront_log, "frontV2", fContainer_log, false, 0, fCheckOverlaps);
// ----------------------------------------------- //
@@ -152,17 +152,17 @@ VKBetaDetectorV2::VKBetaDetectorV2(G4RotationMatrix *pRot,
// ----------------------------------------------- //
// solid
- fSignal_sol = new G4Tubs("signal_sol", 0., fSignalDiameter/2,
+ fSignal_sol = new G4Tubs("signalV2_sol", 0., fSignalDiameter/2,
fSignalLength/2, 0.*deg, 360.*deg);
// logical
- fSignal_log = new G4LogicalVolume(fSignal_sol, fEJ200, "signal_log");
+ fSignal_log = new G4LogicalVolume(fSignal_sol, fEJ200, "signalV2_log");
// physical
displ = -(fContainerLength - fSignalLength)/2 + fReflectorThickness
+ fFrontLength; //+ fReflectorThickness;
pos = G4ThreeVector(0, 0, displ);
- new G4PVPlacement(0, pos, fSignal_log, "signal", fContainer_log, false, 0, fCheckOverlaps);
+ new G4PVPlacement(0, pos, fSignal_log, "signalV2", fContainer_log, false, 0, fCheckOverlaps);
// ----------------------------------------------- //
@@ -173,10 +173,10 @@ VKBetaDetectorV2::VKBetaDetectorV2(G4RotationMatrix *pRot,
outerRadius = innerRadius + fVetoThickness;
// tube
- G4Tubs* vetoTube_sol = new G4Tubs("vetoTube_sol", innerRadius, outerRadius,
+ G4Tubs* vetoTube_sol = new G4Tubs("vetoTubeV2_sol", innerRadius, outerRadius,
fVetoLength/2, 0.*deg, 360.*deg);
// end cap
- G4Tubs* vetoEndcap_sol = new G4Tubs("vetoEndcap_sol", 0., innerRadius,
+ G4Tubs* vetoEndcap_sol = new G4Tubs("vetoEndcapV2_sol", 0., innerRadius,
fVetoEndcapLength/2, 0.*deg, 360.*deg);
// relative placement
displ = fVetoLength/2 - fVetoEndcapLength/2;
@@ -184,15 +184,15 @@ VKBetaDetectorV2::VKBetaDetectorV2(G4RotationMatrix *pRot,
transform = G4Transform3D(doNotRotate, pos);
// create union
- fVeto_sol = new G4UnionSolid("veto_sol", vetoTube_sol, vetoEndcap_sol, transform);
+ fVeto_sol = new G4UnionSolid("vetoV2_sol", vetoTube_sol, vetoEndcap_sol, transform);
// logical
- fVeto_log = new G4LogicalVolume(fVeto_sol, fEJ200, "veto_log");
+ fVeto_log = new G4LogicalVolume(fVeto_sol, fEJ200, "vetoV2_log");
// physical
displ = -(fContainerLength - fVetoLength)/2 + fReflectorThickness;
pos = G4ThreeVector(0, 0, displ);
- new G4PVPlacement(0, pos, fVeto_log, "veto", fContainer_log, false, 0, fCheckOverlaps);
+ new G4PVPlacement(0, pos, fVeto_log, "vetoV2", fContainer_log, false, 0, fCheckOverlaps);
// ----------------------------------------------- //
@@ -200,16 +200,16 @@ VKBetaDetectorV2::VKBetaDetectorV2(G4RotationMatrix *pRot,
// ----------------------------------------------- //
// solid
- fWindow_sol = new G4Tubs("window_sol", 0., fSignalDiameter/2,
+ fWindow_sol = new G4Tubs("windowV2_sol", 0., fSignalDiameter/2,
fWindowThickness/2, 0.*deg, 360.*deg);
// logical
- fWindow_log = new G4LogicalVolume(fWindow_sol, fMylar, "window_log");
+ fWindow_log = new G4LogicalVolume(fWindow_sol, fMylar, "windowV2_log");
// physical
displ = -(fContainerLength + fWindowThickness)/2 + fReflectorThickness;
pos = G4ThreeVector(0, 0, displ);
- new G4PVPlacement(0, pos, fWindow_log, "window", fContainer_log, false, 0, fCheckOverlaps);
+ new G4PVPlacement(0, pos, fWindow_log, "windowV2", fContainer_log, false, 0, fCheckOverlaps);
// ----------------------------------------------- //
@@ -220,11 +220,11 @@ VKBetaDetectorV2::VKBetaDetectorV2(G4RotationMatrix *pRot,
ConstructCoatingSolid(displ);
// logical
- fCoating_log = new G4LogicalVolume(fCoating_sol, fTeflon, "coating_log");
+ fCoating_log = new G4LogicalVolume(fCoating_sol, fTeflon, "coatingV2_log");
// physical
pos = G4ThreeVector(0, 0, displ);
- new G4PVPlacement(0, pos, fCoating_log, "coating", fContainer_log, false, 0, fCheckOverlaps);
+ new G4PVPlacement(0, pos, fCoating_log, "coatingV2", fContainer_log, false, 0, fCheckOverlaps);
// ----------------------------------------------- //
@@ -232,8 +232,8 @@ VKBetaDetectorV2::VKBetaDetectorV2(G4RotationMatrix *pRot,
// ----------------------------------------------- //
// SiPM solid and logical volumes
- fSiPM_sol = new G4Tubs("SiPM_sol", 0., fSiPMDiameter/2, fSiPMLength/2, 0.*deg, 360.*deg);
- fSiPM_log = new G4LogicalVolume(fSiPM_sol, fGlass, "SiPM_log");
+ fSiPM_sol = new G4Tubs("SiPMV2_sol", 0., fSiPMDiameter/2, fSiPMLength/2, 0.*deg, 360.*deg);
+ fSiPM_log = new G4LogicalVolume(fSiPM_sol, fGlass, "SiPMV2_log");
// Photocathode volumes
//
@@ -241,23 +241,23 @@ VKBetaDetectorV2::VKBetaDetectorV2(G4RotationMatrix *pRot,
// is only a very rough approximation of the real thing since it only
// absorbs or detects the photons based on the efficiency set below.
//
- fPhotocathode_sol = new G4Tubs("photocathode_sol", 0., fSiPMDiameter/2, fSiPMLength/4, 0.*deg, 360.*deg);
- fPhotocathode_log = new G4LogicalVolume(fPhotocathode_sol, fAluminium, "photocathode_log");
+ fPhotocathode_sol = new G4Tubs("photocathodeV2_sol", 0., fSiPMDiameter/2, fSiPMLength/4, 0.*deg, 360.*deg);
+ fPhotocathode_log = new G4LogicalVolume(fPhotocathode_sol, fAluminium, "photocathodeV2_log");
pos = G4ThreeVector(0, 0, fSiPMLength/4);
- new G4PVPlacement(0, pos, fPhotocathode_log, "photocathode", fSiPM_log, false, fCheckOverlaps);
+ new G4PVPlacement(0, pos, fPhotocathode_log, "photocathodeV2", fSiPM_log, false, fCheckOverlaps);
// Place one SiPM on the back side of the SIGNAL scintillator
displ = -(fContainerLength - fSiPMLength)/2 + fReflectorThickness + fSignalLength;
pos = G4ThreeVector(0, 0, displ);
copyNo = 0;
- new G4PVPlacement(0, pos, fSiPM_log, "SiPM", fContainer_log, false, copyNo, fCheckOverlaps);
+ new G4PVPlacement(0, pos, fSiPM_log, "SiPMV2", fContainer_log, false, copyNo, fCheckOverlaps);
fSiPMPositions.push_back(pos);
// Place another SiPM on the back side of the VETO scintillator
displ = -(fContainerLength - fSiPMLength)/2 + fReflectorThickness + fVetoLength;
pos = G4ThreeVector(0, 0, displ);
copyNo = 1;
- new G4PVPlacement(0, pos, fSiPM_log, "SiPM", fContainer_log, false, copyNo, fCheckOverlaps);
+ new G4PVPlacement(0, pos, fSiPM_log, "SiPMV2", fContainer_log, false, copyNo, fCheckOverlaps);
fSiPMPositions.push_back(pos);
// ***IMPORTANT***
@@ -280,7 +280,9 @@ VKBetaDetectorV2::VKBetaDetectorV2(G4RotationMatrix *pRot,
// Move the detector so that 'tlate' coincides with
// the front of the detector (including coating)
- pos = tlate + G4ThreeVector (0, 0, fContainerLength/2);
+ G4ThreeVector dx(0, 0, fContainerLength/2);
+ dx = (*pRot).inverse() * dx;
+ pos = tlate + dx;
this -> SetTranslation(pos);
// Check for overlaps
@@ -367,7 +369,7 @@ void VKBetaDetectorV2::SurfaceProperties()
photocath_mpt -> AddProperty("IMAGINARYRINDEX", ephoton, photocath_ImR, num);
// Create optical surface
G4OpticalSurface* photocathOptSurf =
- new G4OpticalSurface("photocathOptSurf", glisur, polished, dielectric_metal);
+ new G4OpticalSurface("photocathOptSurfV2", glisur, polished, dielectric_metal);
photocathOptSurf -> SetMaterialPropertiesTable(photocath_mpt);
//
@@ -380,7 +382,7 @@ void VKBetaDetectorV2::SurfaceProperties()
mylarWindow_mpt -> AddProperty("EFFICIENCY", ephoton, mylarEff, num);
// Create optical surface
G4OpticalSurface* mylarWindowOptSurf =
- new G4OpticalSurface("mylarWindowOptSurf", unified, polished, dielectric_metal);
+ new G4OpticalSurface("mylarWindowOptSurfV2", unified, polished, dielectric_metal);
mylarWindowOptSurf -> SetMaterialPropertiesTable(mylarWindow_mpt);
@@ -389,9 +391,9 @@ void VKBetaDetectorV2::SurfaceProperties()
//
// *** NB: Mylar window gets same properties as reflective coating
//
- new G4LogicalSkinSurface("coating_surf", fCoating_log, reflCoatOptSurf);
- new G4LogicalSkinSurface("window_surf", fWindow_log, mylarWindowOptSurf);
- new G4LogicalSkinSurface("photocathode_surf", fPhotocathode_log, photocathOptSurf);
+ new G4LogicalSkinSurface("coatingV2_surf", fCoating_log, reflCoatOptSurf);
+ new G4LogicalSkinSurface("windowV2_surf", fWindow_log, mylarWindowOptSurf);
+ new G4LogicalSkinSurface("photocathodeV2_surf", fPhotocathode_log, photocathOptSurf);
}
@@ -408,7 +410,7 @@ void VKBetaDetectorV2::ConstructCoatingSolid( G4double& displ )
G4double r1min = fSignalDiameter/2;
G4double r1max = fSignalDiameter/2 + fReflectorThickness;
G4double l1 = fReflectorThickness + fVetoLength - fVetoEndcapLength;
- G4Tubs* c1 = new G4Tubs("coating1_sol", r1min, r1max, l1/2, 0.*deg, 360.*deg);
+ G4Tubs* c1 = new G4Tubs("coating1V2_sol", r1min, r1max, l1/2, 0.*deg, 360.*deg);
// this will determine placement of coating volumes relative to other volumes
displ = -(fContainerLength - l1)/2;
@@ -417,7 +419,7 @@ void VKBetaDetectorV2::ConstructCoatingSolid( G4double& displ )
G4double r2min = fSiPMDiameter/2;
G4double r2max = r1min;
G4double l2 = fReflectorThickness;
- G4Tubs* c2 = new G4Tubs("coating2_sol", r2min, r2max, l2/2, 0.*deg, 360.*deg);
+ G4Tubs* c2 = new G4Tubs("coating2V2_sol", r2min, r2max, l2/2, 0.*deg, 360.*deg);
// relative placement of 1 and 2
dz = (l2-l1)/2 + fSignalLength + fFrontLength + l2;
@@ -425,13 +427,13 @@ void VKBetaDetectorV2::ConstructCoatingSolid( G4double& displ )
transform = G4Transform3D(doNotRotate, translation);
// create union of 1 and 2
- fCoating_sol = new G4UnionSolid("coating12_sol", c1, c2, transform);
+ fCoating_sol = new G4UnionSolid("coating12V2_sol", c1, c2, transform);
// coating 3
G4double r3min = 0.;
G4double r3max = r1min;
G4double l3 = fReflectorThickness;
- G4Tubs* c3 = new G4Tubs("coating3_sol", r3min, r3max, l3/2, 0.*deg, 360.*deg);
+ G4Tubs* c3 = new G4Tubs("coating3V2_sol", r3min, r3max, l3/2, 0.*deg, 360.*deg);
// relative placement of 1 and 3
dz = (l1-l3)/2;
@@ -439,13 +441,13 @@ void VKBetaDetectorV2::ConstructCoatingSolid( G4double& displ )
transform = G4Transform3D(doNotRotate, translation);
// create union of 12 and 3
- fCoating_sol = new G4UnionSolid("coating123_sol", fCoating_sol, c3, transform);
+ fCoating_sol = new G4UnionSolid("coating123V2_sol", fCoating_sol, c3, transform);
// coating 4
G4double r4min = r1max;
G4double r4max = r4min + fVetoThickness;
G4double l4 = fReflectorThickness;
- G4Tubs* c4 = new G4Tubs("coating4_sol", r4min, r4max, l4/2, 0.*deg, 360.*deg);
+ G4Tubs* c4 = new G4Tubs("coating4V2_sol", r4min, r4max, l4/2, 0.*deg, 360.*deg);
// relative placement of 1 and 4
dz = (l4-l1)/2;
@@ -453,13 +455,13 @@ void VKBetaDetectorV2::ConstructCoatingSolid( G4double& displ )
transform = G4Transform3D(doNotRotate, translation);
// create union of 123 and 4
- fCoating_sol = new G4UnionSolid("coating1234_sol", fCoating_sol, c4, transform);
+ fCoating_sol = new G4UnionSolid("coating1234V2_sol", fCoating_sol, c4, transform);
// coating 5
G4double r5min = r4max;
G4double r5max = r5min + fReflectorThickness;
G4double l5 = fVetoLength + 2*fReflectorThickness;
- G4Tubs* c5 = new G4Tubs("coating5_sol", r5min, r5max, l5/2, 0.*deg, 360.*deg);
+ G4Tubs* c5 = new G4Tubs("coating5V2_sol", r5min, r5max, l5/2, 0.*deg, 360.*deg);
// relative placement of 1 and 5
dz = (l5-l1)/2;
@@ -467,13 +469,13 @@ void VKBetaDetectorV2::ConstructCoatingSolid( G4double& displ )
transform = G4Transform3D(doNotRotate, translation);
// create union of 1234 and 5
- fCoating_sol = new G4UnionSolid("coating12345_sol", fCoating_sol, c5, transform);
+ fCoating_sol = new G4UnionSolid("coating12345V2_sol", fCoating_sol, c5, transform);
// coating 6
G4double r6min = fSiPMDiameter/2;
G4double r6max = r5min;
G4double l6 = fReflectorThickness;
- G4Tubs* c6 = new G4Tubs("coating6_sol", r6min, r6max, l6/2, 0.*deg, 360.*deg);
+ G4Tubs* c6 = new G4Tubs("coating6V2_sol", r6min, r6max, l6/2, 0.*deg, 360.*deg);
// relative placement of 1 and 6
dz = (l1-l6)/2 + l5-l1;
@@ -481,7 +483,7 @@ void VKBetaDetectorV2::ConstructCoatingSolid( G4double& displ )
transform = G4Transform3D(doNotRotate, translation);
// create union of 12345 and 6
- fCoating_sol = new G4UnionSolid("coating_sol", fCoating_sol, c6, transform);
+ fCoating_sol = new G4UnionSolid("coatingV2_sol", fCoating_sol, c6, transform);
}
@@ -489,7 +491,7 @@ const G4VPhysicalVolume* VKBetaDetectorV2::GetFrontPhysVol() const
{
const G4PhysicalVolumeStore *physVolStore =
static_cast<const G4PhysicalVolumeStore*>(G4PhysicalVolumeStore::GetInstance());
- return physVolStore -> GetVolume("front", false);
+ return physVolStore -> GetVolume("frontV2", false);
}
@@ -497,7 +499,7 @@ const G4VPhysicalVolume* VKBetaDetectorV2::GetSignalPhysVol() const
{
const G4PhysicalVolumeStore *physVolStore =
static_cast<const G4PhysicalVolumeStore*>(G4PhysicalVolumeStore::GetInstance());
- return physVolStore -> GetVolume("signal", false);
+ return physVolStore -> GetVolume("signalV2", false);
}
@@ -505,5 +507,5 @@ const G4VPhysicalVolume* VKBetaDetectorV2::GetVetoPhysVol() const
{
const G4PhysicalVolumeStore *physVolStore =
static_cast<const G4PhysicalVolumeStore*>(G4PhysicalVolumeStore::GetInstance());
- return physVolStore -> GetVolume("veto", false);
+ return physVolStore -> GetVolume("vetoV2", false);
}
diff --git a/src/VKDetectorConstruction.cpp b/src/VKDetectorConstruction.cpp
index 4295ff7f961884bfab4c0572b48c104eaf3b7742..f4f6e2c23ec37f3238848ef4baf4e674e3d110fc 100644
--- a/src/VKDetectorConstruction.cpp
+++ b/src/VKDetectorConstruction.cpp
@@ -78,6 +78,7 @@ VKDetectorConstruction::VKDetectorConstruction()
fMagnet = NULL;
fBetaDetector = NULL;
+ fBetaDetectorV2 = NULL;
fCatcherFoil = NULL;
fShield = NULL;
fLeadBlock = NULL;
@@ -89,10 +90,15 @@ VKDetectorConstruction::VKDetectorConstruction()
fBField.Put(NULL);
fSiPM_SD.Put(NULL);
+ fSiPMV2_SD.Put(NULL);
fBetaDetSignal_SD.Put(NULL);
fBetaDetVeto_SD.Put(NULL);
+ fBetaDetV2Front_SD.Put(NULL);
+ fBetaDetV2Signal_SD.Put(NULL);
+ fBetaDetV2Veto_SD.Put(NULL);
fSetSDs = true;
+ fSetV2SDs = true;
SetDefaults();
@@ -150,7 +156,12 @@ G4VPhysicalVolume* VKDetectorConstruction::ConstructDetector()
// Create beta detector
if (fBetaDetectorOn) {
- fBetaDetector = new VKBetaDetector(0, fBetaDetectorPosition, mother_log, false, 0, this);
+ fBetaDetector = new VKBetaDetector(&fBetaDetectorRotation, fBetaDetectorPosition, mother_log, false, 0, this);
+ }
+
+ // Create beta detector V2
+ if (fBetaDetectorV2On) {
+ fBetaDetectorV2 = new VKBetaDetectorV2(&fBetaDetectorV2Rotation, fBetaDetectorV2Position, mother_log, false, 0, this);
}
// Create catcher foil
@@ -293,6 +304,74 @@ void VKDetectorConstruction::ConstructSDandField()
SetSensitiveDetector(fBetaDetector->GetSignalLogVol(), fBetaDetSignal_SD.Get());
SetSensitiveDetector(fBetaDetector->GetVetoLogVol(), fBetaDetVeto_SD.Get());
}
+
+
+ //
+ // Create sensitive detectors only if beta detector V2 is part of setup
+ //
+ if (fBetaDetectorV2On)
+ {
+ // SiPMs
+ if (!fSiPMV2_SD.Get()) { // create only once
+ // Created here so it exists as SiPMs are being placed
+ G4cout << " ---> Constructing SiPMV2_SD" << G4endl;
+ VKSiPMSensitiveDetector* sipm_SD = new VKSiPMSensitiveDetector("SiPMV2_SD");
+ fSiPMV2_SD.Put(sipm_SD);
+ const std::vector<G4ThreeVector> positions = fBetaDetectorV2->GetSiPMPositions();
+ sipm_SD -> InitSiPMs(positions.size()); // number of SiPMs
+ sipm_SD -> SetSiPMPositions(positions); // their positions
+ }
+ // Sensitive detector is not actually on the photocathode.
+ // processHits gets done manually by the stepping action.
+ // It is used to detect when photons hit and get absorbed&detected at the
+ // boundary to the photocathode (which doesnt get done by attaching it to a
+ // logical volume.
+ // It does however need to be attached to something or else it doesnt get
+ // reset at the begining of events.
+
+ // FRONT scintillator detector
+ if (!fBetaDetV2Front_SD.Get()) { // create only once
+ G4cout << " ---> Constructing FrontV2_SD" << G4endl;
+ VKScintillatorSensitiveDetector* front_SD = new VKScintillatorSensitiveDetector("FrontV2_SD");
+ fBetaDetV2Front_SD.Put(front_SD);
+ }
+
+ // SIGNAL scintillator detector
+ if (!fBetaDetV2Signal_SD.Get()) { // create only once
+ G4cout << " ---> Constructing SignalV2_SD" << G4endl;
+ VKScintillatorSensitiveDetector* signal_SD = new VKScintillatorSensitiveDetector("SignalV2_SD");
+ fBetaDetV2Signal_SD.Put(signal_SD);
+ }
+
+ // VETO scintillator detector
+ if (!fBetaDetV2Veto_SD.Get()) { // create only once
+ G4cout << " ---> Constructing VetoV2_SD" << G4endl;
+ VKScintillatorSensitiveDetector* veto_SD = new VKScintillatorSensitiveDetector("VetoV2_SD");
+ fBetaDetV2Veto_SD.Put(veto_SD);
+ }
+
+ if (fSetV2SDs) {
+ // set verbosity
+ G4SDManager::GetSDMpointer()->SetVerboseLevel(0);
+
+ // register sensitive detectors with SD manager
+ G4SDManager* SDman = G4SDManager::GetSDMpointer();
+ SDman -> AddNewDetector(fSiPMV2_SD.Get());
+ SDman -> AddNewDetector(fBetaDetV2Front_SD.Get());
+ SDman -> AddNewDetector(fBetaDetV2Signal_SD.Get());
+ SDman -> AddNewDetector(fBetaDetV2Veto_SD.Get());
+
+ // Only add detectors once
+ fSetV2SDs = false;
+ }
+
+ // register sensitive detectors to logical volumes
+ // (repeat every time geometry is reinitialized)
+ SetSensitiveDetector(fBetaDetectorV2->GetPhotocathodeLogVol(), fSiPMV2_SD.Get());
+ SetSensitiveDetector(fBetaDetectorV2->GetFrontLogVol(), fBetaDetV2Front_SD.Get());
+ SetSensitiveDetector(fBetaDetectorV2->GetSignalLogVol(), fBetaDetV2Signal_SD.Get());
+ SetSensitiveDetector(fBetaDetectorV2->GetVetoLogVol(), fBetaDetV2Veto_SD.Get());
+ }
}
@@ -304,9 +383,15 @@ void VKDetectorConstruction::SetMagnetOn(G4bool b)
}
-void VKDetectorConstruction::SetBetaDetectorOn(G4bool b)
+void VKDetectorConstruction::SetBetaDetectorOn(G4bool b, bool defaultPos)
{
+ if (b && defaultPos) {
+ fBetaDetectorPosition = fBetaDetectorDefaultPosition;
+ fBetaDetectorRotation = fBetaDetectorDefaultRotation;
+ }
+
if (!fBetaDetectorOn && b) fSetSDs = true;
+
fBetaDetectorOn = b;
// G4SDManager::GetSDMpointer()->Activate("/VKDet", b); // activate/deactivate sensitive detectors
G4RunManager::GetRunManager()->ReinitializeGeometry();
@@ -314,15 +399,33 @@ void VKDetectorConstruction::SetBetaDetectorOn(G4bool b)
}
-void VKDetectorConstruction::SetBetaDetectorV2On(G4bool b)
+void VKDetectorConstruction::SetBetaDetectorV2On(G4bool b, bool defaultPos)
{
- // not implemented yet !!!
+ if (b && defaultPos) {
+ fBetaDetectorV2Position = fBetaDetectorDefaultPosition;
+ fBetaDetectorV2Rotation = fBetaDetectorDefaultRotation;
+ }
+
+ if (!fBetaDetectorV2On && b) fSetV2SDs = true;
+
+ fBetaDetectorV2On = b;
+ G4RunManager::GetRunManager()->ReinitializeGeometry();
+ G4RunManager::GetRunManager()->Initialize();
}
void VKDetectorConstruction::SetMuonTelescopeOn(G4bool b)
{
- // not implemented yet !!!
+ // detector positions
+ fBetaDetectorPosition = G4ThreeVector(0, -8.5*cm, 0);
+ fBetaDetectorV2Position = G4ThreeVector(0, 28.5*cm, 0);
+
+ // rotate detectors to vertical
+ fBetaDetectorRotation.rotateX(-90*deg);
+ fBetaDetectorV2Rotation.rotateX(-90*deg);
+
+ SetBetaDetectorOn(b, false);
+ SetBetaDetectorV2On(b, false);
}
@@ -410,11 +513,23 @@ void VKDetectorConstruction::SetDefaults()
fSpectrometerLength = 80.*cm;
fIronEndcapThickness = 4.5*cm;
fIronConeLength = 5.5*cm;
-
- // Beta detector ON/OFF and position
- fBetaDetectorOn = true;
- fBetaDetectorPosition = G4ThreeVector(0, 0, 30.0*cm);
+
+ // Coating reflectivity for beta detectors
fRefl = 0.99; // from http://dx.doi.org/10.1109/TNS.2008.20032531.0
+
+ // default beta detector detector position
+ fBetaDetectorDefaultPosition = G4ThreeVector(0, 0, 30.0*cm);
+ fBetaDetectorDefaultRotation.rotateX(0.);
+
+ // Beta detector v1.0 ON/OFF and position
+ fBetaDetectorOn = true;
+ fBetaDetectorPosition = fBetaDetectorDefaultPosition;
+ fBetaDetectorRotation = fBetaDetectorDefaultRotation;
+
+ // Beta detector v2.0 ON/OFF and position
+ fBetaDetectorV2On = false;
+ fBetaDetectorV2Position = fBetaDetectorDefaultPosition;
+ fBetaDetectorV2Rotation = fBetaDetectorDefaultRotation;
// Catcher foil ON/OFF and position
fCatcherFoilOn = true;
diff --git a/src/VKEventAction.cpp b/src/VKEventAction.cpp
index 96bc38babaf6b96e575d7c587ffc296896f914f0..470255584bb2cf5738ae212871a2f3e6709719a4 100644
--- a/src/VKEventAction.cpp
+++ b/src/VKEventAction.cpp
@@ -56,7 +56,11 @@ VKEventAction::VKEventAction(VKRecorderBase* r)
fScintCollID(-1),
fVetoCollID(-1),
fSiPMCollID(-1),
- fVerbose(0),
+ fFrontV2CollID(-1),
+ fScintV2CollID(-1),
+ fVetoV2CollID(-1),
+ fSiPMV2CollID(-1),
+ fVerbose(1),
fSiPMThreshold(1),
fForcedrawphotons(false),
fForcenophotons(false),
@@ -74,13 +78,25 @@ VKEventAction::~VKEventAction()
void VKEventAction::BeginOfEventAction(const G4Event* anEvent)
{
+ // current event no.
+ double evt = G4RunManager::GetRunManager()
+ ->GetCurrentEvent()->GetEventID();
+
// New event, add the user information object
G4EventManager::GetEventManager() -> SetUserInformation(new VKUserEventInformation);
G4SDManager* SDman = G4SDManager::GetSDMpointer();
- if (fScintCollID < 0) fScintCollID = SDman -> GetCollectionID("Signal_SD_HitCollection");
- if (fVetoCollID < 0) fVetoCollID = SDman -> GetCollectionID("Veto_SD_HitCollection");
- if (fSiPMCollID < 0) fSiPMCollID = SDman -> GetCollectionID("SiPM_SD_HitCollection");
+
+ if (evt == 0) { // only get collection IDs once
+ if (fScintCollID < 0) fScintCollID = SDman -> GetCollectionID("Signal_SD_HitCollection");
+ if (fVetoCollID < 0) fVetoCollID = SDman -> GetCollectionID("Veto_SD_HitCollection");
+ if (fSiPMCollID < 0) fSiPMCollID = SDman -> GetCollectionID("SiPM_SD_HitCollection");
+
+ if (fFrontV2CollID < 0) fFrontV2CollID = SDman -> GetCollectionID("FrontV2_SD_HitCollection");
+ if (fScintV2CollID < 0) fScintV2CollID = SDman -> GetCollectionID("SignalV2_SD_HitCollection");
+ if (fVetoV2CollID < 0) fVetoV2CollID = SDman -> GetCollectionID("VetoV2_SD_HitCollection");
+ if (fSiPMV2CollID < 0) fSiPMV2CollID = SDman -> GetCollectionID("SiPMV2_SD_HitCollection");
+ }
if (fRecorder) fRecorder -> RecordBeginOfEvent(anEvent);
}
@@ -114,6 +130,12 @@ void VKEventAction::EndOfEventAction(const G4Event* anEvent)
VKScintillatorHitsCollection * scintHC = 0;
VKScintillatorHitsCollection * vetoHC = 0;
VKSiPMHitsCollection* sipmHC = 0;
+
+ VKScintillatorHitsCollection * frontV2HC = 0;
+ VKScintillatorHitsCollection * scintV2HC = 0;
+ VKScintillatorHitsCollection * vetoV2HC = 0;
+ VKSiPMHitsCollection* sipmV2HC = 0;
+
G4HCofThisEvent* hitsCE = anEvent -> GetHCofThisEvent();
// Get hit collections
@@ -121,10 +143,15 @@ void VKEventAction::EndOfEventAction(const G4Event* anEvent)
if (fScintCollID >= 0) scintHC = (VKScintillatorHitsCollection*)(hitsCE->GetHC(fScintCollID));
if (fVetoCollID >= 0) vetoHC = (VKScintillatorHitsCollection*)(hitsCE->GetHC(fVetoCollID));
if (fSiPMCollID >= 0) sipmHC = (VKSiPMHitsCollection*)(hitsCE->GetHC(fSiPMCollID));
+
+ if (fFrontV2CollID >= 0) frontV2HC = (VKScintillatorHitsCollection*)(hitsCE->GetHC(fFrontV2CollID));
+ if (fScintV2CollID >= 0) scintV2HC = (VKScintillatorHitsCollection*)(hitsCE->GetHC(fScintV2CollID));
+ if (fVetoV2CollID >= 0) vetoV2HC = (VKScintillatorHitsCollection*)(hitsCE->GetHC(fVetoV2CollID));
+ if (fSiPMV2CollID >= 0) sipmV2HC = (VKSiPMHitsCollection*)(hitsCE->GetHC(fSiPMV2CollID));
}
// Hits in scintillators (signal and veto)
- vector<VKScintillatorHitsCollection*> scintHCs = {scintHC, vetoHC};
+ vector<VKScintillatorHitsCollection*> scintHCs = {scintHC, vetoHC, frontV2HC, scintV2HC, vetoV2HC};
G4int id = 0;
for (auto& hc : scintHCs) {
if (hc) {
@@ -162,8 +189,9 @@ void VKEventAction::EndOfEventAction(const G4Event* anEvent)
G4cout << "\tTotal energy deposition in " << hc->GetName() << " : "
<< eventInfo->GetEDep(id) / keV << " (keV)" << G4endl;
}
+
+ id++;
}
- id++;
}
G4String scintName[2] = {"Signal", "Veto"};
@@ -228,6 +256,71 @@ void VKEventAction::EndOfEventAction(const G4Event* anEvent)
sipmHC -> DrawAllHits();
}
+
+
+ scintName[0] = "SignalV2";
+ scintName[1] = "VetoV2";
+
+ // Hits in SiPMV2
+ if (sipmV2HC) {
+ G4ThreeVector reconPos(0.,0.,0.);
+ G4int sipms = sipmV2HC->entries(); // number of hits
+
+ //Gather info from all SiPMs
+ for (G4int i=0; i<sipms; i++) {
+ G4int id = (*sipmV2HC)[i]->GetSiPMNumber();
+ G4String name = (id < 2)? scintName[id] : "Unknown";
+ if(fVerbose > 0) {
+ G4cout << "\t======= SiPM data for: " << name << " =======" << G4endl;
+ }
+
+ G4int photons = (*sipmV2HC)[i]->GetPhotonCount();
+ eventInfo -> IncHitCount(id, photons);
+ reconPos += (*sipmV2HC)[i]->GetSiPMPos() * photons;
+ if (photons >= fSiPMThreshold) {
+ eventInfo -> IncSiPMsAboveThreshold(id);
+ }
+ else { // if it wasnt above the threshold, turn it back off
+ (*sipmV2HC)[i]->SetDrawit(false);
+ }
+
+ // if there were any hits, determine position
+ if (eventInfo->GetHitCount(id) > 0) {
+ reconPos /= eventInfo->GetHitCount(id);
+ if(fVerbose > 0) {
+ G4cout << "\tReconstructed position of hits in SiPMs: "
+ << reconPos/mm << G4endl;
+ }
+ eventInfo -> SetReconPos(id, reconPos);
+ }
+
+ //End of event output. later to be controlled by a verbose level
+ if (fVerbose > 0) {
+ G4cout << "\tNumber of photons that hit SiPMs in this event : "
+ << eventInfo->GetHitCount(id) << G4endl;
+ G4cout << "\tNumber of SiPMs above threshold ("<<fSiPMThreshold<<") : "
+ << eventInfo->GetSiPMsAboveThreshold(id) << G4endl;
+ G4cout << "\tNumber of photons produced by scintillation in this event : "
+ << eventInfo->GetPhotonCount_Scint(id) << G4endl;
+ G4cout << "\tNumber of photons produced by cerenkov in this event : "
+ << eventInfo->GetPhotonCount_Ceren(id) << G4endl;
+ G4cout << "\tNumber of photons absorbed (OpAbsorption) in this event : "
+ << eventInfo->GetAbsorptionCount(id) << G4endl;
+ G4cout << "\tNumber of photons absorbed at boundaries (OpBoundary) in "
+ << "this event : " << eventInfo->GetBoundaryAbsorptionCount(id)
+ << G4endl;
+ G4cout << "Unacounted for photons in this event : "
+ << (eventInfo->GetPhotonCount_Scint(id) +
+ eventInfo->GetPhotonCount_Ceren(id) -
+ eventInfo->GetAbsorptionCount(id) -
+ eventInfo->GetHitCount(id) -
+ eventInfo->GetBoundaryAbsorptionCount(id))
+ << G4endl;
+ }
+ }
+
+ sipmV2HC -> DrawAllHits();
+ }
// If we have set the flag to save 'special' events, save here
if (fSaveThreshold && eventInfo->GetPhotonCount() <= fSaveThreshold)
diff --git a/src/VKPhysicsList.cpp b/src/VKPhysicsList.cpp
index 2e904f75989e0b59ee226ba15dd326031f9ebcab..5b49c2cb07a5e9beb5fbf713db0da6c8d2c9205a 100644
--- a/src/VKPhysicsList.cpp
+++ b/src/VKPhysicsList.cpp
@@ -57,8 +57,9 @@ VKPhysicsList::VKPhysicsList()
RegisterPhysics(new VKGeneralPhysics("general"));
// EM Physics
- // RegisterPhysics(new VKEMPhysics("standard EM"));
- RegisterPhysics(new G4EmStandardPhysics_option4());
+ RegisterPhysics(new VKEMPhysics("standard EM"));
+ // RegisterPhysics(new G4EmStandardPhysics_option4()); // overpredicts muon energy loss by a factor of 2!
+
// RegisterPhysics(new G4EmLivermorePhysics()); // overpredicts muon energy loss by a factor of 2!
// RegisterPhysics(new G4EmLowEPPhysics());
// RegisterPhysics(new G4EmPenelopePhysics()); // not yet tested
diff --git a/src/VKRecorderRoot.cpp b/src/VKRecorderRoot.cpp
index 8cb1e113f4eda56580736f92e3c326dcd3fec819..2bb4d5568e5f0f4defd0a1de888bcb51d19106d9 100644
--- a/src/VKRecorderRoot.cpp
+++ b/src/VKRecorderRoot.cpp
@@ -27,6 +27,7 @@
#include "VKDetectorConstruction.hh"
#include "VKRecorderRoot.hh"
#include "VKBetaDetector.hh"
+#include "VKBetaDetectorV2.hh"
#include "VKRecorderRootMessenger.hh"
#include "VKPrimaryGeneratorAction.hh"
#include "VKUserEventInformation.hh"
@@ -56,6 +57,7 @@ VKRecorderRoot::VKRecorderRoot(G4String name)
fFocalPlane = fMidwayPlane = NULL;
fSignalPV = fVetoPV = NULL;
+ fFrontV2PV = fSignalV2PV = fVetoV2PV = NULL;
// create messenger class
fRecorderRoot_msg = new VKRecorderRootMessenger(this);
@@ -82,13 +84,21 @@ void VKRecorderRoot::RecordBeginOfRun(const G4Run*)
// which detectors are in the setup?
fBetaOn = fDetector -> GetBetaDetectorOn();
+ fBetaV2On = fDetector -> GetBetaDetectorV2On();
fVirtualOn = fDetector -> GetVirtualDetectorsOn();
- // get pointers to beta-detector phys volumes
+ // get pointers to beta-detector v1.0 phys volumes
if (fBetaOn) {
fSignalPV = fDetector->GetBetaDetector()->GetSignalPhysVol();
fVetoPV = fDetector->GetBetaDetector()->GetVetoPhysVol();
}
+
+ // get pointers to beta-detector v2.0 phys volumes
+ if (fBetaV2On) {
+ fFrontV2PV = fDetector->GetBetaDetectorV2()->GetFrontPhysVol();
+ fSignalV2PV = fDetector->GetBetaDetectorV2()->GetSignalPhysVol();
+ fVetoV2PV = fDetector->GetBetaDetectorV2()->GetVetoPhysVol();
+ }
// get pointers to virtual detectors
if (fVirtualOn) {
@@ -119,6 +129,13 @@ void VKRecorderRoot::RecordBeginOfEvent(const G4Event* anEvent)
fEDepSignal = 0;
fEDepVeto = 0;
+
+ fEDepV2Front = 0;
+ fEDepV2Signal = 0;
+ fEDepV2Veto = 0;
+
+ fPhotonsSignal = 0;
+ fPhotonsVeto = 0;
// clear Root branches
ClearRootBranches();
@@ -126,7 +143,7 @@ void VKRecorderRoot::RecordBeginOfEvent(const G4Event* anEvent)
// Record data about particle source:
// vertex
- fPrimVertex = anEvent -> GetPrimaryVertex(0);
+ fPrimVertex = anEvent -> GetPrimaryVertex(0);
// NB: if several vertices exist, this will only get the first one
// particle
@@ -156,14 +173,28 @@ void VKRecorderRoot::RecordEndOfEvent(const G4Event* anEvent)
VKUserEventInformation* eventInfo
= (VKUserEventInformation*)anEvent -> GetUserInformation();
- // deposited energy
- fEDepSignal = eventInfo->GetEDep(0);
- fEDepVeto = eventInfo->GetEDep(1);
+ int id = -1;
- // detected optical photons
- fPhotonsSignal = eventInfo->GetHitCount(0);
- fPhotonsVeto = eventInfo->GetHitCount(1);
+ // deposited energy
+ if (fBetaOn) {
+ id++;
+ fEDepSignal = eventInfo->GetEDep(id);
+ fPhotonsSignal = eventInfo->GetHitCount(id);
+ id++;
+ fEDepVeto = eventInfo->GetEDep(id);
+ fPhotonsVeto = eventInfo->GetHitCount(id);
+ }
+ // deposited energy (v2.0)
+ if (fBetaV2On) {
+ id++;
+ fEDepV2Front = eventInfo->GetEDep(id);
+ id++;
+ fEDepV2Signal = eventInfo->GetEDep(id);
+ id++;
+ fEDepV2Veto = eventInfo->GetEDep(id);
+ }
+
// fill ROOT tree
FillTree();
}
@@ -231,10 +262,19 @@ void VKRecorderRoot::CreateDetBranches()
{
if (fTFile == NULL) return;
- fDetTree -> Branch("EDepSignal", &_eDepSignal, "_eDepSignal/F");
- fDetTree -> Branch("EDepVeto" , &_eDepVeto, "_eDepVeto/F");
- fDetTree -> Branch("PhotonsSignal", &_photonsSignal, "_photonsSignal/i");
- fDetTree -> Branch("PhotonsVeto" , &_photonsVeto , "_photonsVeto/i");
+ if (fBetaOn) {
+ fDetTree -> Branch("EDepSignal", &_eDepSignal, "_eDepSignal/F");
+ fDetTree -> Branch("EDepVeto" , &_eDepVeto, "_eDepVeto/F");
+
+ fDetTree -> Branch("PhotonsSignal", &_photonsSignal, "_photonsSignal/i");
+ fDetTree -> Branch("PhotonsVeto" , &_photonsVeto , "_photonsVeto/i");
+ }
+
+ if (fBetaV2On) {
+ fDetTree -> Branch("EDepV2Front", &_eDepV2Front, "_eDepV2Front/F");
+ fDetTree -> Branch("EDepV2Signal", &_eDepV2Signal, "_eDepV2Signal/F");
+ fDetTree -> Branch("EDepV2Veto" , &_eDepV2Veto, "_eDepV2Veto/F");
+ }
}
@@ -272,18 +312,22 @@ void VKRecorderRoot::FillTree()
if (fTFile == NULL) return;
// zero suppression
- G4bool skip = fZeroSuppression && fEDepSignal==0 && fEDepVeto==0;
-
+ G4bool zero = fZeroSuppression && fEDepSignal==0 && fEDepVeto==0;
+ G4bool zeroV2 = fZeroSuppression && fEDepV2Front==0 && fEDepV2Signal==0 && fEDepV2Veto==0;
// beta detector
- _eDepSignal = fEDepSignal /keV;
- _eDepVeto = fEDepVeto /keV;
+ _eDepSignal = fEDepSignal /keV;
+ _eDepVeto = fEDepVeto /keV;
+
+ _eDepV2Front = fEDepV2Front /keV;
+ _eDepV2Signal = fEDepV2Signal /keV;
+ _eDepV2Veto = fEDepV2Veto /keV;
+
_photonsSignal = fPhotonsSignal;
_photonsVeto = fPhotonsVeto;
- if (fBetaOn && !skip) fDetTree -> Fill();
-
+ if ((fBetaOn && !zero) || (fBetaV2On && !zeroV2)) fDetTree -> Fill();
// virtual detectors
@@ -318,6 +362,7 @@ void VKRecorderRoot::FillTree()
void VKRecorderRoot::ClearRootBranches()
{
_eDepSignal = _eDepVeto = 0;
+ _eDepV2Front = _eDepV2Signal = _eDepV2Veto = 0;
_photonsSignal = _photonsVeto = 0;
_Xfp = _Yfp = _Zfp = _Efp = _THfp = _PHIfp = 0;
@@ -339,7 +384,7 @@ void VKRecorderRoot::OpenRootFile(G4String filename)
G4cout << " ---> Open file: " << fFilename << G4endl;
// Create trees
- if (fBetaOn) {
+ if (fBetaOn || fBetaV2On) {
G4cout << " ---> Create Detector tree" << G4endl;
fDetTree = new TTree("Detector", "Detector data");
CreateDetBranches();
@@ -358,7 +403,7 @@ void VKRecorderRoot::WriteToRootFileAndClose()
G4cout << " ---> Write to file: " << fFilename << G4endl;
- if (fBetaOn) fDetTree -> Write();
+ if (fBetaOn || fBetaV2On) fDetTree -> Write();
if (fVirtualOn) fTraTree -> Write();
CloseRootFile();
@@ -412,7 +457,7 @@ void VKRecorderRoot::WriteMetaDataToRootFile()
TParameter<double>* pNorm = new TParameter<double>("NORMALIZATION", normalization);
// Write meta data to detector tree
- if (fBetaOn) {
+ if (fBetaOn || fBetaV2On) {
fDetTree -> GetUserInfo() -> Add(pArea);
fDetTree -> GetUserInfo() -> Add(pCoverage);
fDetTree -> GetUserInfo() -> Add(pEvents);
diff --git a/src/VKSteppingAction.cpp b/src/VKSteppingAction.cpp
index 44ef59487d01aa20353ff476fd7ed9eafc859f76..4c41c1d8a81eeb382824004e6c64aaf445038229 100644
--- a/src/VKSteppingAction.cpp
+++ b/src/VKSteppingAction.cpp
@@ -143,7 +143,12 @@ void VKSteppingAction::UserSteppingAction(const G4Step * theStep)
// If user has selected to only allow primaries to go one step
// before being killed, kill track now
- if (fOneStepPrimaries && (thePrePV->GetName()=="signal" || thePrePV->GetName()=="veto")) {
+ if (fOneStepPrimaries && (thePrePV->GetName()=="signal"
+ || thePrePV->GetName()=="veto"
+ || thePrePV->GetName()=="frontV2"
+ || thePrePV->GetName()=="signalV2"
+ || thePrePV->GetName()=="vetoV2")) {
+
theTrack -> SetTrackStatus(fStopAndKill);
}
}
@@ -165,8 +170,10 @@ void VKSteppingAction::UserSteppingAction(const G4Step * theStep)
auto preName = thePrePV->GetName();
auto preCopyNo = thePrePV->GetCopyNo();
G4int SiPMScintId = 0;
- if (preName == "signal" || (preName == "SiPM" && preCopyNo == 0)) SiPMScintId = 0;
- else if (preName == "veto" || (preName == "SiPM" && preCopyNo == 1)) SiPMScintId = 1;
+ if (preName == "signal" || (preName == "SiPM" && preCopyNo == 0)) SiPMScintId = 0;
+ else if (preName == "veto" || (preName == "SiPM" && preCopyNo == 1)) SiPMScintId = 1;
+ else if (preName == "signalV2" || (preName == "SiPMV2" && preCopyNo == 0)) SiPMScintId = 0;
+ else if (preName == "vetoV2" || (preName == "SiPMV2" && preCopyNo == 1)) SiPMScintId = 1;
// If photon was absorbed by the absorption process
G4String procName = thePostPoint->GetProcessDefinedStep()->GetProcessName();
diff --git a/test/VKDetectorTest.cpp b/test/VKDetectorTest.cpp
index deb55666ecfe2b5cd88ef9e6312f797693eda725..1db904f579ef89ac89fc374c9886e2dfe36bc819 100644
--- a/test/VKDetectorTest.cpp
+++ b/test/VKDetectorTest.cpp
@@ -106,6 +106,13 @@ SUITE(VKDetectorTest) {
CHECK_EQUAL(true, logVolExists);
}
+ TEST_FIXTURE(SetupFixture, TestFrontV2Exists) {
+ UImanager->ApplyCommand("/VK/detector/volumes/betaDetectorV2 true");
+ G4bool verbose = false;
+ G4bool frontExists = physVolStore->GetVolume("frontV2", verbose) != nullptr;
+ CHECK_EQUAL(true, frontExists);
+ }
+
TEST_FIXTURE(SetupFixture, TestSignalSuccessfullyRemoved) {
UImanager->ApplyCommand("/VK/detector/volumes/betaDetector false");
G4bool verbose = false;
diff --git a/test/VKEnergyDepositionTest.cpp b/test/VKEnergyDepositionTest.cpp
index 0d16b70f48ecc966e3a8b1ef191ba619068a996a..754a3a5f17ea0a8e8131f64961c73634b7695e79 100644
--- a/test/VKEnergyDepositionTest.cpp
+++ b/test/VKEnergyDepositionTest.cpp
@@ -49,22 +49,16 @@ SUITE(VKEnergyDepositionTest) {
CHECK_EQUAL(1, 1);
}
- TEST_FIXTURE(SetupFixture, Test3MeVElectronInsideSignalDetector) {
+ TEST_FIXTURE(SetupFixture, Test1MeVElectronInsideSignalDetector) {
G4double edep=0;
UImanager->ApplyCommand("/gps/direction 0 0 1");
- UImanager->ApplyCommand("/gps/position 0 0 30.5 cm");
- // check energy deposition for 3 MeV electron
+ UImanager->ApplyCommand("/gps/position 0 0 32 cm");
+ // check energy deposition for 1 MeV electron
UImanager->ApplyCommand("/gps/particle e-");
- UImanager->ApplyCommand("/gps/energy 3 MeV");
+ UImanager->ApplyCommand("/gps/energy 1 MeV");
UImanager->ApplyCommand("/run/beamOn 1");
edep = recorder -> GetEDepSignal();
- CHECK_CLOSE(3000., edep/keV, 1e-3);
- // check for 1.5 MeV electron
- UImanager->ApplyCommand("/gps/particle e-");
- UImanager->ApplyCommand("/gps/energy 1.5 MeV");
- UImanager->ApplyCommand("/run/beamOn 1");
- edep = recorder -> GetEDepSignal();
- CHECK_CLOSE(1500., edep/keV, 1e-3);
+ CHECK_CLOSE(1000., edep/keV, 1);
}
TEST_FIXTURE(SetupFixture, Test1GeVVerticalMuon) {