diff --git a/.gitignore b/.gitignore index 521c740..193cfae 100644 --- a/.gitignore +++ b/.gitignore @@ -10,6 +10,7 @@ EventBuilder* *.png Mapper AnasenMS +Armory/anasenMS data/ data_proton/ diff --git a/Armory/Makefile b/Armory/Makefile index 6693ef3..14092ca 100644 --- a/Armory/Makefile +++ b/Armory/Makefile @@ -8,7 +8,11 @@ CC = g++ #COPTS = -fPIC -DLINUX -O2 -std=c++17 -lpthread COPTS = -fPIC -DLINUX -g -O0 -Wall -std=c++17 -lpthread -ROOTLIBS = `root-config --cflags --glibs` +ROOTCFLAGS := $(shell root-config --cflags) +ROOTLIBS := $(shell root-config --libs) + +CXX := g++ +CXXFLAGS := -O2 $(ROOTCFLAGS) ALL = Mapper EventBuilder#AnasenMS @@ -29,5 +33,8 @@ AnasenMS : constant.h Isotope.h ClassTransfer.h ClassSX3.h ClassPW.h ClassAnasen EventBuilder : EventBuilder.cpp ClassData.h fsuReader.h Hit.h @echo "--------- making EventBuilder" - $(CC) $(COPTS) -o EventBuilder EventBuilder.cpp $(ROOTLIBS) + $(CC) $(COPTS) -o EventBuilder EventBuilder.cpp $(ROOTLIBS) + +anasenMS: anasenMS.cpp + $(CXX) $(CXXFLAGS) anasenMS.cpp -o anasenMS $(ROOTLIBS) diff --git a/Armory/README_anasenMS.md b/Armory/README_anasenMS.md new file mode 100644 index 0000000..fd0543e --- /dev/null +++ b/Armory/README_anasenMS.md @@ -0,0 +1,63 @@ +# ANASEN Monte Carlo (anasenMS.cpp) + +## Overview + +`anasenMS.cpp` is a standalone Monte Carlo simulation for an ANASEN-style detector setup. It generates transfer reaction kinematics, propagates products to a wire chamber (PW) and a silicon array (SX3), reconstructs tracks, and writes output to a ROOT tree. + +## Requirements + +- ROOT (e.g. `root-config` for compile flags) +- C++ compiler (gcc/g++) +- Project includes: `ClassTransfer.h`, `ClassAnasen.h`, plus their dependent implementation files. + +## Build + +In `Armory` directory: + +```bash +g++ -O2 -o anasenMS anasenMS.cpp ClassTransfer.cpp ClassAnasen.cpp ... `root-config --cflags --libs` +``` + +(Adjust source file list based on actual project layout.) + +## Run + +```bash +./anasenMS [numEvents] +``` + +- `numEvents`: optional integer, default `1000000` +- Outputs: `SimAnasen1.root` containing `tree` + +## What the code does + +- Initializes reaction: `TransferReaction transfer` + - `SetA(24,12,0)` target + - `SetIncidentEnergyAngle(10,0,0)` beam energy and direction + - `Seta`, `Setb` reaction fragment indices +- Sets excitation lists: `ExAList`, `ExList` +- Vertex and resolution settings: + - `vertexX/Y/Z` ranges + - `sigmaSX3_W`, `sigmaSX3_L`, `sigmaPW_A`, `sigmaPW_C` +- Prepares ROOT output tree and branches for truth/reconstructed +- Loop over events: + - Sample excitation and CM direction + - `transfer.Event(thetaCM, phiCM)` outputs `TLorentzVector` products + - Compute lab angles/energies + - Random vertex inside target volume + - Run detector response: + - `pw->FindWireID(...)` + - `sx3->FindSX3Pos(...)` + - Read out wire hits and SX3 channel + depth + - Apply position smearing for SX3 + - Reconstruct track via `pw->CalTrack` and `pw->CalTrack2` + - Fill ROOT tree +- At end: write tree, close file, clean up + +## Notes + +- The code now has expanded inline comments explaining each stage. +- Important methods are from: + - `ClassTransfer` (`SetA`, `SetIncidentEnergyAngle`, `Seta`, `Setb`, `SetExA`, `SetExB`, `CalReactionConstant`, `Event`) + - `ClassAnasen` / `SX3` / `PW` (`FindWireID`, `FindSX3Pos`, `GetHitInfo`, `CalTrack`, `CalTrack2`, `GetTrackTheta`, `GetTrackPhi`, `GetZ0`, `GetHitPosWithSigma`, `GetID`, etc.) +- Optional: change excitation lists, vertex spread, and sigma values to mimic different beam/target conditions. diff --git a/Armory/anasenMS.cpp b/Armory/anasenMS.cpp index 162e64f..ec347be 100644 --- a/Armory/anasenMS.cpp +++ b/Armory/anasenMS.cpp @@ -1,52 +1,54 @@ -#include "TRandom.h" -#include "TFile.h" -#include "TTree.h" -#include "TH1.h" -#include "TH2.h" -#include "TStyle.h" -#include "TCanvas.h" -#include "TBenchmark.h" +#include "TRandom.h" // ROOT random number generators, gRandom +#include "TFile.h" // ROOT file I/O +#include "TTree.h" // ROOT tree storage +#include "TH1.h" // 1D histograms (not directly used here but common in analyzers) +#include "TH2.h" // 2D histograms +#include "TStyle.h" // ROOT plotting style controls +#include "TCanvas.h" // ROOT canvas drawing +#include "TBenchmark.h" // timing measurement -#include "ClassTransfer.h" -#include "ClassAnasen.h" +#include "ClassTransfer.h" // Reaction kinematics and MC event generation +#include "ClassAnasen.h" // ANASEN detector model classes (SX3, PW, etc.) -//======== Gerneate light particle based on reaction -// find out the CalTrack and the real track -// find out the Q-value uncertaintly +//======== Generate light particle based on reaction +// calculate real and reconstructed tracks and Q-value uncertainty int main(int argc, char **argv){ printf("=========================================\n"); printf("=== ANASEN Monte Carlo ===\n"); - printf("=========================================\n"); - + printf("=========================================\n"); + + // number of events can be overridden from command line int numEvent = 1000000; if( argc >= 2 ) numEvent = atoi(argv[1]); - //Reaction + // Reaction setup: projectile + target configuration, energy, and product IDs TransferReaction transfer; - transfer.SetA(24,12, 0); - transfer.SetIncidentEnergyAngle(10, 0, 0); - transfer.Seta( 4, 2); - transfer.Setb( 1, 1); + transfer.SetA(24,12, 0); // e.g., 24Mg (Z=12) with 0 excitation + transfer.SetIncidentEnergyAngle(10, 0, 0); // 10 MeV beam, 0 polar and azimuthal angle + transfer.Seta( 4, 2); // identify reaction product a in internal indexing + transfer.Setb( 1, 1); // identify reaction product b - //TODO add alpha source + // TODO add alpha source or alternative reaction channel selection + // Excited state lists (target and projectile/excited products) std::vector ExAList = {0}; std::vector ExList = {0, 1, 2}; + // define vertex position uniform distribution ranges (mm) double vertexXRange[2] = { -5, 5}; // mm - double vertexYRange[2] = { -5, 5}; - double vertexZRange[2] = { -100, 100}; + double vertexYRange[2] = { -5, 5}; + double vertexZRange[2] = { -100, 100}; - double sigmaSX3_W = -1; // mm, < 0 use mid-point - double sigmaSX3_L = 3; // mm, < 0 use mid-point - double sigmaPW_A = 0; // from 0 to 1. - double sigmaPW_C = 0; // from 0 to 1. - - //################################################### + // detector resolution / uncertainty parameters + double sigmaSX3_W = -1; // mm, if < 0 use mid-point (no spread in SX3 horizontal dimension) + double sigmaSX3_L = 3; // mm, vertical spread for SX3 + double sigmaPW_A = 0; // normalized anode uncertainty term (0-1) + double sigmaPW_C = 0; // normalized cathode uncertainty term (0-1) + // status printout printf("------------ Vertex :\n"); printf("X : %7.2f - %7.2f mm\n", vertexXRange[0], vertexXRange[1]); printf("Y : %7.2f - %7.2f mm\n", vertexYRange[0], vertexYRange[1]); @@ -57,20 +59,25 @@ int main(int argc, char **argv){ printf(" Anode : %.1f mm\n", sigmaPW_A); printf(" Cathode : %.1f mm\n", sigmaPW_C); printf(" num_eve : %d \n",numEvent); + + // calculates energy/momentum/kinematics constants for transfer reaction transfer.CalReactionConstant(); int nExA = ExAList.size(); int nEx = ExList.size(); - ANASEN * anasen = new ANASEN(); - SX3 * sx3 = anasen->GetSX3(); - PW * pw = anasen->GetPW(); + // create detector representation in memory + ANASEN * anasen = new ANASEN(); // top-level detector object + SX3 * sx3 = anasen->GetSX3(); // silicon array part + PW * pw = anasen->GetPW(); // proportional wire chamber part + // output file + tree TString saveFileName = "SimAnasen1.root"; printf("\e[32m#################################### building Tree in %s\e[0m\n", saveFileName.Data()); TFile * saveFile = new TFile(saveFileName, "recreate"); TTree * tree = new TTree("tree", "tree"); + // beam and CM variables saved in tree double KEA; tree->Branch("beamKEA", &KEA, "beamKEA/D"); @@ -78,6 +85,7 @@ int main(int argc, char **argv){ tree->Branch("thetaCM", &thetaCM, "thetaCM/D"); tree->Branch("phiCM", &phiCM, "phiCM/D"); + // outgoing particles in lab frame (light/heavy) double thetab, phib, Tb; double thetaB, phiB, TB; tree->Branch("thetab", &thetab, "thetab/D"); @@ -87,6 +95,7 @@ int main(int argc, char **argv){ tree->Branch("phiB", &phiB, "phiB/D"); tree->Branch("TB", &TB, "TB/D"); + // excitation state identifiers int ExAID; double ExA; tree->Branch("ExAID", &ExAID, "ExAID/I"); @@ -97,24 +106,29 @@ int main(int argc, char **argv){ tree->Branch("ExID", &ExID, "ExID/I"); tree->Branch("Ex", &Ex, "Ex/D"); + // true vertex position in target volume double vertexX, vertexY, vertexZ; tree->Branch("vX", &vertexX, "VertexX/D"); tree->Branch("vY", &vertexY, "VertexY/D"); tree->Branch("vZ", &vertexZ, "VertexZ/D"); + // reconstructed SX3 hit position double sx3X, sx3Y, sx3Z; tree->Branch("sx3X", &sx3X, "sx3X/D"); tree->Branch("sx3Y", &sx3Y, "sx3Y/D"); tree->Branch("sx3Z", &sx3Z, "sx3Z/D"); + // PW nearest and next nearest wires int anodeID[2], cathodeID[2]; tree->Branch("aID", anodeID, "anodeID/I"); tree->Branch("cID", cathodeID, "cathodeID/I"); + // distances to nearest wires double anodeDist[2], cathodeDist[2]; tree->Branch("aDist", anodeDist, "anodeDist/D"); tree->Branch("cDist", cathodeDist, "cathodeDist/D"); + // SX3 channel assignment and Z fraction (depth) information int sx3ID, sx3Up, sx3Dn, sx3Bk; double sx3ZFrac; tree->Branch("sx3ID", &sx3ID, "sx3ID/I"); @@ -123,6 +137,7 @@ int main(int argc, char **argv){ tree->Branch("sx3Bk", &sx3Bk, "sx3Bk/I"); tree->Branch("sx3ZFrac", &sx3ZFrac, "sx3ZFrac/D"); + // reconstructed angles from PW track fit, method 1 and 2 double reTheta, rePhi; tree->Branch("reTheta", &reTheta, "reconstucted_theta/D"); tree->Branch("rePhi", &rePhi, "reconstucted_phi/D"); @@ -131,20 +146,22 @@ int main(int argc, char **argv){ tree->Branch("reTheta1", &reTheta1, "reconstucted_theta1/D"); tree->Branch("rePhi1", &rePhi1, "reconstucted_phi1/D"); + // reconstructed vertex Z from PW fit double z0; tree->Branch("z0", &z0, "reconstucted_Z/D"); //========timer TBenchmark clock; - bool shown ; + bool shown ; clock.Reset(); clock.Start("timer"); shown = false; - //================================= Calculate event + //================================= Calculate event loop for( int i = 0; i < numEvent ; i++){ + // randomly sample target/projectile excitations ExAID = gRandom->Integer(nExA); ExA = ExAList[ExAID]; transfer.SetExA(ExA); @@ -153,15 +170,17 @@ int main(int argc, char **argv){ Ex = ExList[ExID]; transfer.SetExB(Ex); + // recalc kinematic constants for chosen states transfer.CalReactionConstant(); - thetaCM = TMath::ACos(2 * gRandom->Rndm() - 1) ; + // isotropic CM direction + thetaCM = TMath::ACos(2 * gRandom->Rndm() - 1) ; phiCM = (gRandom->Rndm() - 0.5) * TMath::TwoPi(); - //==== Calculate reaction - TLorentzVector * output = transfer.Event(thetaCM, phiCM); - TLorentzVector Pb = output[2]; - TLorentzVector PB = output[3]; + //==== Calculate reaction kinematics in lab frame + TLorentzVector * output = transfer.Event(thetaCM, phiCM); // returns array of outputs + TLorentzVector Pb = output[2]; // light particle or product A + TLorentzVector PB = output[3]; // heavy particle or product B thetab = Pb.Theta() * TMath::RadToDeg(); thetaB = PB.Theta() * TMath::RadToDeg(); @@ -172,22 +191,25 @@ int main(int argc, char **argv){ phib = Pb.Phi() * TMath::RadToDeg(); phiB = PB.Phi() * TMath::RadToDeg(); + // vertex position in target volume vertexX = (vertexXRange[1]- vertexXRange[0])*gRandom->Rndm() + vertexXRange[0]; vertexY = (vertexYRange[1]- vertexYRange[0])*gRandom->Rndm() + vertexYRange[0]; vertexZ = (vertexZRange[1]- vertexZRange[0])*gRandom->Rndm() + vertexZRange[0]; TVector3 vertex(vertexX, vertexY, vertexZ); + // set direction vector from lab angle TVector3 dir(1, 0, 0); dir.SetTheta(thetab * TMath::DegToRad()); dir.SetPhi(phib * TMath::DegToRad()); - + // run detector response models for PW and SX3 pw->FindWireID(vertex, dir, false); - sx3->FindSX3Pos(vertex, dir, false); + sx3->FindSX3Pos(vertex, dir, false); PWHitInfo hitInfo = pw->GetHitInfo(); + // store nearest/next closest wire data anodeID[0] = hitInfo.nearestWire.first; cathodeID[0] = hitInfo.nearestWire.second; anodeID[1] = hitInfo.nextNearestWire.first; @@ -198,6 +220,7 @@ int main(int argc, char **argv){ anodeDist[1] = hitInfo.nextNearestDist.first; cathodeDist[1] = hitInfo.nextNearestDist.second; + // SX3 hit channel info and depth fraction sx3ID = sx3->GetID(); if( sx3ID >= 0 ){ sx3Up = sx3->GetChUp(); @@ -205,51 +228,48 @@ int main(int argc, char **argv){ sx3Bk = sx3->GetChBk(); sx3ZFrac = sx3->GetZFrac(); - //Introduce uncertaity - // TVector3 hitPos = sx3->GetHitPos(); + // apply intrinsic detector resolution to true SX3 hit position + // for no smearing comment out and use GetHitPos(); TVector3 hitPos = sx3->GetHitPosWithSigma(sigmaSX3_W, sigmaSX3_L); sx3X = hitPos.X(); sx3Y = hitPos.Y(); sx3Z = hitPos.Z(); - + + // reconstruct track from PW readings + SX3 hit pw->CalTrack(hitPos, anodeID[0], cathodeID[0], false); reTheta = pw->GetTrackTheta() * TMath::RadToDeg(); - rePhi = pw->GetTrackPhi() * TMath::RadToDeg(); + rePhi = pw->GetTrackPhi() * TMath::RadToDeg(); + // alternative track algorithm with uncertainty parameters pw->CalTrack2(hitPos, hitInfo, sigmaPW_A, sigmaPW_C, false); reTheta1 = pw->GetTrackTheta() * TMath::RadToDeg(); - rePhi1 = pw->GetTrackPhi() * TMath::RadToDeg(); + rePhi1 = pw->GetTrackPhi() * TMath::RadToDeg(); z0 = pw->GetZ0(); }else{ + // no valid SX3 hit: mark clearly invalid sx3Up = -1; sx3Dn = -1; sx3Bk = -1; sx3ZFrac = TMath::QuietNaN(); - + sx3X = TMath::QuietNaN(); sx3Y = TMath::QuietNaN(); sx3Z = TMath::QuietNaN(); - // for( int i = 0; i < 12; i++){ - // sx3Index[i] = -1; - // } - - reTheta = TMath::QuietNaN(); - rePhi = TMath::QuietNaN(); - + reTheta = TMath::QuietNaN(); + rePhi = TMath::QuietNaN(); reTheta1 = TMath::QuietNaN(); - rePhi1 = TMath::QuietNaN(); - - z0 = TMath::QuietNaN(); - + rePhi1 = TMath::QuietNaN(); + z0 = TMath::QuietNaN(); } tree->Fill(); - //#################################################################### Timer + //#################################################################### Timer + // measure elapsed real time and print progress roughly every 10 sec clock.Stop("timer"); Double_t time = clock.GetRealTime("timer"); clock.Start("timer"); @@ -259,7 +279,7 @@ int main(int argc, char **argv){ printf( "%10d[%2d%%]| %8.2f sec | expect: %5.1f min \n", i, TMath::Nint((i+1)*100./numEvent), time , numEvent*time/(i+1)/60); shown = 1; } - }else{ + } else { if (fmod(time, 10) > 9 ){ shown = 0; } @@ -267,6 +287,7 @@ int main(int argc, char **argv){ } + // write results to ROOT file and close tree->Write(); int count = tree->GetEntries(); saveFile->Close(); diff --git a/DetectorConstruction.cc b/DetectorConstruction.cc new file mode 100644 index 0000000..42be93b --- /dev/null +++ b/DetectorConstruction.cc @@ -0,0 +1,131 @@ +#include "DetectorConstruction.hh" + +#include "G4RunManager.hh" +#include "G4NistManager.hh" +#include "G4Box.hh" +#include "G4Tubs.hh" +#include "G4LogicalVolume.hh" +#include "G4PVPlacement.hh" +#include "G4SystemOfUnits.hh" +#include "G4RotationMatrix.hh" +#include "G4ThreeVector.hh" + +DetectorConstruction::DetectorConstruction() +: G4VUserDetectorConstruction(), + fScoringVolume(nullptr) +{} + +DetectorConstruction::~DetectorConstruction() +{} + +G4VPhysicalVolume* DetectorConstruction::Construct() +{ + // Get nist material manager + G4NistManager* nist = G4NistManager::Instance(); + + // Materials + G4Material* vacuum = nist->FindOrBuildMaterial("G4_Galactic"); + G4Material* al = nist->FindOrBuildMaterial("G4_Al"); + + // World + G4double worldx = 200.*mm; + G4double worldy = 200.*mm; + G4double worldz = 200.*mm; + G4Box* solidWorld = new G4Box("World", worldx, worldy, worldz); + G4LogicalVolume* logicWorld = new G4LogicalVolume(solidWorld, vacuum, "World"); + G4VPhysicalVolume* physWorld = new G4PVPlacement(0, G4ThreeVector(), logicWorld, "World", 0, false, 0, true); + + // Axes (optional, for visualization) + G4Tubs* solidAxis = new G4Tubs("Axis", 0, 0.1*mm, 5.*mm, 0, 360*deg); + G4LogicalVolume* logicAxisX = new G4LogicalVolume(solidAxis, al, "AxisX"); + G4LogicalVolume* logicAxisY = new G4LogicalVolume(solidAxis, al, "AxisY"); + G4LogicalVolume* logicAxisZ = new G4LogicalVolume(solidAxis, al, "AxisZ"); + + G4RotationMatrix* rotX = new G4RotationMatrix(); + rotX->rotateY(90*deg); + rotX->rotateZ(90*deg); + new G4PVPlacement(rotX, G4ThreeVector(5*mm, 0, 0), logicAxisX, "AxisX", logicWorld, false, 0); + + G4RotationMatrix* rotY = new G4RotationMatrix(); + rotY->rotateX(90*deg); + new G4PVPlacement(rotY, G4ThreeVector(0, 5*mm, 0), logicAxisY, "AxisY", logicWorld, false, 0); + + new G4PVPlacement(0, G4ThreeVector(0, 0, 5*mm), logicAxisZ, "AxisZ", logicWorld, false, 0); + + // ANASEN geometry + const int nWire = 24; + const int wireShift = 3; + const G4double zLen = 350.*mm; + const G4double radiusA = 38.*mm; + const G4double radiusC = 43.*mm; + + G4double dAngle = wireShift * 2 * CLHEP::pi / nWire; + G4double radiusAnew = radiusA * cos(dAngle / 2.); + G4double wireALength = sqrt(zLen*zLen + pow(2 * radiusA * sin(dAngle/2.), 2)); + G4double wireATheta = atan2(2 * radiusA * sin(dAngle / 2.), zLen); + + G4Tubs* solidPC_A = new G4Tubs("PC_A", 0, 0.01*mm, wireALength/2., 0, 360*deg); + G4LogicalVolume* logicPC_A = new G4LogicalVolume(solidPC_A, al, "PC_A"); + + for(int i = 0; i < nWire; i++){ + G4double phi = 2 * CLHEP::pi / nWire * i + dAngle / 2.; + G4ThreeVector pos(radiusAnew * cos(phi), radiusAnew * sin(phi), 0); + G4RotationMatrix* rot = new G4RotationMatrix(); + rot->rotateZ(360./nWire * (i + wireShift/2.) * deg); + rot->rotateY(wireATheta); + new G4PVPlacement(rot, pos, logicPC_A, "PC_A", logicWorld, false, i); + } + + G4double radiusCnew = radiusC * cos(dAngle / 2.); + G4double wireCLength = sqrt(zLen*zLen + pow(2 * radiusC * sin(dAngle/2.), 2)); + G4double wireCTheta = atan2(2 * radiusC * sin(dAngle / 2.), zLen); + + G4Tubs* solidPC_C = new G4Tubs("PC_C", 0, 0.01*mm, wireCLength/2., 0, 360*deg); + G4LogicalVolume* logicPC_C = new G4LogicalVolume(solidPC_C, al, "PC_C"); + + for(int i = 0; i < nWire; i++){ + G4double phi = 2 * CLHEP::pi / nWire * i - dAngle/2.; + G4ThreeVector pos(radiusCnew * cos(phi), radiusCnew * sin(phi), 0); + G4RotationMatrix* rot = new G4RotationMatrix(); + rot->rotateZ(360./nWire * (i - wireShift/2.) * deg); + rot->rotateY(-wireCTheta); + new G4PVPlacement(rot, pos, logicPC_C, "PC_C", logicWorld, false, i); + } + + const int nSX3 = 12; + const G4double sx3Radius = 88.*mm; + const G4double sx3Width = 40.*mm; + const G4double sx3Length = 75.*mm; + const G4double sx3Gap = 5.*mm; + + G4Box* solidSX3 = new G4Box("SX3", 0.1*mm, sx3Width/2., sx3Length/2.); + G4LogicalVolume* logicSX3 = new G4LogicalVolume(solidSX3, al, "SX3"); + fScoringVolume = logicSX3; + + for(int i = 0; i < nSX3; i++){ + G4double phi = 2 * CLHEP::pi / nSX3 * (i + 0.5); + G4ThreeVector pos1(sx3Radius * cos(phi), sx3Radius * sin(phi), sx3Length/2. + sx3Gap); + G4RotationMatrix* rot1 = new G4RotationMatrix(); + rot1->rotateZ(360./nSX3 * (i + 0.5) * deg); + new G4PVPlacement(rot1, pos1, logicSX3, "SX3_front", logicWorld, false, 2*i); + + G4ThreeVector pos2(sx3Radius * cos(phi), sx3Radius * sin(phi), -sx3Length/2. - sx3Gap); + G4RotationMatrix* rot2 = new G4RotationMatrix(); + rot2->rotateZ(360./nSX3 * (i + 0.5) * deg); + new G4PVPlacement(rot2, pos2, logicSX3, "SX3_back", logicWorld, false, 2*i+1); + } + + const G4double qqqR1 = 50.*mm; + const G4double qqqR2 = 100.*mm; + G4Tubs* solidQQQ = new G4Tubs("QQQ", qqqR1, qqqR2, 0.5*mm, 5*deg, 85*deg); + G4LogicalVolume* logicQQQ = new G4LogicalVolume(solidQQQ, al, "QQQ"); + + for(int i = 0; i < 4; i++){ + G4ThreeVector pos(0, 0, 100.*mm); + G4RotationMatrix* rot = new G4RotationMatrix(); + rot->rotateZ(360./4 * i * deg); + new G4PVPlacement(rot, pos, logicQQQ, "QQQ", logicWorld, false, i); + } + + return physWorld; +} \ No newline at end of file diff --git a/DetectorConstruction.hh b/DetectorConstruction.hh new file mode 100644 index 0000000..bb64854 --- /dev/null +++ b/DetectorConstruction.hh @@ -0,0 +1,20 @@ +#ifndef DetectorConstruction_h +#define DetectorConstruction_h 1 + +#include "G4VUserDetectorConstruction.hh" + +class DetectorConstruction : public G4VUserDetectorConstruction +{ + public: + DetectorConstruction(); + virtual ~DetectorConstruction(); + + virtual G4VPhysicalVolume* Construct(); + + G4LogicalVolume* GetScoringVolume() const { return fScoringVolume; } + + private: + G4LogicalVolume* fScoringVolume; +}; + +#endif \ No newline at end of file