From 1d5492b1cbba4b8fe23b3f94c4791938c705a9ab Mon Sep 17 00:00:00 2001 From: AnalitikGnid Date: Wed, 28 May 2025 12:42:54 +0200 Subject: [PATCH 01/25] prestrawdetector added --- CMakeLists.txt | 2 + prestrawdetector/.prestrawdetector.cxx.swo | 0 prestrawdetector/CMakeLists.txt | 30 +++ prestrawdetector/prestrawdetector.cxx | 219 +++++++++++++++++++++ prestrawdetector/prestrawdetector.h | 75 +++++++ prestrawdetector/prestrawdetectorLinkDef.h | 10 + prestrawdetector/prestrawdetectorPoint.cxx | 42 ++++ prestrawdetector/prestrawdetectorPoint.h | 48 +++++ 8 files changed, 426 insertions(+) create mode 100644 prestrawdetector/.prestrawdetector.cxx.swo create mode 100644 prestrawdetector/CMakeLists.txt create mode 100644 prestrawdetector/prestrawdetector.cxx create mode 100644 prestrawdetector/prestrawdetector.h create mode 100644 prestrawdetector/prestrawdetectorLinkDef.h create mode 100644 prestrawdetector/prestrawdetectorPoint.cxx create mode 100644 prestrawdetector/prestrawdetectorPoint.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 43b666beab..d1e4136d87 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -220,11 +220,13 @@ add_subdirectory(veto) add_subdirectory(TimeDet) add_subdirectory(UpstreamTagger) add_subdirectory(strawtubes) +add_subdirectory (prestrawdetector) add_subdirectory(shipgen) add_subdirectory(field) add_subdirectory(muonShieldOptimization) add_subdirectory(tests) + add_custom_target( geometry.link ALL COMMAND ${CMAKE_COMMAND} -E create_symlink ${CMAKE_SOURCE_DIR}/geometry diff --git a/prestrawdetector/.prestrawdetector.cxx.swo b/prestrawdetector/.prestrawdetector.cxx.swo new file mode 100644 index 0000000000..e69de29bb2 diff --git a/prestrawdetector/CMakeLists.txt b/prestrawdetector/CMakeLists.txt new file mode 100644 index 0000000000..078844f895 --- /dev/null +++ b/prestrawdetector/CMakeLists.txt @@ -0,0 +1,30 @@ +set(INCLUDE_DIRECTORIES +${CMAKE_SOURCE_DIR}/shipdata +${CMAKE_SOURCE_DIR}/prestrawdetector +${genfit2_INCDIR} +) + +include_directories(${INCLUDE_DIRECTORIES} ${VMC_INCLUDE_DIRS} ${FAIRROOT_INCLUDE_DIR}) +include_directories(SYSTEM ${SYSTEM_INCLUDE_DIRECTORIES}) + +set(LINK_DIRECTORIES +${ROOT_LIBRARY_DIR} +${FAIRROOT_LIBRARY_DIR} +${genfit2_LIBDIR} +) + +link_directories( ${LINK_DIRECTORIES}) + +set(SRCS +#Put here your sourcefiles +prestrawdetector.cxx +prestrawdetectorPoint.cxx +) + +Set(LINKDEF prestrawdetectorLinkDef.h) +Set(LIBRARY_NAME prestrawdetector) +Set(DEPENDENCIES + Base ShipData FairLogger::FairLogger +) + +GENERATE_LIBRARY() diff --git a/prestrawdetector/prestrawdetector.cxx b/prestrawdetector/prestrawdetector.cxx new file mode 100644 index 0000000000..1e03832405 --- /dev/null +++ b/prestrawdetector/prestrawdetector.cxx @@ -0,0 +1,219 @@ +#include "prestrawdetector.h" +#include "prestrawdetectorPoint.h" + +#include "FairVolume.h" +#include "FairGeoVolume.h" +#include "FairGeoNode.h" +#include "FairRootManager.h" +#include "FairGeoLoader.h" +#include "FairGeoInterface.h" +#include "FairGeoMedia.h" +#include "FairGeoBuilder.h" +#include "FairRun.h" +#include "FairRuntimeDb.h" +#include "ShipDetectorList.h" +#include "ShipStack.h" + +#include "TClonesArray.h" +#include "TVirtualMC.h" +#include "TGeoManager.h" +#include "TGeoBBox.h" +#include "TGeoCompositeShape.h" +#include "TGeoTube.h" +#include "TGeoMaterial.h" +#include "TGeoMedium.h" +#include "TMath.h" +#include "TParticle.h" +#include "TVector3.h" + + +prestrawdetector::prestrawdetector() + : FairDetector("prestrawdetector", kTRUE) + , fTrackID(-1) + , fVolumeID(-1) + , fPos() + , fMom() + , fTime(-1.) + , fLength(-1.) + , fELoss(-1) + , fMedium("air") + , fprestrawdetectorPointCollection(new TClonesArray("prestrawdetectorPoint")) +{} + +prestrawdetector::prestrawdetector(const char* name, Bool_t active) + : FairDetector(name, active) + , fTrackID(-1) + , fVolumeID(-1) + , fPos() + , fMom() + , fTime(-1.) + , fLength(-1.) + , fELoss(-1) + , fMedium("air") + , fprestrawdetectorPointCollection(new TClonesArray("prestrawdetectorPoint")) +{} + +prestrawdetector::~prestrawdetector() { + if (fprestrawdetectorPointCollection) { + fprestrawdetectorPointCollection->Delete(); + delete fprestrawdetectorPointCollection; + } +} + +void prestrawdetector::Initialize() +{ + FairDetector::Initialize(); +// FairRuntimeDb* rtdb= FairRun::Instance()->GetRuntimeDb(); +// vetoGeoPar* par=(vetoGeoPar*)(rtdb->getContainer("vetoGeoPar")); +} + +// ----- Private method InitMedium +Int_t prestrawdetector::InitMedium(const char* name) +{ + static FairGeoLoader *geoLoad=FairGeoLoader::Instance(); + static FairGeoInterface *geoFace=geoLoad->getGeoInterface(); + static FairGeoMedia *media=geoFace->getMedia(); + static FairGeoBuilder *geoBuild=geoLoad->getGeoBuilder(); + + FairGeoMedium *ShipMedium=media->getMedium(name); + + if (!ShipMedium) + { + Fatal("InitMedium","Material %s not defined in media file.", name); + return -1111; + } + TGeoMedium* medium=gGeoManager->GetMedium(name); + if (medium!=NULL) + return ShipMedium->getMediumIndex(); + + return geoBuild->createMedium(ShipMedium); +} + +Bool_t prestrawdetector::ProcessHits(FairVolume* vol) +{ + + std::cout << "In PH\n"; + /** This method is called from the MC stepping */ + //Set parameters at entrance of volume. Reset ELoss. + if ( gMC->IsTrackEntering() ) { + fELoss = 0.; + fTime = gMC->TrackTime() * 1.0e09; + fLength = gMC->TrackLength(); + gMC->TrackPosition(fPos); + gMC->TrackMomentum(fMom); + } + // Sum energy loss for all steps in the active volume + fELoss += gMC->Edep(); + //std::cout << gMC->IsTrackExiting()<IsTrackStop()<IsTrackDisappeared(); + // Create strawtubesPoint at exit of active volume + if ( gMC->IsTrackEntering() ) { + //if (fELoss == 0. ) { return kFALSE; } + std::cout << "Entering track was found\n"; + TParticle* p=gMC->GetStack()->GetCurrentTrack(); + Int_t pdgCode = p->GetPdgCode(); + fTrackID = gMC->GetStack()->GetCurrentTrackNumber(); + Int_t det_uniqueId; + gMC->CurrentVolID(det_uniqueId); + /*if (fVolumeID == det_uniqueId) { + //std::cout << pdgCode<< " same volume again ? "<< straw_uniqueId << " exit:" << gMC->IsTrackExiting() << " stop:" << gMC->IsTrackStop() << " disappeared:" << gMC->IsTrackDisappeared()<< std::endl; + return kTRUE; }*/ + fVolumeID = det_uniqueId; + // # d = |pq . u x v|/|u x v| + TLorentzVector Pos; + gMC->TrackPosition(Pos); + Double_t xmean = (fPos.X()+Pos.X())/2. ; + Double_t ymean = (fPos.Y()+Pos.Y())/2. ; + Double_t zmean = (fPos.Z()+Pos.Z())/2. ; + Double_t deltaTrackLength = gMC->TrackLength() - fLength; + AddHit(fTrackID, det_uniqueId, TVector3(xmean, ymean, zmean), + TVector3(fMom.Px(), fMom.Py(), fMom.Pz()), fTime, deltaTrackLength, + fELoss,pdgCode); + // Increment number of strawtubes det points in TParticle + /*ShipStack* stack = (ShipStack*) gMC->GetStack(); + stack->AddPoint(kStraw);*/ + } + return kTRUE; +} + +prestrawdetectorPoint* prestrawdetector::AddHit(Int_t trackID, Int_t detID, + TVector3 pos, TVector3 mom, + Double_t time, Double_t length, + Double_t eLoss, Int_t pdgCode) +{ + TClonesArray& clref = *fprestrawdetectorPointCollection; + Int_t size = clref.GetEntriesFast(); + std::cout << "adding hit detid " << detID <Register("prestrawdetectorPoint", "prestrawdetector", + fprestrawdetectorPointCollection, kTRUE); + std::cout << "reg2"; + +} + +TClonesArray* prestrawdetector::GetCollection(Int_t iColl) const +{ + if (iColl == 0) { return fprestrawdetectorPointCollection; } + else { return NULL; } +} + +void prestrawdetector::Reset() +{ + fprestrawdetectorPointCollection->Clear(); +} + +void prestrawdetector::SetZ(Double_t z0) +{ + fPSDz = z0; //! z-position of veto station +} + +void prestrawdetector::ConstructGeometry() +{ + TGeoVolume *top = gGeoManager->GetTopVolume(); + InitMedium("air"); + TGeoMedium *air = gGeoManager->GetMedium("air"); + InitMedium("ShipSens"); + TGeoMedium *Se = gGeoManager->GetMedium("ShipSens"); + InitMedium("aluminium"); + TGeoMedium *Al = gGeoManager->GetMedium("aluminium"); + InitMedium("mylar"); + TGeoMedium *mylar = gGeoManager->GetMedium("mylar"); + InitMedium("STTmix9010_2bar"); + TGeoMedium *sttmix9010_2bar = gGeoManager->GetMedium("STTmix9010_2bar"); + InitMedium("tungsten"); + TGeoMedium *tungsten = gGeoManager->GetMedium("tungsten"); + + + gGeoManager->SetVisLevel(4); + gGeoManager->SetTopVisible(); + Double_t eps=0.0001; + + TGeoBBox* myBox = new TGeoBBox("myBox", 300, 300, 1); + TGeoVolume* myVol = new TGeoVolume("myVol", myBox, air); + myVol->SetLineColor(kYellow); + top->AddNode(myVol, 0, new TGeoTranslation(0, 0, fPSDz)); + + AddSensitiveVolume(myVol); + + // Position it somewhere + // top->AddNode(myVol, 1, new TGeoTranslation(0, 0, 500)); +} + +void prestrawdetector::EndOfEvent() +{ + fprestrawdetectorPointCollection->Clear(); +} + + + diff --git a/prestrawdetector/prestrawdetector.h b/prestrawdetector/prestrawdetector.h new file mode 100644 index 0000000000..9df4a3e6d8 --- /dev/null +++ b/prestrawdetector/prestrawdetector.h @@ -0,0 +1,75 @@ +#ifndef PRESTRAWDETECOR_H +#define PRESTRAWDETECOR_H + +#include "FairDetector.h" +#include "FairVolume.h" +#include "TVector3.h" +#include "TClonesArray.h" +#include "TLorentzVector.h" + +class prestrawdetectorPoint; +class FairVolume; +class TClonesArray; + +class prestrawdetector : public FairDetector +{ + public: + prestrawdetector(const char* name, Bool_t active); + prestrawdetector(); + virtual ~prestrawdetector(); + + /** Initialization of the detector is done here */ + virtual void Initialize(); //done + + virtual Bool_t ProcessHits( FairVolume* v=0);//done + + prestrawdetectorPoint* AddHit(Int_t trackID, Int_t detID, + TVector3 pos, TVector3 mom, + Double_t time, Double_t length, + Double_t eLoss, Int_t pdgCode); + + /** Registers the produced collections in FAIRRootManager. */ + virtual void Register(); + + /** Gets the produced collections */ + virtual TClonesArray* GetCollection(Int_t iColl) const ; //done + + /** has to be called after each event to reset the containers */ + virtual void Reset(); + + void SetZ(Double_t z0); + + void ConstructGeometry(); + + virtual void CopyClones( TClonesArray* cl1, TClonesArray* cl2 , + Int_t offset) {;} + virtual void SetSpecialPhysicsCuts() {;} + virtual void EndOfEvent(); + virtual void FinishPrimary() {;} + virtual void FinishRun() {;} + virtual void BeginPrimary() {;} + virtual void PostTrack() {;} + virtual void PreTrack() {;} + virtual void BeginEvent() {;} + + + private: + Double_t fPSDz; + std::string fMedium; + Int_t fTrackID; //! track index + Int_t fVolumeID; //! volume id + TLorentzVector fPos; //! position at entrance + TLorentzVector fMom; //! momentum at entrance + Double_t fTime; //! time + Double_t fLength; //! length + Double_t fELoss; + TClonesArray* fprestrawdetectorPointCollection; // collection of hit points + + prestrawdetector(const prestrawdetector&); + prestrawdetector& operator=(const prestrawdetector&); + + Int_t InitMedium(const char* name); + ClassDef(prestrawdetector,2) +}; + +#endif diff --git a/prestrawdetector/prestrawdetectorLinkDef.h b/prestrawdetector/prestrawdetectorLinkDef.h new file mode 100644 index 0000000000..997d4e11c0 --- /dev/null +++ b/prestrawdetector/prestrawdetectorLinkDef.h @@ -0,0 +1,10 @@ +#ifdef __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class prestrawdetector+; +#pragma link C++ class prestrawdetectorPoint+; + +#endif diff --git a/prestrawdetector/prestrawdetectorPoint.cxx b/prestrawdetector/prestrawdetectorPoint.cxx new file mode 100644 index 0000000000..db609a81e5 --- /dev/null +++ b/prestrawdetector/prestrawdetectorPoint.cxx @@ -0,0 +1,42 @@ +#include "prestrawdetectorPoint.h" + +#include +#include +using std::cout; +using std::endl; + + +// ----- Default constructor ------------------------------------------- +prestrawdetectorPoint::prestrawdetectorPoint() + : FairMCPoint() +{ +} +// ------------------------------------------------------------------------- + +// ----- Standard constructor ------------------------------------------ +prestrawdetectorPoint::prestrawdetectorPoint(Int_t trackID, Int_t detID, + TVector3 pos, TVector3 mom, + Double_t tof, Double_t length, + Double_t eLoss,Int_t pdgcode) + : FairMCPoint(trackID, detID, pos, mom, tof, length, eLoss), fPdgCode(pdgcode) +{ +} +// ------------------------------------------------------------------------- + +// ----- Destructor ---------------------------------------------------- +prestrawdetectorPoint::~prestrawdetectorPoint() { } +// ------------------------------------------------------------------------- + +// ----- Public method Print ------------------------------------------- +void prestrawdetectorPoint::Print() const +{ + cout << "-I- prestrawdetectorPoint: prestrawdetector point for track " << fTrackID + << " in detector " << fDetectorID << endl; + cout << " Position (" << fX << ", " << fY << ", " << fZ + << ") cm" << endl; + cout << " Momentum (" << fPx << ", " << fPy << ", " << fPz + << ") GeV" << endl; + cout << " Time " << fTime << " ns, Length " << fLength << "pdgid" << fPdgCode + << " cm, Energy loss " << fELoss*1.0e06 << " keV" << endl; +} +// ------------------------------------------------------------------------- diff --git a/prestrawdetector/prestrawdetectorPoint.h b/prestrawdetector/prestrawdetectorPoint.h new file mode 100644 index 0000000000..1bb1d988a1 --- /dev/null +++ b/prestrawdetector/prestrawdetectorPoint.h @@ -0,0 +1,48 @@ +#ifndef PRESTRAWDETECTORPOINT_H +#define PRESTRAWDETECTORPOINT_H 1 + + +#include "FairMCPoint.h" + +#include "TObject.h" +#include "TVector3.h" + +class prestrawdetectorPoint : public FairMCPoint +{ + public: + + /** Default constructor **/ + prestrawdetectorPoint(); + + /** Constructor with arguments + *@param trackID Index of MCTrack + *@param detID Detector ID + *@param pos Ccoordinates at entrance to active volume [cm] + *@param mom Momentum of track at entrance [GeV] + *@param tof Time since event start [ns] + *@param length Track length since creation [cm] + *@param eLoss Energy deposit [GeV] + **/ + prestrawdetectorPoint(Int_t trackID, Int_t detID, TVector3 pos, TVector3 mom, + Double_t tof, Double_t length, Double_t eLoss, Int_t pdgcode); + + /** Destructor **/ + virtual ~prestrawdetectorPoint(); + + /** Output to screen **/ + virtual void Print() const; + Int_t PdgCode() const {return fPdgCode;} + + + private: + /** Copy constructor **/ + prestrawdetectorPoint(const prestrawdetectorPoint& point); + prestrawdetectorPoint operator=(const prestrawdetectorPoint& point); + + Int_t fPdgCode; + ClassDef(prestrawdetectorPoint,2); + + +}; + +#endif From 2977ceee4363449477c27b2aa86714b6128be6fb Mon Sep 17 00:00:00 2001 From: AnalitikGnid Date: Wed, 28 May 2025 12:44:28 +0200 Subject: [PATCH 02/25] prestrawdetector added --- prestrawdetector/.prestrawdetector.cxx.swo | 0 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 prestrawdetector/.prestrawdetector.cxx.swo diff --git a/prestrawdetector/.prestrawdetector.cxx.swo b/prestrawdetector/.prestrawdetector.cxx.swo deleted file mode 100644 index e69de29bb2..0000000000 From 815c18fdb792d32c231394bdb13d2cd83217e116 Mon Sep 17 00:00:00 2001 From: AnalitikGnid Date: Thu, 29 May 2025 12:24:17 +0200 Subject: [PATCH 03/25] MyGen added, NtupleGen updated --- shipgen/CMakeLists.txt | 3 +- shipgen/LinkDef.h | 1 + shipgen/MyGenerator.cxx | 80 +++++++++++++++++++++++++++++++++++++++ shipgen/MyGenerator.h | 80 +++++++++++++++++++++++++++++++++++++++ shipgen/NtupleGenerator.h | 4 +- 5 files changed, 165 insertions(+), 3 deletions(-) create mode 100644 shipgen/MyGenerator.cxx create mode 100644 shipgen/MyGenerator.h diff --git a/shipgen/CMakeLists.txt b/shipgen/CMakeLists.txt index 63a6a7e6d5..86d607cffc 100644 --- a/shipgen/CMakeLists.txt +++ b/shipgen/CMakeLists.txt @@ -77,7 +77,8 @@ set(SRCS EvtCalcGenerator.cxx TEvtGenDecayer.cxx BeamSmearingUtils.cxx - MeanMaterialBudget.cxx) + MeanMaterialBudget.cxx + MyGenerator.cxx) set(LINKDEF LinkDef.h) set(LIBRARY_NAME ShipGen) diff --git a/shipgen/LinkDef.h b/shipgen/LinkDef.h index ed34b9f276..97cb940099 100644 --- a/shipgen/LinkDef.h +++ b/shipgen/LinkDef.h @@ -13,6 +13,7 @@ #pragma link C++ class DPPythia8Generator+; #pragma link C++ class GenieGenerator+; #pragma link C++ class NtupleGenerator+; +#pragma link C++ class MyGenerator+; #pragma link C++ class MuonBackGenerator+; #pragma link C++ class CosmicsGenerator+; #pragma link C++ class MuDISGenerator+; diff --git a/shipgen/MyGenerator.cxx b/shipgen/MyGenerator.cxx new file mode 100644 index 0000000000..628a4cb4a9 --- /dev/null +++ b/shipgen/MyGenerator.cxx @@ -0,0 +1,80 @@ +#include +#include "TROOT.h" +#include "TFile.h" +#include "FairPrimaryGenerator.h" +#include "MyGenerator.h" +#include "TDatabasePDG.h" // for TDatabasePDG +#include "TMath.h" // for Sqrt + +using std::cout; +using std::endl; +// read events from ntuples produced + +// ----- Default constructor ------------------------------------------- +MyGenerator::MyGenerator() {} +// ------------------------------------------------------------------------- +// ----- Default constructor ------------------------------------------- +Bool_t MyGenerator::Init(const char* fileName) { + return Init(fileName, 0); +} +// ----- Default constructor ------------------------------------------- +Bool_t MyGenerator::Init(const char* fileName, const int firstEvent) { + cout << "Info MyGenerator: Opening input file " << fileName << endl; + fInputFile = new TFile(fileName); + if (fInputFile->IsZombie()) { + cout << "-E MyGenerator: Error opening the Signal file" << fileName << endl; + } + fTree = (TTree *)fInputFile->Get("mytree"); + fNevents = fTree->GetEntries(); + fn = firstEvent; + fTree->SetBranchAddress("pdgcode",&id); // particle id + fTree->SetBranchAddress("vx",&vx); // position + fTree->SetBranchAddress("vy",&vy); + fTree->SetBranchAddress("vz",&vz); + fTree->SetBranchAddress("px",&px); // momentum + fTree->SetBranchAddress("py",&py); + fTree->SetBranchAddress("pz",&pz); + cout << "Init OK" << fn << " " << fNevents << endl; + return kTRUE; +} +// ------------------------------------------------------------------------- + + +// ----- Destructor ---------------------------------------------------- +MyGenerator::~MyGenerator() +{ + // cout << "destroy My" << endl; + fInputFile->Close(); + fInputFile->Delete(); + delete fInputFile; +} +// ------------------------------------------------------------------------- + +// ----- Passing the event --------------------------------------------- +Bool_t MyGenerator::ReadEvent(FairPrimaryGenerator* cpg) +{ + cout <GetEntry(fn); + fn++; + cout << "reading event "<=fNevents) { + cout << "No more input events"<GetParticle(id)->Mass(); + Double_t e = TMath::Sqrt( px[0]*px[0]+py[0]*py[0]+pz[0]*pz[0]+ mass*mass ); + cpg->AddTrack(id,px[0],py[0],pz[0],0.,0.,0.); + cout << vx[0]<< " " << vy[0]<< " " << vz[0]< +#include "TROOT.h" +#include "TFile.h" +#include "FairPrimaryGenerator.h" +#include "MyGenerator.h" +#include "TDatabasePDG.h" // for TDatabasePDG +#include "TMath.h" // for Sqrt + +using std::cout; +using std::endl; +// read events from ntuples produced + +// ----- Default constructor ------------------------------------------- +MyGenerator::MyGenerator() {} +// ------------------------------------------------------------------------- +// ----- Default constructor ------------------------------------------- +Bool_t MyGenerator::Init(const char* fileName) { + return Init(fileName, 0); +} +// ----- Default constructor ------------------------------------------- +Bool_t MyGenerator::Init(const char* fileName, const int firstEvent) { + cout << "Info MyGenerator: Opening input file " << fileName << endl; + fInputFile = new TFile(fileName); + if (fInputFile->IsZombie()) { + cout << "-E MyGenerator: Error opening the Signal file" << fileName << endl; + } + fTree = (TTree *)fInputFile->Get("mytree"); + fNevents = fTree->GetEntries(); + fn = firstEvent; + fTree->SetBranchAddress("pdgcode",&id); // particle id + fTree->SetBranchAddress("vx",&vx); // position + fTree->SetBranchAddress("vy",&vy); + fTree->SetBranchAddress("vz",&vz); + fTree->SetBranchAddress("px",&px); // momentum + fTree->SetBranchAddress("py",&py); + fTree->SetBranchAddress("pz",&pz); + cout << "Init OK" << fn << " " << fNevents << endl; + return kTRUE; +} +// ------------------------------------------------------------------------- + + +// ----- Destructor ---------------------------------------------------- +MyGenerator::~MyGenerator() +{ + // cout << "destroy My" << endl; + fInputFile->Close(); + fInputFile->Delete(); + delete fInputFile; +} +// ------------------------------------------------------------------------- + +// ----- Passing the event --------------------------------------------- +Bool_t MyGenerator::ReadEvent(FairPrimaryGenerator* cpg) +{ + cout <GetEntry(fn); + fn++; + cout << "reading event "<=fNevents) { + cout << "No more input events"<GetParticle(id)->Mass(); + Double_t e = TMath::Sqrt( px[0]*px[0]+py[0]*py[0]+pz[0]*pz[0]+ mass*mass ); + cpg->AddTrack(id,px[0],py[0],pz[0],0.,0.,0.); + cout << vx[0]<< " " << vy[0]<< " " << vz[0]< Date: Thu, 29 May 2025 12:37:44 +0200 Subject: [PATCH 04/25] Geometry tuning via csv added --- sstDecouplingTools/inputfile.root | 3 + sstDecouplingTools/produce_ttree.C | 83 ++++++++++++++++++++++++++ sstDecouplingTools/sst.csv | 8 +++ sstDecouplingTools/update_from_file.py | 27 +++++++++ 4 files changed, 121 insertions(+) create mode 100644 sstDecouplingTools/inputfile.root create mode 100644 sstDecouplingTools/produce_ttree.C create mode 100644 sstDecouplingTools/sst.csv create mode 100644 sstDecouplingTools/update_from_file.py diff --git a/sstDecouplingTools/inputfile.root b/sstDecouplingTools/inputfile.root new file mode 100644 index 0000000000..0ec1cadbf1 --- /dev/null +++ b/sstDecouplingTools/inputfile.root @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:0d0968b48fc3c6363c92dcca2bb5ec427fd5930d5927a25ef6b4e1a46e7b5eab +size 8256 diff --git a/sstDecouplingTools/produce_ttree.C b/sstDecouplingTools/produce_ttree.C new file mode 100644 index 0000000000..a52975488d --- /dev/null +++ b/sstDecouplingTools/produce_ttree.C @@ -0,0 +1,83 @@ +#include +#include +#include +#include +#include +#include "FairRootManager.h" +#include "ShipMCTrack.h" +#include "TLorentzVector.h" + +// Include class definition for strawtubesPoint if available +#include "strawtubesPoint.h" +#include "prestrawdetectorPoint.h" + +int produce_ttree() { + const char* input_file_name = "ship.conical.Pythia8-TGeant4.root"; + const char* output_file_name = "inputfile.root"; + const char* tree_name = "cbmsim"; + const char* subtree_name = "prestrawdetectorPoint"; + + // Open input file + TFile* input_file = TFile::Open(input_file_name, "READ"); + if (!input_file || input_file->IsZombie()) { + std::cerr << "Cannot open input file: " << input_file_name << std::endl; + return 1; + } + + // Get main tree + TTree* tree = dynamic_cast(input_file->Get(tree_name)); + if (!tree) { + std::cerr << "Cannot find tree named " << tree_name << " in input file." << std::endl; + return 1; + } + + // Pointer to TClonesArray of strawtubesPoint objects + TClonesArray* points = nullptr; + tree->SetBranchAddress(subtree_name, &points); + + // Output variables + double px, py, pz, vx, vy, vz; + int pdgcode; + + // Create output file and new tree + TFile* output_file = new TFile(output_file_name, "RECREATE"); + TTree* new_tree = new TTree("mytree", "Tree with PrestrawdetectorPoint"); + + new_tree->Branch("px", &px, "px/D"); + new_tree->Branch("py", &py, "py/D"); + new_tree->Branch("pz", &pz, "pz/D"); + new_tree->Branch("vx", &vx, "vx/D"); + new_tree->Branch("vy", &vy, "vy/D"); + new_tree->Branch("vz", &vz, "vz/D"); + new_tree->Branch("pdgcode", &pdgcode, "pdgcode/I"); + + Long64_t nEntries = tree->GetEntries(); + std::cout << "Total entries: " << nEntries << std::endl; + + for (Long64_t i = 0; i < nEntries; ++i) { + tree->GetEntry(i); + int nPoints = points->GetEntriesFast(); + + for (int j = 0; j < nPoints; ++j) { + auto* p = dynamic_cast(points->At(j)); + if (!p) continue; + + vx = p->GetX(); + vy = p->GetY(); + vz = p->GetZ(); + px = p->GetPx(); + py = p->GetPy(); + pz = p->GetPz(); + pdgcode = p->PdgCode(); // or p->GetPdgCode() depending on definition + + new_tree->Fill(); + } + } + new_tree->Print(); + new_tree->Write(); + output_file->Close(); + input_file->Close(); + + std::cout << "TTree '" << subtree_name << "' successfully copied to " << output_file_name << std::endl; + return 0; +} diff --git a/sstDecouplingTools/sst.csv b/sstDecouplingTools/sst.csv new file mode 100644 index 0000000000..f233af4dad --- /dev/null +++ b/sstDecouplingTools/sst.csv @@ -0,0 +1,8 @@ +key,value +options.theSeed,0 +options.nEvents,10 +ship_geo.TrackStation1.z,2598.0 +ship_geo.TrackStation2.z,2698.0 +ship_geo.TrackStation3.z,3438.0 +ship_geo.TrackStation4.z,3538.0 +ship_geo.strawtubes.ViewAngle,4.57 diff --git a/sstDecouplingTools/update_from_file.py b/sstDecouplingTools/update_from_file.py new file mode 100644 index 0000000000..149a2524ee --- /dev/null +++ b/sstDecouplingTools/update_from_file.py @@ -0,0 +1,27 @@ +import pandas as pd +import numpy as np + +def update_from_file(ship_geo, options, file): + + data = pd.read_csv(file) + + keys = data['key'].tolist() + values = data['value'].tolist() + + for key, value in zip(keys, values): + parts = key.split(".") + + target = locals()[parts[0]] + for part in parts[1:-1]: + target = getattr(target, part) + + if hasattr(target, parts[-1]): + if parts[0] == 'options': + value = int(value) + setattr(target, parts[-1], value) + print(f'{key} was set to {value}') + + # update dependent values + ship_geo.z = ship_geo.TrackStation2.z + 0.5 * (ship_geo.TrackStation3.z - ship_geo.TrackStation2.z) + ship_geo.Bfield.z = ship_geo.z + From a3be5d1fe537fe2465b43c66e9deb7b59a594571 Mon Sep 17 00:00:00 2001 From: AnalitikGnid Date: Thu, 29 May 2025 12:42:16 +0200 Subject: [PATCH 05/25] Geometry tuning via csv added --- {sstDecouplingTools => macro}/update_from_file.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename {sstDecouplingTools => macro}/update_from_file.py (100%) diff --git a/sstDecouplingTools/update_from_file.py b/macro/update_from_file.py similarity index 100% rename from sstDecouplingTools/update_from_file.py rename to macro/update_from_file.py From 097e89c94e1b295a80d72d8fd89d88c1ead045e9 Mon Sep 17 00:00:00 2001 From: AnalitikGnid Date: Thu, 29 May 2025 13:35:25 +0200 Subject: [PATCH 06/25] MyGen fixed --- shipgen/MyGenerator.h | 112 ++++++++++++++---------------------------- 1 file changed, 36 insertions(+), 76 deletions(-) diff --git a/shipgen/MyGenerator.h b/shipgen/MyGenerator.h index 628a4cb4a9..d85df70688 100644 --- a/shipgen/MyGenerator.h +++ b/shipgen/MyGenerator.h @@ -1,80 +1,40 @@ -#include -#include "TROOT.h" -#include "TFile.h" -#include "FairPrimaryGenerator.h" -#include "MyGenerator.h" -#include "TDatabasePDG.h" // for TDatabasePDG -#include "TMath.h" // for Sqrt - -using std::cout; -using std::endl; -// read events from ntuples produced - -// ----- Default constructor ------------------------------------------- -MyGenerator::MyGenerator() {} -// ------------------------------------------------------------------------- -// ----- Default constructor ------------------------------------------- -Bool_t MyGenerator::Init(const char* fileName) { - return Init(fileName, 0); -} -// ----- Default constructor ------------------------------------------- -Bool_t MyGenerator::Init(const char* fileName, const int firstEvent) { - cout << "Info MyGenerator: Opening input file " << fileName << endl; - fInputFile = new TFile(fileName); - if (fInputFile->IsZombie()) { - cout << "-E MyGenerator: Error opening the Signal file" << fileName << endl; - } - fTree = (TTree *)fInputFile->Get("mytree"); - fNevents = fTree->GetEntries(); - fn = firstEvent; - fTree->SetBranchAddress("pdgcode",&id); // particle id - fTree->SetBranchAddress("vx",&vx); // position - fTree->SetBranchAddress("vy",&vy); - fTree->SetBranchAddress("vz",&vz); - fTree->SetBranchAddress("px",&px); // momentum - fTree->SetBranchAddress("py",&py); - fTree->SetBranchAddress("pz",&pz); - cout << "Init OK" << fn << " " << fNevents << endl; - return kTRUE; -} -// ------------------------------------------------------------------------- +#ifndef MYGENERATOR_H +#define MYGENERATOR_H 1 +#include "TROOT.h" +#include "FairGenerator.h" +#include "TTree.h" // for TTree +#include "FairLogger.h" // for FairLogger, MESSAGE_ORIGIN -// ----- Destructor ---------------------------------------------------- -MyGenerator::~MyGenerator() -{ - // cout << "destroy My" << endl; - fInputFile->Close(); - fInputFile->Delete(); - delete fInputFile; -} -// ------------------------------------------------------------------------- - -// ----- Passing the event --------------------------------------------- -Bool_t MyGenerator::ReadEvent(FairPrimaryGenerator* cpg) -{ - cout <GetEntry(fn); - fn++; - cout << "reading event "<=fNevents) { - cout << "No more input events"<GetParticle(id)->Mass(); - Double_t e = TMath::Sqrt( px[0]*px[0]+py[0]*py[0]+pz[0]*pz[0]+ mass*mass ); - cpg->AddTrack(id,px[0],py[0],pz[0],0.,0.,0.); - cout << vx[0]<< " " << vy[0]<< " " << vz[0]< Date: Thu, 29 May 2025 17:36:46 +0200 Subject: [PATCH 07/25] NtupleGen fix --- shipgen/NtupleGenerator.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/shipgen/NtupleGenerator.h b/shipgen/NtupleGenerator.h index 89f0031780..de0765809e 100644 --- a/shipgen/NtupleGenerator.h +++ b/shipgen/NtupleGenerator.h @@ -29,9 +29,9 @@ class NtupleGenerator : public FairGenerator private: protected: - Float_t id,Nmeas,volid[500],procid[500],parentid; + Int_t id,Nmeas,volid[500],procid[500],parentid; Float_t Ezero,tof; - Float_t w; + Double_t w; Float_t px[500], py[500], pz[500],vx[500], vy[500], vz[500]; TFile* fInputFile; TTree* fTree; From f4847b77ec95bb9fcb2a794e3ef8c82531fd3fc5 Mon Sep 17 00:00:00 2001 From: AnalitikGnid Date: Fri, 6 Jun 2025 11:15:50 +0200 Subject: [PATCH 08/25] decoupling implementation, generator fix --- prestrawdetector/prestrawdetector.cxx | 3 +++ python/shipDet_conf.py | 4 +++- shipgen/MyGenerator.cxx | 10 +++++----- shipgen/MyGenerator.h | 2 +- sstDecouplingTools/sst.csv | 1 + 5 files changed, 13 insertions(+), 7 deletions(-) diff --git a/prestrawdetector/prestrawdetector.cxx b/prestrawdetector/prestrawdetector.cxx index 1e03832405..2a567bef75 100644 --- a/prestrawdetector/prestrawdetector.cxx +++ b/prestrawdetector/prestrawdetector.cxx @@ -183,6 +183,9 @@ void prestrawdetector::ConstructGeometry() TGeoVolume *top = gGeoManager->GetTopVolume(); InitMedium("air"); TGeoMedium *air = gGeoManager->GetMedium("air"); + InitMedium("vacuum"); + TGeoMedium *vacuum = gGeoManager->GetMedium("vacuum"); + InitMedium("ShipSens"); TGeoMedium *Se = gGeoManager->GetMedium("ShipSens"); InitMedium("aluminium"); diff --git a/python/shipDet_conf.py b/python/shipDet_conf.py index 160c95d59a..8738e84f59 100644 --- a/python/shipDet_conf.py +++ b/python/shipDet_conf.py @@ -257,7 +257,9 @@ def configure_strawtubes(yaml_file, ship_geo): ship_geo.strawtubesDigi.v_drift, ship_geo.strawtubesDigi.sigma_spatial, ) - + Prestrawdetector = ROOT.prestrawdetector('Prestrawdetector', True) + Prestrawdetector.SetZ(ship_geo.psd) + detectorList.append(Prestrawdetector) detectorList.append(strawtubes) diff --git a/shipgen/MyGenerator.cxx b/shipgen/MyGenerator.cxx index 628a4cb4a9..ba80f94e95 100644 --- a/shipgen/MyGenerator.cxx +++ b/shipgen/MyGenerator.cxx @@ -54,21 +54,21 @@ MyGenerator::~MyGenerator() Bool_t MyGenerator::ReadEvent(FairPrimaryGenerator* cpg) { cout <GetEntry(fn); fn++; cout << "reading event "<=fNevents) { cout << "No more input events"<GetParticle(id)->Mass(); - Double_t e = TMath::Sqrt( px[0]*px[0]+py[0]*py[0]+pz[0]*pz[0]+ mass*mass ); - cpg->AddTrack(id,px[0],py[0],pz[0],0.,0.,0.); - cout << vx[0]<< " " << vy[0]<< " " << vz[0]<AddTrack(id, px, py, pz, vx, vy, vz); + cout << vx<< " " << vy<< " " << vz< Date: Fri, 6 Jun 2025 12:25:09 +0200 Subject: [PATCH 09/25] README added --- sstDecouplingTools/README.md | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) create mode 100644 sstDecouplingTools/README.md diff --git a/sstDecouplingTools/README.md b/sstDecouplingTools/README.md new file mode 100644 index 0000000000..4d1363d978 --- /dev/null +++ b/sstDecouplingTools/README.md @@ -0,0 +1,27 @@ +### Geometry tuning + +To set geometrical parameters from user csv file run 'python run\_simScript.py --useCSV' + +### + +To run simulation with events from file 'python run\_simScript.py --MyGen -f _filename_'. +'_filename_' can be produced from 'ship.conical.Pythia8-TGeant4.root' with sstDecouplingTools/produce\_ttree.C + +### SST decoupling + +To run simulation on the decoupled SST (with magnet and prestrawdetector) 'python run\_simScript.py --decouple' + + +So, the pipeline is the following: + +1. Generate events using Pythia and store them with prestrawdetector +'python $FAIRSHIP/macro/run\_simScript.py --useCSV' + +2. Make input file to MyGenerator +'root $FAIRSHIP/sstDecouplingTools/produce\_ttree.C' + +3. Run simulation on the decoupled SST using stored events +'python $FAIRSHIP/macro/run\_simScript.py --useCSV --decouple --MyGen -f ./inputfile.root' + + + From 1680e39543898665b5fe4ac5fa49d41876ca2dd1 Mon Sep 17 00:00:00 2001 From: AnalitikGnid Date: Fri, 6 Jun 2025 12:40:38 +0200 Subject: [PATCH 10/25] README updated --- sstDecouplingTools/README.md | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/sstDecouplingTools/README.md b/sstDecouplingTools/README.md index 4d1363d978..8ae5ffbec9 100644 --- a/sstDecouplingTools/README.md +++ b/sstDecouplingTools/README.md @@ -1,27 +1,27 @@ ### Geometry tuning -To set geometrical parameters from user csv file run 'python run\_simScript.py --useCSV' +To set geometrical parameters from user csv file run `python run\_simScript.py --useCSV` ### -To run simulation with events from file 'python run\_simScript.py --MyGen -f _filename_'. -'_filename_' can be produced from 'ship.conical.Pythia8-TGeant4.root' with sstDecouplingTools/produce\_ttree.C +To run simulation with events from file `python run\_simScript.py --MyGen -f _filename_`. +`_filename_` can be produced from `ship.conical.Pythia8-TGeant4.root` with sstDecouplingTools/produce\_ttree.C ### SST decoupling -To run simulation on the decoupled SST (with magnet and prestrawdetector) 'python run\_simScript.py --decouple' +To run simulation on the decoupled SST (with magnet and prestrawdetector) `python run\_simScript.py --decouple` So, the pipeline is the following: 1. Generate events using Pythia and store them with prestrawdetector -'python $FAIRSHIP/macro/run\_simScript.py --useCSV' +`python $FAIRSHIP/macro/run\_simScript.py --useCSV` 2. Make input file to MyGenerator -'root $FAIRSHIP/sstDecouplingTools/produce\_ttree.C' +`root $FAIRSHIP/sstDecouplingTools/produce\_ttree.C` 3. Run simulation on the decoupled SST using stored events -'python $FAIRSHIP/macro/run\_simScript.py --useCSV --decouple --MyGen -f ./inputfile.root' +`python $FAIRSHIP/macro/run\_simScript.py --useCSV --decouple --MyGen -f ./inputfile.root` From 2b00878ee106a3c46df60a9fadbe4bc312ccca12 Mon Sep 17 00:00:00 2001 From: AnalitikGnid Date: Fri, 6 Jun 2025 12:42:40 +0200 Subject: [PATCH 11/25] README updated --- sstDecouplingTools/README.md | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/sstDecouplingTools/README.md b/sstDecouplingTools/README.md index 8ae5ffbec9..ba24cff626 100644 --- a/sstDecouplingTools/README.md +++ b/sstDecouplingTools/README.md @@ -1,27 +1,27 @@ ### Geometry tuning -To set geometrical parameters from user csv file run `python run\_simScript.py --useCSV` +To set geometrical parameters from user csv file run `python run_simScript.py --useCSV` ### -To run simulation with events from file `python run\_simScript.py --MyGen -f _filename_`. -`_filename_` can be produced from `ship.conical.Pythia8-TGeant4.root` with sstDecouplingTools/produce\_ttree.C +To run simulation with events from file `python run_simScript.py --MyGen -f _filename_`. +`_filename_` can be produced from `ship.conical.Pythia8-TGeant4.root` with sstDecouplingTools/produce_ttree.C ### SST decoupling -To run simulation on the decoupled SST (with magnet and prestrawdetector) `python run\_simScript.py --decouple` +To run simulation on the decoupled SST (also with magnet and prestrawdetector) `python run_simScript.py --decouple` So, the pipeline is the following: 1. Generate events using Pythia and store them with prestrawdetector -`python $FAIRSHIP/macro/run\_simScript.py --useCSV` +`python $FAIRSHIP/macro/run_simScript.py --useCSV` 2. Make input file to MyGenerator -`root $FAIRSHIP/sstDecouplingTools/produce\_ttree.C` +`root $FAIRSHIP/sstDecouplingTools/produce_ttree.C` 3. Run simulation on the decoupled SST using stored events -`python $FAIRSHIP/macro/run\_simScript.py --useCSV --decouple --MyGen -f ./inputfile.root` +`python $FAIRSHIP/macro/run_simScript.py --useCSV --decouple --MyGen -f ./inputfile.root` From 97cd7e0284ee0521d519a8269bc663e78459b449 Mon Sep 17 00:00:00 2001 From: AnalitikGnid Date: Fri, 6 Jun 2025 12:46:28 +0200 Subject: [PATCH 12/25] README updated --- sstDecouplingTools/README.md | 20 +++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/sstDecouplingTools/README.md b/sstDecouplingTools/README.md index ba24cff626..7ff09d9e62 100644 --- a/sstDecouplingTools/README.md +++ b/sstDecouplingTools/README.md @@ -1,27 +1,33 @@ ### Geometry tuning -To set geometrical parameters from user csv file run `python run_simScript.py --useCSV` +To set geometrical parameters from user csv file run the macro with `--useCSV` ### -To run simulation with events from file `python run_simScript.py --MyGen -f _filename_`. -`_filename_` can be produced from `ship.conical.Pythia8-TGeant4.root` with sstDecouplingTools/produce_ttree.C +To run simulation with events from file `--MyGen -f _filename_`. +`_filename_` can be produced from `ship.conical.Pythia8-TGeant4.root` with `sstDecouplingTools/produce_ttree.C` ### SST decoupling -To run simulation on the decoupled SST (also with magnet and prestrawdetector) `python run_simScript.py --decouple` +To run simulation on the decoupled SST (also with magnet and prestrawdetector) `--decouple` So, the pipeline is the following: 1. Generate events using Pythia and store them with prestrawdetector -`python $FAIRSHIP/macro/run_simScript.py --useCSV` +``` +python $FAIRSHIP/macro/run_simScript.py --useCSV +``` 2. Make input file to MyGenerator -`root $FAIRSHIP/sstDecouplingTools/produce_ttree.C` +``` +root $FAIRSHIP/sstDecouplingTools/produce_ttree.C +``` 3. Run simulation on the decoupled SST using stored events -`python $FAIRSHIP/macro/run_simScript.py --useCSV --decouple --MyGen -f ./inputfile.root` +``` +python $FAIRSHIP/macro/run_simScript.py --useCSV --decouple --MyGen -f ./inputfile.root +``` From f43d98da95107770f1fc2819857a7b6b88408667 Mon Sep 17 00:00:00 2001 From: AnalitikGnid Date: Fri, 6 Jun 2025 12:49:17 +0200 Subject: [PATCH 13/25] README updated --- sstDecouplingTools/README.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sstDecouplingTools/README.md b/sstDecouplingTools/README.md index 7ff09d9e62..8a77d0f2a4 100644 --- a/sstDecouplingTools/README.md +++ b/sstDecouplingTools/README.md @@ -4,8 +4,8 @@ To set geometrical parameters from user csv file run the macro with `--useCSV` ### -To run simulation with events from file `--MyGen -f _filename_`. -`_filename_` can be produced from `ship.conical.Pythia8-TGeant4.root` with `sstDecouplingTools/produce_ttree.C` +To run simulation with events from file `--MyGen -f filename`. +`filename` can be produced from `ship.conical.Pythia8-TGeant4.root` with `sstDecouplingTools/produce_ttree.C` ### SST decoupling From 459c6d4d2d8031a7cc4f674c01d6ede830e334b2 Mon Sep 17 00:00:00 2001 From: AnalitikGnid Date: Tue, 5 Aug 2025 11:39:12 +0200 Subject: [PATCH 14/25] csv -> json --- macro/update_from_file.py | 37 ++++++++++++++++--------------------- sstDecouplingTools/sst.csv | 2 +- sstDecouplingTools/sst.json | 10 ++++++++++ 3 files changed, 27 insertions(+), 22 deletions(-) create mode 100644 sstDecouplingTools/sst.json diff --git a/macro/update_from_file.py b/macro/update_from_file.py index 149a2524ee..de40289cb2 100644 --- a/macro/update_from_file.py +++ b/macro/update_from_file.py @@ -1,27 +1,22 @@ -import pandas as pd -import numpy as np +import json def update_from_file(ship_geo, options, file): - - data = pd.read_csv(file) + with open(file, 'r') as f: + data = json.load(f) - keys = data['key'].tolist() - values = data['value'].tolist() + if 'positions' in data: + stations = data['positions'] + ship_geo.TrackStation1.z=stations[0] + ship_geo.TrackStation2.z=stations[1] + ship_geo.TrackStation3.z=stations[2] + ship_geo.TrackStation4.z=stations[3] + print(f'Station Z-positions {stations}') + if 'angles' in data: + ship_geo.strawtubes.ViewAngle = data['angles'] + print(f"ship_geo.VA was set to {data['angles']}") + # if 'magnetic_strength' in data: + # ship_geo.magnetic_strength = data['magnetic_strength'] + # print(f'ship_geo.magnetic_strength was set to {data["magnetic_strength"]}') - for key, value in zip(keys, values): - parts = key.split(".") - - target = locals()[parts[0]] - for part in parts[1:-1]: - target = getattr(target, part) - - if hasattr(target, parts[-1]): - if parts[0] == 'options': - value = int(value) - setattr(target, parts[-1], value) - print(f'{key} was set to {value}') - - # update dependent values ship_geo.z = ship_geo.TrackStation2.z + 0.5 * (ship_geo.TrackStation3.z - ship_geo.TrackStation2.z) ship_geo.Bfield.z = ship_geo.z - diff --git a/sstDecouplingTools/sst.csv b/sstDecouplingTools/sst.csv index 6ef6203e88..2926f3f0d5 100644 --- a/sstDecouplingTools/sst.csv +++ b/sstDecouplingTools/sst.csv @@ -1,5 +1,5 @@ key,value -options.theSeed,0 +options.theSeed,1234 options.nEvents,10 ship_geo.TrackStation1.z,2598.0 ship_geo.TrackStation2.z,2698.0 diff --git a/sstDecouplingTools/sst.json b/sstDecouplingTools/sst.json new file mode 100644 index 0000000000..32494fa242 --- /dev/null +++ b/sstDecouplingTools/sst.json @@ -0,0 +1,10 @@ +{ + "positions": [ + 2598.0, + 2698.0, + 3498.0, + 3538.0 + ], + "angles": 4.57, + "magnetic_strength": 0.45 +} From 434b995ad9ac70fcb276c97f9a084584ca8827fb Mon Sep 17 00:00:00 2001 From: AnalitikGnid Date: Tue, 5 Aug 2025 11:41:04 +0200 Subject: [PATCH 15/25] HNL production info --- shipgen/HNLPythia8Generator.cxx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/shipgen/HNLPythia8Generator.cxx b/shipgen/HNLPythia8Generator.cxx index 602ecebe27..034dd39b7a 100644 --- a/shipgen/HNLPythia8Generator.cxx +++ b/shipgen/HNLPythia8Generator.cxx @@ -171,6 +171,8 @@ Bool_t HNLPythia8Generator::ReadEvent(FairPrimaryGenerator* cpg) xS = xp + lam * px; yS = yp + lam * py; zS = zp + lam * pz; + std::cout << "HNL generated at z = " << zp << " mm, decays at z = " << zS << " mm" << std::endl; + Double_t gam = e/TMath::Sqrt(e*e-p*p); Double_t beta = p/e; tS = tp + LS / beta; // units ? [mm/c] + [mm/beta] (beta is dimensionless speed, and c=1 here) From 662b4a66547b0098134ac2d37252fc1b799ca6c8 Mon Sep 17 00:00:00 2001 From: AnalitikGnid Date: Thu, 18 Sep 2025 10:06:16 +0200 Subject: [PATCH 16/25] technical commit --- prestrawdetector/prestrawdetector.cxx | 29 +++++++++++++++++++-------- 1 file changed, 21 insertions(+), 8 deletions(-) diff --git a/prestrawdetector/prestrawdetector.cxx b/prestrawdetector/prestrawdetector.cxx index 2a567bef75..882ceb1d1b 100644 --- a/prestrawdetector/prestrawdetector.cxx +++ b/prestrawdetector/prestrawdetector.cxx @@ -13,6 +13,8 @@ #include "FairRuntimeDb.h" #include "ShipDetectorList.h" #include "ShipStack.h" +#include "FairRunSim.h" +#include "FairMCEventHeader.h" #include "TClonesArray.h" #include "TVirtualMC.h" @@ -91,8 +93,14 @@ Int_t prestrawdetector::InitMedium(const char* name) Bool_t prestrawdetector::ProcessHits(FairVolume* vol) { - - std::cout << "In PH\n"; + float fEventT0 = 0.; + if (auto* run = FairRunSim::Instance()) { + if (auto* hdr = run->GetMCEventHeader()) { + fEventT0 = hdr->GetT(); + std::cout << "\nto = " << fEventT0<< "\n"; + } + } + /** This method is called from the MC stepping */ //Set parameters at entrance of volume. Reset ELoss. if ( gMC->IsTrackEntering() ) { @@ -106,12 +114,18 @@ Bool_t prestrawdetector::ProcessHits(FairVolume* vol) fELoss += gMC->Edep(); //std::cout << gMC->IsTrackExiting()<IsTrackStop()<IsTrackDisappeared(); // Create strawtubesPoint at exit of active volume - if ( gMC->IsTrackEntering() ) { + if (gMC->IsTrackExiting() || + gMC->IsTrackStop() || + gMC->IsTrackDisappeared() ) { //if (fELoss == 0. ) { return kFALSE; } + //if (!gMC->IsTrackPrimary()) return kTRUE; + //if (gMC->GetStack()->GetCurrentTrackNumber() != 0) return kTRUE; std::cout << "Entering track was found\n"; TParticle* p=gMC->GetStack()->GetCurrentTrack(); Int_t pdgCode = p->GetPdgCode(); fTrackID = gMC->GetStack()->GetCurrentTrackNumber(); + std::cout << "track = " << fTrackID << "\n"; + std::cout << "pdg = " << pdgCode << "\n"; Int_t det_uniqueId; gMC->CurrentVolID(det_uniqueId); /*if (fVolumeID == det_uniqueId) { @@ -128,6 +142,8 @@ Bool_t prestrawdetector::ProcessHits(FairVolume* vol) AddHit(fTrackID, det_uniqueId, TVector3(xmean, ymean, zmean), TVector3(fMom.Px(), fMom.Py(), fMom.Pz()), fTime, deltaTrackLength, fELoss,pdgCode); + std::cout << "Stored \n"; + std::cout << "\n fTime = " << fTime << "\n"; // Increment number of strawtubes det points in TParticle /*ShipStack* stack = (ShipStack*) gMC->GetStack(); stack->AddPoint(kStraw);*/ @@ -202,8 +218,8 @@ void prestrawdetector::ConstructGeometry() gGeoManager->SetTopVisible(); Double_t eps=0.0001; - TGeoBBox* myBox = new TGeoBBox("myBox", 300, 300, 1); - TGeoVolume* myVol = new TGeoVolume("myVol", myBox, air); + TGeoBBox* myBox = new TGeoBBox("myBox", 300, 300, 1e-6); + TGeoVolume* myVol = new TGeoVolume("myVol", myBox, vacuum); myVol->SetLineColor(kYellow); top->AddNode(myVol, 0, new TGeoTranslation(0, 0, fPSDz)); @@ -217,6 +233,3 @@ void prestrawdetector::EndOfEvent() { fprestrawdetectorPointCollection->Clear(); } - - - From 0731465a428e4801d69305de2d86bf523f8bc715 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 18 Sep 2025 09:09:38 +0000 Subject: [PATCH 17/25] style: pre-commit fixes --- macro/update_from_file.py | 2 +- prestrawdetector/prestrawdetector.cxx | 259 +++++++++++---------- prestrawdetector/prestrawdetector.h | 68 +++--- prestrawdetector/prestrawdetectorPoint.cxx | 40 ++-- prestrawdetector/prestrawdetectorPoint.h | 20 +- shipgen/MyGenerator.cxx | 2 +- sstDecouplingTools/README.md | 5 +- sstDecouplingTools/produce_ttree.C | 21 +- sstDecouplingTools/sst.json | 6 +- 9 files changed, 215 insertions(+), 208 deletions(-) diff --git a/macro/update_from_file.py b/macro/update_from_file.py index de40289cb2..28bad486fd 100644 --- a/macro/update_from_file.py +++ b/macro/update_from_file.py @@ -1,7 +1,7 @@ import json def update_from_file(ship_geo, options, file): - with open(file, 'r') as f: + with open(file) as f: data = json.load(f) if 'positions' in data: diff --git a/prestrawdetector/prestrawdetector.cxx b/prestrawdetector/prestrawdetector.cxx index 882ceb1d1b..76f2a6e796 100644 --- a/prestrawdetector/prestrawdetector.cxx +++ b/prestrawdetector/prestrawdetector.cxx @@ -1,33 +1,31 @@ #include "prestrawdetector.h" -#include "prestrawdetectorPoint.h" -#include "FairVolume.h" -#include "FairGeoVolume.h" -#include "FairGeoNode.h" -#include "FairRootManager.h" -#include "FairGeoLoader.h" +#include "FairGeoBuilder.h" #include "FairGeoInterface.h" +#include "FairGeoLoader.h" #include "FairGeoMedia.h" -#include "FairGeoBuilder.h" +#include "FairGeoNode.h" +#include "FairGeoVolume.h" +#include "FairMCEventHeader.h" +#include "FairRootManager.h" #include "FairRun.h" +#include "FairRunSim.h" #include "FairRuntimeDb.h" +#include "FairVolume.h" #include "ShipDetectorList.h" #include "ShipStack.h" -#include "FairRunSim.h" -#include "FairMCEventHeader.h" - #include "TClonesArray.h" -#include "TVirtualMC.h" -#include "TGeoManager.h" #include "TGeoBBox.h" #include "TGeoCompositeShape.h" -#include "TGeoTube.h" +#include "TGeoManager.h" #include "TGeoMaterial.h" #include "TGeoMedium.h" +#include "TGeoTube.h" #include "TMath.h" #include "TParticle.h" #include "TVector3.h" - +#include "TVirtualMC.h" +#include "prestrawdetectorPoint.h" prestrawdetector::prestrawdetector() : FairDetector("prestrawdetector", kTRUE) @@ -55,175 +53,182 @@ prestrawdetector::prestrawdetector(const char* name, Bool_t active) , fprestrawdetectorPointCollection(new TClonesArray("prestrawdetectorPoint")) {} -prestrawdetector::~prestrawdetector() { +prestrawdetector::~prestrawdetector() +{ if (fprestrawdetectorPointCollection) { - fprestrawdetectorPointCollection->Delete(); - delete fprestrawdetectorPointCollection; - } + fprestrawdetectorPointCollection->Delete(); + delete fprestrawdetectorPointCollection; + } } void prestrawdetector::Initialize() { - FairDetector::Initialize(); -// FairRuntimeDb* rtdb= FairRun::Instance()->GetRuntimeDb(); -// vetoGeoPar* par=(vetoGeoPar*)(rtdb->getContainer("vetoGeoPar")); + FairDetector::Initialize(); + // FairRuntimeDb* rtdb= FairRun::Instance()->GetRuntimeDb(); + // vetoGeoPar* par=(vetoGeoPar*)(rtdb->getContainer("vetoGeoPar")); } // ----- Private method InitMedium Int_t prestrawdetector::InitMedium(const char* name) { - static FairGeoLoader *geoLoad=FairGeoLoader::Instance(); - static FairGeoInterface *geoFace=geoLoad->getGeoInterface(); - static FairGeoMedia *media=geoFace->getMedia(); - static FairGeoBuilder *geoBuild=geoLoad->getGeoBuilder(); - - FairGeoMedium *ShipMedium=media->getMedium(name); - - if (!ShipMedium) - { - Fatal("InitMedium","Material %s not defined in media file.", name); - return -1111; - } - TGeoMedium* medium=gGeoManager->GetMedium(name); - if (medium!=NULL) - return ShipMedium->getMediumIndex(); - - return geoBuild->createMedium(ShipMedium); + static FairGeoLoader* geoLoad = FairGeoLoader::Instance(); + static FairGeoInterface* geoFace = geoLoad->getGeoInterface(); + static FairGeoMedia* media = geoFace->getMedia(); + static FairGeoBuilder* geoBuild = geoLoad->getGeoBuilder(); + + FairGeoMedium* ShipMedium = media->getMedium(name); + + if (!ShipMedium) { + Fatal("InitMedium", "Material %s not defined in media file.", name); + return -1111; + } + TGeoMedium* medium = gGeoManager->GetMedium(name); + if (medium != NULL) + return ShipMedium->getMediumIndex(); + + return geoBuild->createMedium(ShipMedium); } Bool_t prestrawdetector::ProcessHits(FairVolume* vol) { - float fEventT0 = 0.; - if (auto* run = FairRunSim::Instance()) { - if (auto* hdr = run->GetMCEventHeader()) { - fEventT0 = hdr->GetT(); - std::cout << "\nto = " << fEventT0<< "\n"; - } - } - + float fEventT0 = 0.; + if (auto* run = FairRunSim::Instance()) { + if (auto* hdr = run->GetMCEventHeader()) { + fEventT0 = hdr->GetT(); + std::cout << "\nto = " << fEventT0 << "\n"; + } + } + /** This method is called from the MC stepping */ - //Set parameters at entrance of volume. Reset ELoss. - if ( gMC->IsTrackEntering() ) { - fELoss = 0.; - fTime = gMC->TrackTime() * 1.0e09; - fLength = gMC->TrackLength(); - gMC->TrackPosition(fPos); - gMC->TrackMomentum(fMom); - } - // Sum energy loss for all steps in the active volume - fELoss += gMC->Edep(); - //std::cout << gMC->IsTrackExiting()<IsTrackStop()<IsTrackDisappeared(); - // Create strawtubesPoint at exit of active volume - if (gMC->IsTrackExiting() || - gMC->IsTrackStop() || - gMC->IsTrackDisappeared() ) { - //if (fELoss == 0. ) { return kFALSE; } - //if (!gMC->IsTrackPrimary()) return kTRUE; - //if (gMC->GetStack()->GetCurrentTrackNumber() != 0) return kTRUE; - std::cout << "Entering track was found\n"; - TParticle* p=gMC->GetStack()->GetCurrentTrack(); - Int_t pdgCode = p->GetPdgCode(); - fTrackID = gMC->GetStack()->GetCurrentTrackNumber(); - std::cout << "track = " << fTrackID << "\n"; - std::cout << "pdg = " << pdgCode << "\n"; - Int_t det_uniqueId; - gMC->CurrentVolID(det_uniqueId); - /*if (fVolumeID == det_uniqueId) { - //std::cout << pdgCode<< " same volume again ? "<< straw_uniqueId << " exit:" << gMC->IsTrackExiting() << " stop:" << gMC->IsTrackStop() << " disappeared:" << gMC->IsTrackDisappeared()<< std::endl; - return kTRUE; }*/ - fVolumeID = det_uniqueId; - // # d = |pq . u x v|/|u x v| - TLorentzVector Pos; - gMC->TrackPosition(Pos); - Double_t xmean = (fPos.X()+Pos.X())/2. ; - Double_t ymean = (fPos.Y()+Pos.Y())/2. ; - Double_t zmean = (fPos.Z()+Pos.Z())/2. ; - Double_t deltaTrackLength = gMC->TrackLength() - fLength; - AddHit(fTrackID, det_uniqueId, TVector3(xmean, ymean, zmean), - TVector3(fMom.Px(), fMom.Py(), fMom.Pz()), fTime, deltaTrackLength, - fELoss,pdgCode); - std::cout << "Stored \n"; - std::cout << "\n fTime = " << fTime << "\n"; - // Increment number of strawtubes det points in TParticle - /*ShipStack* stack = (ShipStack*) gMC->GetStack(); - stack->AddPoint(kStraw);*/ - } + // Set parameters at entrance of volume. Reset ELoss. + if (gMC->IsTrackEntering()) { + fELoss = 0.; + fTime = gMC->TrackTime() * 1.0e09; + fLength = gMC->TrackLength(); + gMC->TrackPosition(fPos); + gMC->TrackMomentum(fMom); + } + // Sum energy loss for all steps in the active volume + fELoss += gMC->Edep(); + // std::cout << gMC->IsTrackExiting()<IsTrackStop()<IsTrackDisappeared(); + // Create strawtubesPoint at exit of active volume + if (gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()) { + // if (fELoss == 0. ) { return kFALSE; } + // if (!gMC->IsTrackPrimary()) return kTRUE; + // if (gMC->GetStack()->GetCurrentTrackNumber() != 0) return kTRUE; + std::cout << "Entering track was found\n"; + TParticle* p = gMC->GetStack()->GetCurrentTrack(); + Int_t pdgCode = p->GetPdgCode(); + fTrackID = gMC->GetStack()->GetCurrentTrackNumber(); + std::cout << "track = " << fTrackID << "\n"; + std::cout << "pdg = " << pdgCode << "\n"; + Int_t det_uniqueId; + gMC->CurrentVolID(det_uniqueId); + /*if (fVolumeID == det_uniqueId) { + //std::cout << pdgCode<< " same volume again ? "<< straw_uniqueId << " exit:" << gMC->IsTrackExiting() << " + stop:" << gMC->IsTrackStop() << " disappeared:" << gMC->IsTrackDisappeared()<< std::endl; return kTRUE; }*/ + fVolumeID = det_uniqueId; + // # d = |pq . u x v|/|u x v| + TLorentzVector Pos; + gMC->TrackPosition(Pos); + Double_t xmean = (fPos.X() + Pos.X()) / 2.; + Double_t ymean = (fPos.Y() + Pos.Y()) / 2.; + Double_t zmean = (fPos.Z() + Pos.Z()) / 2.; + Double_t deltaTrackLength = gMC->TrackLength() - fLength; + AddHit(fTrackID, + det_uniqueId, + TVector3(xmean, ymean, zmean), + TVector3(fMom.Px(), fMom.Py(), fMom.Pz()), + fTime, + deltaTrackLength, + fELoss, + pdgCode); + std::cout << "Stored \n"; + std::cout << "\n fTime = " << fTime << "\n"; + // Increment number of strawtubes det points in TParticle + /*ShipStack* stack = (ShipStack*) gMC->GetStack(); + stack->AddPoint(kStraw);*/ + } return kTRUE; } -prestrawdetectorPoint* prestrawdetector::AddHit(Int_t trackID, Int_t detID, - TVector3 pos, TVector3 mom, - Double_t time, Double_t length, - Double_t eLoss, Int_t pdgCode) +prestrawdetectorPoint* prestrawdetector::AddHit(Int_t trackID, + Int_t detID, + TVector3 pos, + TVector3 mom, + Double_t time, + Double_t length, + Double_t eLoss, + Int_t pdgCode) { - TClonesArray& clref = *fprestrawdetectorPointCollection; - Int_t size = clref.GetEntriesFast(); - std::cout << "adding hit detid " << detID <Register("prestrawdetectorPoint", "prestrawdetector", - fprestrawdetectorPointCollection, kTRUE); - std::cout << "reg2"; - + /** This will create a branch in the output tree called + strawtubesPoint, setting the last parameter to kFALSE means: + this collection will not be written to the file, it will exist + only during the simulation. + */ + std::cout << "reg1"; + FairRootManager::Instance()->Register( + "prestrawdetectorPoint", "prestrawdetector", fprestrawdetectorPointCollection, kTRUE); + std::cout << "reg2"; } TClonesArray* prestrawdetector::GetCollection(Int_t iColl) const { - if (iColl == 0) { return fprestrawdetectorPointCollection; } - else { return NULL; } + if (iColl == 0) { + return fprestrawdetectorPointCollection; + } else { + return NULL; + } } void prestrawdetector::Reset() { - fprestrawdetectorPointCollection->Clear(); + fprestrawdetectorPointCollection->Clear(); } void prestrawdetector::SetZ(Double_t z0) { - fPSDz = z0; //! z-position of veto station + fPSDz = z0; //! z-position of veto station } void prestrawdetector::ConstructGeometry() { - TGeoVolume *top = gGeoManager->GetTopVolume(); + TGeoVolume* top = gGeoManager->GetTopVolume(); InitMedium("air"); - TGeoMedium *air = gGeoManager->GetMedium("air"); + TGeoMedium* air = gGeoManager->GetMedium("air"); InitMedium("vacuum"); - TGeoMedium *vacuum = gGeoManager->GetMedium("vacuum"); + TGeoMedium* vacuum = gGeoManager->GetMedium("vacuum"); InitMedium("ShipSens"); - TGeoMedium *Se = gGeoManager->GetMedium("ShipSens"); + TGeoMedium* Se = gGeoManager->GetMedium("ShipSens"); InitMedium("aluminium"); - TGeoMedium *Al = gGeoManager->GetMedium("aluminium"); + TGeoMedium* Al = gGeoManager->GetMedium("aluminium"); InitMedium("mylar"); - TGeoMedium *mylar = gGeoManager->GetMedium("mylar"); + TGeoMedium* mylar = gGeoManager->GetMedium("mylar"); InitMedium("STTmix9010_2bar"); - TGeoMedium *sttmix9010_2bar = gGeoManager->GetMedium("STTmix9010_2bar"); + TGeoMedium* sttmix9010_2bar = gGeoManager->GetMedium("STTmix9010_2bar"); InitMedium("tungsten"); - TGeoMedium *tungsten = gGeoManager->GetMedium("tungsten"); - + TGeoMedium* tungsten = gGeoManager->GetMedium("tungsten"); gGeoManager->SetVisLevel(4); gGeoManager->SetTopVisible(); - Double_t eps=0.0001; + Double_t eps = 0.0001; TGeoBBox* myBox = new TGeoBBox("myBox", 300, 300, 1e-6); TGeoVolume* myVol = new TGeoVolume("myVol", myBox, vacuum); myVol->SetLineColor(kYellow); top->AddNode(myVol, 0, new TGeoTranslation(0, 0, fPSDz)); - AddSensitiveVolume(myVol); + AddSensitiveVolume(myVol); // Position it somewhere // top->AddNode(myVol, 1, new TGeoTranslation(0, 0, 500)); @@ -231,5 +236,5 @@ void prestrawdetector::ConstructGeometry() void prestrawdetector::EndOfEvent() { - fprestrawdetectorPointCollection->Clear(); + fprestrawdetectorPointCollection->Clear(); } diff --git a/prestrawdetector/prestrawdetector.h b/prestrawdetector/prestrawdetector.h index 9df4a3e6d8..4df7f9aed8 100644 --- a/prestrawdetector/prestrawdetector.h +++ b/prestrawdetector/prestrawdetector.h @@ -3,9 +3,9 @@ #include "FairDetector.h" #include "FairVolume.h" -#include "TVector3.h" #include "TClonesArray.h" #include "TLorentzVector.h" +#include "TVector3.h" class prestrawdetectorPoint; class FairVolume; @@ -19,57 +19,59 @@ class prestrawdetector : public FairDetector virtual ~prestrawdetector(); /** Initialization of the detector is done here */ - virtual void Initialize(); //done + virtual void Initialize(); // done + + virtual Bool_t ProcessHits(FairVolume* v = 0); // done - virtual Bool_t ProcessHits( FairVolume* v=0);//done + prestrawdetectorPoint* AddHit(Int_t trackID, + Int_t detID, + TVector3 pos, + TVector3 mom, + Double_t time, + Double_t length, + Double_t eLoss, + Int_t pdgCode); - prestrawdetectorPoint* AddHit(Int_t trackID, Int_t detID, - TVector3 pos, TVector3 mom, - Double_t time, Double_t length, - Double_t eLoss, Int_t pdgCode); - /** Registers the produced collections in FAIRRootManager. */ - virtual void Register(); + virtual void Register(); /** Gets the produced collections */ - virtual TClonesArray* GetCollection(Int_t iColl) const ; //done + virtual TClonesArray* GetCollection(Int_t iColl) const; // done /** has to be called after each event to reset the containers */ - virtual void Reset(); + virtual void Reset(); void SetZ(Double_t z0); void ConstructGeometry(); - - virtual void CopyClones( TClonesArray* cl1, TClonesArray* cl2 , - Int_t offset) {;} - virtual void SetSpecialPhysicsCuts() {;} - virtual void EndOfEvent(); - virtual void FinishPrimary() {;} - virtual void FinishRun() {;} - virtual void BeginPrimary() {;} - virtual void PostTrack() {;} - virtual void PreTrack() {;} - virtual void BeginEvent() {;} + virtual void CopyClones(TClonesArray* cl1, TClonesArray* cl2, Int_t offset) { ; } + virtual void SetSpecialPhysicsCuts() { ; } + virtual void EndOfEvent(); + virtual void FinishPrimary() { ; } + virtual void FinishRun() { ; } + virtual void BeginPrimary() { ; } + virtual void PostTrack() { ; } + virtual void PreTrack() { ; } + virtual void BeginEvent() { ; } private: - Double_t fPSDz; - std::string fMedium; - Int_t fTrackID; //! track index - Int_t fVolumeID; //! volume id - TLorentzVector fPos; //! position at entrance - TLorentzVector fMom; //! momentum at entrance - Double_t fTime; //! time - Double_t fLength; //! length - Double_t fELoss; - TClonesArray* fprestrawdetectorPointCollection; // collection of hit points + Double_t fPSDz; + std::string fMedium; + Int_t fTrackID; //! track index + Int_t fVolumeID; //! volume id + TLorentzVector fPos; //! position at entrance + TLorentzVector fMom; //! momentum at entrance + Double_t fTime; //! time + Double_t fLength; //! length + Double_t fELoss; + TClonesArray* fprestrawdetectorPointCollection; // collection of hit points prestrawdetector(const prestrawdetector&); prestrawdetector& operator=(const prestrawdetector&); Int_t InitMedium(const char* name); - ClassDef(prestrawdetector,2) + ClassDef(prestrawdetector, 2) }; #endif diff --git a/prestrawdetector/prestrawdetectorPoint.cxx b/prestrawdetector/prestrawdetectorPoint.cxx index db609a81e5..383452c409 100644 --- a/prestrawdetector/prestrawdetectorPoint.cxx +++ b/prestrawdetector/prestrawdetectorPoint.cxx @@ -5,38 +5,38 @@ using std::cout; using std::endl; - // ----- Default constructor ------------------------------------------- prestrawdetectorPoint::prestrawdetectorPoint() - : FairMCPoint() -{ -} + : FairMCPoint() +{} // ------------------------------------------------------------------------- // ----- Standard constructor ------------------------------------------ -prestrawdetectorPoint::prestrawdetectorPoint(Int_t trackID, Int_t detID, - TVector3 pos, TVector3 mom, - Double_t tof, Double_t length, - Double_t eLoss,Int_t pdgcode) - : FairMCPoint(trackID, detID, pos, mom, tof, length, eLoss), fPdgCode(pdgcode) -{ -} +prestrawdetectorPoint::prestrawdetectorPoint(Int_t trackID, + Int_t detID, + TVector3 pos, + TVector3 mom, + Double_t tof, + Double_t length, + Double_t eLoss, + Int_t pdgcode) + : FairMCPoint(trackID, detID, pos, mom, tof, length, eLoss) + , fPdgCode(pdgcode) +{} // ------------------------------------------------------------------------- // ----- Destructor ---------------------------------------------------- -prestrawdetectorPoint::~prestrawdetectorPoint() { } +prestrawdetectorPoint::~prestrawdetectorPoint() {} // ------------------------------------------------------------------------- // ----- Public method Print ------------------------------------------- void prestrawdetectorPoint::Print() const { - cout << "-I- prestrawdetectorPoint: prestrawdetector point for track " << fTrackID - << " in detector " << fDetectorID << endl; - cout << " Position (" << fX << ", " << fY << ", " << fZ - << ") cm" << endl; - cout << " Momentum (" << fPx << ", " << fPy << ", " << fPz - << ") GeV" << endl; - cout << " Time " << fTime << " ns, Length " << fLength << "pdgid" << fPdgCode - << " cm, Energy loss " << fELoss*1.0e06 << " keV" << endl; + cout << "-I- prestrawdetectorPoint: prestrawdetector point for track " << fTrackID << " in detector " << fDetectorID + << endl; + cout << " Position (" << fX << ", " << fY << ", " << fZ << ") cm" << endl; + cout << " Momentum (" << fPx << ", " << fPy << ", " << fPz << ") GeV" << endl; + cout << " Time " << fTime << " ns, Length " << fLength << "pdgid" << fPdgCode << " cm, Energy loss " + << fELoss * 1.0e06 << " keV" << endl; } // ------------------------------------------------------------------------- diff --git a/prestrawdetector/prestrawdetectorPoint.h b/prestrawdetector/prestrawdetectorPoint.h index 1bb1d988a1..88ba5ff6b7 100644 --- a/prestrawdetector/prestrawdetectorPoint.h +++ b/prestrawdetector/prestrawdetectorPoint.h @@ -1,16 +1,13 @@ #ifndef PRESTRAWDETECTORPOINT_H #define PRESTRAWDETECTORPOINT_H 1 - #include "FairMCPoint.h" - #include "TObject.h" #include "TVector3.h" class prestrawdetectorPoint : public FairMCPoint { public: - /** Default constructor **/ prestrawdetectorPoint(); @@ -23,16 +20,21 @@ class prestrawdetectorPoint : public FairMCPoint *@param length Track length since creation [cm] *@param eLoss Energy deposit [GeV] **/ - prestrawdetectorPoint(Int_t trackID, Int_t detID, TVector3 pos, TVector3 mom, - Double_t tof, Double_t length, Double_t eLoss, Int_t pdgcode); + prestrawdetectorPoint(Int_t trackID, + Int_t detID, + TVector3 pos, + TVector3 mom, + Double_t tof, + Double_t length, + Double_t eLoss, + Int_t pdgcode); /** Destructor **/ virtual ~prestrawdetectorPoint(); /** Output to screen **/ virtual void Print() const; - Int_t PdgCode() const {return fPdgCode;} - + Int_t PdgCode() const { return fPdgCode; } private: /** Copy constructor **/ @@ -40,9 +42,7 @@ class prestrawdetectorPoint : public FairMCPoint prestrawdetectorPoint operator=(const prestrawdetectorPoint& point); Int_t fPdgCode; - ClassDef(prestrawdetectorPoint,2); - - + ClassDef(prestrawdetectorPoint, 2); }; #endif diff --git a/shipgen/MyGenerator.cxx b/shipgen/MyGenerator.cxx index ba80f94e95..12332d75a3 100644 --- a/shipgen/MyGenerator.cxx +++ b/shipgen/MyGenerator.cxx @@ -61,7 +61,7 @@ Bool_t MyGenerator::ReadEvent(FairPrimaryGenerator* cpg) } if (fn>=fNevents) { cout << "No more input events"< -#include -#include -#include -#include #include "FairRootManager.h" #include "ShipMCTrack.h" #include "TLorentzVector.h" +#include +#include +#include +#include +#include + // Include class definition for strawtubesPoint if available -#include "strawtubesPoint.h" #include "prestrawdetectorPoint.h" +#include "strawtubesPoint.h" -int produce_ttree() { +int produce_ttree() +{ const char* input_file_name = "ship.conical.Pythia8-TGeant4.root"; const char* output_file_name = "inputfile.root"; const char* tree_name = "cbmsim"; @@ -60,7 +62,8 @@ int produce_ttree() { for (int j = 0; j < nPoints; ++j) { auto* p = dynamic_cast(points->At(j)); - if (!p) continue; + if (!p) + continue; vx = p->GetX(); vy = p->GetY(); @@ -68,7 +71,7 @@ int produce_ttree() { px = p->GetPx(); py = p->GetPy(); pz = p->GetPz(); - pdgcode = p->PdgCode(); // or p->GetPdgCode() depending on definition + pdgcode = p->PdgCode(); // or p->GetPdgCode() depending on definition new_tree->Fill(); } diff --git a/sstDecouplingTools/sst.json b/sstDecouplingTools/sst.json index 32494fa242..844372d585 100644 --- a/sstDecouplingTools/sst.json +++ b/sstDecouplingTools/sst.json @@ -1,8 +1,8 @@ { "positions": [ - 2598.0, - 2698.0, - 3498.0, + 2598.0, + 2698.0, + 3498.0, 3538.0 ], "angles": 4.57, From 0cdd5dd8a8cf153a033ab39b3b3d9942f8488a03 Mon Sep 17 00:00:00 2001 From: AnalitikGnid Date: Tue, 23 Sep 2025 12:14:21 +0200 Subject: [PATCH 18/25] minor fixes and edits in prestrawdetector --- prestrawdetector/prestrawdetector.cxx | 276 +++++++++++++------------- prestrawdetector/prestrawdetector.h | 75 +++---- 2 files changed, 183 insertions(+), 168 deletions(-) diff --git a/prestrawdetector/prestrawdetector.cxx b/prestrawdetector/prestrawdetector.cxx index 76f2a6e796..5036837ddb 100644 --- a/prestrawdetector/prestrawdetector.cxx +++ b/prestrawdetector/prestrawdetector.cxx @@ -1,31 +1,33 @@ #include "prestrawdetector.h" +#include "prestrawdetectorPoint.h" -#include "FairGeoBuilder.h" -#include "FairGeoInterface.h" -#include "FairGeoLoader.h" -#include "FairGeoMedia.h" -#include "FairGeoNode.h" +#include "FairVolume.h" #include "FairGeoVolume.h" -#include "FairMCEventHeader.h" +#include "FairGeoNode.h" #include "FairRootManager.h" +#include "FairGeoLoader.h" +#include "FairGeoInterface.h" +#include "FairGeoMedia.h" +#include "FairGeoBuilder.h" #include "FairRun.h" -#include "FairRunSim.h" #include "FairRuntimeDb.h" -#include "FairVolume.h" #include "ShipDetectorList.h" #include "ShipStack.h" +#include "FairRunSim.h" +#include "FairMCEventHeader.h" + #include "TClonesArray.h" +#include "TVirtualMC.h" +#include "TGeoManager.h" #include "TGeoBBox.h" #include "TGeoCompositeShape.h" -#include "TGeoManager.h" +#include "TGeoTube.h" #include "TGeoMaterial.h" #include "TGeoMedium.h" -#include "TGeoTube.h" #include "TMath.h" #include "TParticle.h" #include "TVector3.h" -#include "TVirtualMC.h" -#include "prestrawdetectorPoint.h" + prestrawdetector::prestrawdetector() : FairDetector("prestrawdetector", kTRUE) @@ -53,182 +55,190 @@ prestrawdetector::prestrawdetector(const char* name, Bool_t active) , fprestrawdetectorPointCollection(new TClonesArray("prestrawdetectorPoint")) {} -prestrawdetector::~prestrawdetector() -{ +prestrawdetector::prestrawdetector(std::string medium) + : FairDetector("prestrawdetector", kTRUE, kStraw) + , fTrackID(-1) + , fVolumeID(-1) + , fPos() + , fMom() + , fTime(-1.) + , fLength(-1.) + , fELoss(-1) + , fMedium(medium) + , fprestrawdetectorPointCollection(new TClonesArray("prestrawdetectorPoint")) +{} + +prestrawdetector::~prestrawdetector() { if (fprestrawdetectorPointCollection) { - fprestrawdetectorPointCollection->Delete(); - delete fprestrawdetectorPointCollection; - } + fprestrawdetectorPointCollection->Delete(); + delete fprestrawdetectorPointCollection; + } } void prestrawdetector::Initialize() { - FairDetector::Initialize(); - // FairRuntimeDb* rtdb= FairRun::Instance()->GetRuntimeDb(); - // vetoGeoPar* par=(vetoGeoPar*)(rtdb->getContainer("vetoGeoPar")); + FairDetector::Initialize(); } // ----- Private method InitMedium Int_t prestrawdetector::InitMedium(const char* name) { - static FairGeoLoader* geoLoad = FairGeoLoader::Instance(); - static FairGeoInterface* geoFace = geoLoad->getGeoInterface(); - static FairGeoMedia* media = geoFace->getMedia(); - static FairGeoBuilder* geoBuild = geoLoad->getGeoBuilder(); - - FairGeoMedium* ShipMedium = media->getMedium(name); - - if (!ShipMedium) { - Fatal("InitMedium", "Material %s not defined in media file.", name); - return -1111; - } - TGeoMedium* medium = gGeoManager->GetMedium(name); - if (medium != NULL) - return ShipMedium->getMediumIndex(); - - return geoBuild->createMedium(ShipMedium); + static FairGeoLoader *geoLoad=FairGeoLoader::Instance(); + static FairGeoInterface *geoFace=geoLoad->getGeoInterface(); + static FairGeoMedia *media=geoFace->getMedia(); + static FairGeoBuilder *geoBuild=geoLoad->getGeoBuilder(); + + FairGeoMedium *ShipMedium=media->getMedium(name); + + if (!ShipMedium) + { + Fatal("InitMedium","Material %s not defined in media file.", name); + return -1111; + } + TGeoMedium* medium=gGeoManager->GetMedium(name); + if (medium!=NULL) + return ShipMedium->getMediumIndex(); + + return geoBuild->createMedium(ShipMedium); } Bool_t prestrawdetector::ProcessHits(FairVolume* vol) { - float fEventT0 = 0.; - if (auto* run = FairRunSim::Instance()) { - if (auto* hdr = run->GetMCEventHeader()) { - fEventT0 = hdr->GetT(); - std::cout << "\nto = " << fEventT0 << "\n"; - } - } - + float fEventT0 = 0.; + if (auto* run = FairRunSim::Instance()) { + if (auto* hdr = run->GetMCEventHeader()) { + fEventT0 = hdr->GetT(); + } + } + /** This method is called from the MC stepping */ - // Set parameters at entrance of volume. Reset ELoss. - if (gMC->IsTrackEntering()) { - fELoss = 0.; - fTime = gMC->TrackTime() * 1.0e09; - fLength = gMC->TrackLength(); - gMC->TrackPosition(fPos); - gMC->TrackMomentum(fMom); - } - // Sum energy loss for all steps in the active volume - fELoss += gMC->Edep(); - // std::cout << gMC->IsTrackExiting()<IsTrackStop()<IsTrackDisappeared(); - // Create strawtubesPoint at exit of active volume - if (gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()) { - // if (fELoss == 0. ) { return kFALSE; } - // if (!gMC->IsTrackPrimary()) return kTRUE; - // if (gMC->GetStack()->GetCurrentTrackNumber() != 0) return kTRUE; - std::cout << "Entering track was found\n"; - TParticle* p = gMC->GetStack()->GetCurrentTrack(); - Int_t pdgCode = p->GetPdgCode(); - fTrackID = gMC->GetStack()->GetCurrentTrackNumber(); - std::cout << "track = " << fTrackID << "\n"; - std::cout << "pdg = " << pdgCode << "\n"; - Int_t det_uniqueId; - gMC->CurrentVolID(det_uniqueId); - /*if (fVolumeID == det_uniqueId) { - //std::cout << pdgCode<< " same volume again ? "<< straw_uniqueId << " exit:" << gMC->IsTrackExiting() << " - stop:" << gMC->IsTrackStop() << " disappeared:" << gMC->IsTrackDisappeared()<< std::endl; return kTRUE; }*/ - fVolumeID = det_uniqueId; - // # d = |pq . u x v|/|u x v| - TLorentzVector Pos; - gMC->TrackPosition(Pos); - Double_t xmean = (fPos.X() + Pos.X()) / 2.; - Double_t ymean = (fPos.Y() + Pos.Y()) / 2.; - Double_t zmean = (fPos.Z() + Pos.Z()) / 2.; - Double_t deltaTrackLength = gMC->TrackLength() - fLength; - AddHit(fTrackID, - det_uniqueId, - TVector3(xmean, ymean, zmean), - TVector3(fMom.Px(), fMom.Py(), fMom.Pz()), - fTime, - deltaTrackLength, - fELoss, - pdgCode); - std::cout << "Stored \n"; - std::cout << "\n fTime = " << fTime << "\n"; - // Increment number of strawtubes det points in TParticle - /*ShipStack* stack = (ShipStack*) gMC->GetStack(); - stack->AddPoint(kStraw);*/ - } + //Set parameters at entrance of volume. Reset ELoss. + if ( gMC->IsTrackEntering() ) { + fELoss = 0.; + fTime = gMC->TrackTime() * 1.0e09; + fLength = gMC->TrackLength(); + gMC->TrackPosition(fPos); + gMC->TrackMomentum(fMom); + } + // Sum energy loss for all steps in the active volume + fELoss += gMC->Edep(); + //std::cout << gMC->IsTrackExiting()<IsTrackStop()<IsTrackDisappeared(); + // Create strawtubesPoint at exit of active volume + if (gMC->IsTrackEntering() ) { + //if (fELoss == 0. ) { return kFALSE; } + //if (!gMC->IsTrackPrimary()) return kTRUE; + //if (gMC->GetStack()->GetCurrentTrackNumber() != 0) return kTRUE; + TParticle* p=gMC->GetStack()->GetCurrentTrack(); + Int_t pdgCode = p->GetPdgCode(); + fTrackID = gMC->GetStack()->GetCurrentTrackNumber(); + Int_t det_uniqueId; + gMC->CurrentVolID(det_uniqueId); + /*if (fVolumeID == det_uniqueId) { + //std::cout << pdgCode<< " same volume again ? "<< straw_uniqueId << " exit:" << gMC->IsTrackExiting() << " stop:" << gMC->IsTrackStop() << " disappeared:" << gMC->IsTrackDisappeared()<< std::endl; + return kTRUE; }*/ + fVolumeID = det_uniqueId; + // # d = |pq . u x v|/|u x v| + TLorentzVector Pos; + gMC->TrackPosition(Pos); + Double_t xmean = (fPos.X()+Pos.X())/2. ; + Double_t ymean = (fPos.Y()+Pos.Y())/2. ; + Double_t zmean = (fPos.Z()+Pos.Z())/2. ; + Double_t deltaTrackLength = gMC->TrackLength() - fLength; + AddHit(fTrackID, det_uniqueId, TVector3(xmean, ymean, zmean), + TVector3(fMom.Px(), fMom.Py(), fMom.Pz()), fTime, deltaTrackLength, + fELoss,pdgCode); + + // Increment number of strawtubes det points in TParticle + /*ShipStack* stack = (ShipStack*) gMC->GetStack(); + stack->AddPoint(kStraw);*/ + } return kTRUE; } -prestrawdetectorPoint* prestrawdetector::AddHit(Int_t trackID, - Int_t detID, - TVector3 pos, - TVector3 mom, - Double_t time, - Double_t length, - Double_t eLoss, - Int_t pdgCode) +prestrawdetectorPoint* prestrawdetector::AddHit(Int_t trackID, Int_t detID, + TVector3 pos, TVector3 mom, + Double_t time, Double_t length, + Double_t eLoss, Int_t pdgCode) { - TClonesArray& clref = *fprestrawdetectorPointCollection; - Int_t size = clref.GetEntriesFast(); - std::cout << "adding hit detid " << detID << std::endl; - return new (clref[size]) prestrawdetectorPoint(trackID, detID, pos, mom, time, length, eLoss, pdgCode); + TClonesArray& clref = *fprestrawdetectorPointCollection; + Int_t size = clref.GetEntriesFast(); + return new(clref[size]) prestrawdetectorPoint(trackID, detID, pos, mom, + time, length, eLoss, pdgCode); } void prestrawdetector::Register() { - /** This will create a branch in the output tree called - strawtubesPoint, setting the last parameter to kFALSE means: - this collection will not be written to the file, it will exist - only during the simulation. - */ - std::cout << "reg1"; - FairRootManager::Instance()->Register( - "prestrawdetectorPoint", "prestrawdetector", fprestrawdetectorPointCollection, kTRUE); - std::cout << "reg2"; + /** This will create a branch in the output tree called + strawtubesPoint, setting the last parameter to kFALSE means: + this collection will not be written to the file, it will exist + only during the simulation. + */ + FairRootManager::Instance()->Register("prestrawdetectorPoint", "prestrawdetector", + fprestrawdetectorPointCollection, kTRUE); } TClonesArray* prestrawdetector::GetCollection(Int_t iColl) const { - if (iColl == 0) { - return fprestrawdetectorPointCollection; - } else { - return NULL; - } + if (iColl == 0) { return fprestrawdetectorPointCollection; } + else { return NULL; } } void prestrawdetector::Reset() { - fprestrawdetectorPointCollection->Clear(); + fprestrawdetectorPointCollection->Clear(); } void prestrawdetector::SetZ(Double_t z0) { - fPSDz = z0; //! z-position of veto station + fPSDz = z0; //! z-position of veto station +} + +void prestrawdetector::SetX(Double_t x) +{ + fPSDx = x; //! x-half-size +} + +void prestrawdetector::SetY(Double_t y) +{ + fPSDy = y; //! y-half-size +} + +void prestrawdetector::SetThickness(Double_t thickness) +{ + fPSD_Thickness = thickness; //! thickness } void prestrawdetector::ConstructGeometry() { - TGeoVolume* top = gGeoManager->GetTopVolume(); + TGeoVolume *top = gGeoManager->GetTopVolume(); InitMedium("air"); - TGeoMedium* air = gGeoManager->GetMedium("air"); + TGeoMedium *air = gGeoManager->GetMedium("air"); InitMedium("vacuum"); - TGeoMedium* vacuum = gGeoManager->GetMedium("vacuum"); + TGeoMedium *vacuum = gGeoManager->GetMedium("vacuum"); InitMedium("ShipSens"); - TGeoMedium* Se = gGeoManager->GetMedium("ShipSens"); + TGeoMedium *Se = gGeoManager->GetMedium("ShipSens"); InitMedium("aluminium"); - TGeoMedium* Al = gGeoManager->GetMedium("aluminium"); + TGeoMedium *Al = gGeoManager->GetMedium("aluminium"); InitMedium("mylar"); - TGeoMedium* mylar = gGeoManager->GetMedium("mylar"); + TGeoMedium *mylar = gGeoManager->GetMedium("mylar"); InitMedium("STTmix9010_2bar"); - TGeoMedium* sttmix9010_2bar = gGeoManager->GetMedium("STTmix9010_2bar"); + TGeoMedium *sttmix9010_2bar = gGeoManager->GetMedium("STTmix9010_2bar"); InitMedium("tungsten"); - TGeoMedium* tungsten = gGeoManager->GetMedium("tungsten"); + TGeoMedium *tungsten = gGeoManager->GetMedium("tungsten"); + gGeoManager->SetVisLevel(4); gGeoManager->SetTopVisible(); - Double_t eps = 0.0001; + Double_t eps=0.0001; - TGeoBBox* myBox = new TGeoBBox("myBox", 300, 300, 1e-6); + TGeoBBox* myBox = new TGeoBBox("myBox", fPSDx, fPSDy, fPSD_Thickness); TGeoVolume* myVol = new TGeoVolume("myVol", myBox, vacuum); myVol->SetLineColor(kYellow); top->AddNode(myVol, 0, new TGeoTranslation(0, 0, fPSDz)); - AddSensitiveVolume(myVol); + AddSensitiveVolume(myVol); // Position it somewhere // top->AddNode(myVol, 1, new TGeoTranslation(0, 0, 500)); @@ -236,5 +246,5 @@ void prestrawdetector::ConstructGeometry() void prestrawdetector::EndOfEvent() { - fprestrawdetectorPointCollection->Clear(); + fprestrawdetectorPointCollection->Clear(); } diff --git a/prestrawdetector/prestrawdetector.h b/prestrawdetector/prestrawdetector.h index 4df7f9aed8..718af354df 100644 --- a/prestrawdetector/prestrawdetector.h +++ b/prestrawdetector/prestrawdetector.h @@ -3,9 +3,9 @@ #include "FairDetector.h" #include "FairVolume.h" +#include "TVector3.h" #include "TClonesArray.h" #include "TLorentzVector.h" -#include "TVector3.h" class prestrawdetectorPoint; class FairVolume; @@ -16,62 +16,67 @@ class prestrawdetector : public FairDetector public: prestrawdetector(const char* name, Bool_t active); prestrawdetector(); + prestrawdetector(std::string medium); virtual ~prestrawdetector(); /** Initialization of the detector is done here */ - virtual void Initialize(); // done - - virtual Bool_t ProcessHits(FairVolume* v = 0); // done + virtual void Initialize(); //done - prestrawdetectorPoint* AddHit(Int_t trackID, - Int_t detID, - TVector3 pos, - TVector3 mom, - Double_t time, - Double_t length, - Double_t eLoss, - Int_t pdgCode); + virtual Bool_t ProcessHits( FairVolume* v=0);//done + prestrawdetectorPoint* AddHit(Int_t trackID, Int_t detID, + TVector3 pos, TVector3 mom, + Double_t time, Double_t length, + Double_t eLoss, Int_t pdgCode); + /** Registers the produced collections in FAIRRootManager. */ - virtual void Register(); + virtual void Register(); /** Gets the produced collections */ - virtual TClonesArray* GetCollection(Int_t iColl) const; // done + virtual TClonesArray* GetCollection(Int_t iColl) const ; //done /** has to be called after each event to reset the containers */ - virtual void Reset(); + virtual void Reset(); void SetZ(Double_t z0); + void SetX(Double_t x); + void SetY(Double_t y); + void SetThickness(Double_t thickness); void ConstructGeometry(); + + virtual void CopyClones( TClonesArray* cl1, TClonesArray* cl2 , + Int_t offset) {;} + virtual void SetSpecialPhysicsCuts() {;} + virtual void EndOfEvent(); + virtual void FinishPrimary() {;} + virtual void FinishRun() {;} + virtual void BeginPrimary() {;} + virtual void PostTrack() {;} + virtual void PreTrack() {;} + virtual void BeginEvent() {;} - virtual void CopyClones(TClonesArray* cl1, TClonesArray* cl2, Int_t offset) { ; } - virtual void SetSpecialPhysicsCuts() { ; } - virtual void EndOfEvent(); - virtual void FinishPrimary() { ; } - virtual void FinishRun() { ; } - virtual void BeginPrimary() { ; } - virtual void PostTrack() { ; } - virtual void PreTrack() { ; } - virtual void BeginEvent() { ; } private: - Double_t fPSDz; - std::string fMedium; - Int_t fTrackID; //! track index - Int_t fVolumeID; //! volume id - TLorentzVector fPos; //! position at entrance - TLorentzVector fMom; //! momentum at entrance - Double_t fTime; //! time - Double_t fLength; //! length - Double_t fELoss; - TClonesArray* fprestrawdetectorPointCollection; // collection of hit points + Double_t fPSDz; + Double_t fPSDx; //! x(half of the size) in cm + Double_t fPSDy; //! y(half of the size) in cm + Double_t fPSD_Thickness; //! thickness in cm + std::string fMedium; + Int_t fTrackID; //! track index + Int_t fVolumeID; //! volume id + TLorentzVector fPos; //! position at entrance + TLorentzVector fMom; //! momentum at entrance + Double_t fTime; //! time + Double_t fLength; //! length + Double_t fELoss; + TClonesArray* fprestrawdetectorPointCollection; // collection of hit points prestrawdetector(const prestrawdetector&); prestrawdetector& operator=(const prestrawdetector&); Int_t InitMedium(const char* name); - ClassDef(prestrawdetector, 2) + ClassDef(prestrawdetector,2) }; #endif From aeeae33e121a9982755cc5ec11e79c4c7e1fd788 Mon Sep 17 00:00:00 2001 From: AnalitikGnid Date: Tue, 23 Sep 2025 12:23:45 +0200 Subject: [PATCH 19/25] prestrawdetector configuration updated --- geometry/prestrawdetector_config.yaml | 8 ++++++++ 1 file changed, 8 insertions(+) create mode 100644 geometry/prestrawdetector_config.yaml diff --git a/geometry/prestrawdetector_config.yaml b/geometry/prestrawdetector_config.yaml new file mode 100644 index 0000000000..73c7304077 --- /dev/null +++ b/geometry/prestrawdetector_config.yaml @@ -0,0 +1,8 @@ +# Geometry configuration of PSD (prestrawdetector) in FairShip + +PSD: + width: 300 # Aperture width (x) in cm (half length) + height: 300 # Aperture height (y) in cm (half length) + wall_thickness: 0.001 # Straw wall thickness in cm + z_position: 8350.0 # Z position in cm + material: "air" # material From 41eebf64341f581076ffa8ef015f2d1a8dbe0238 Mon Sep 17 00:00:00 2001 From: AnalitikGnid Date: Tue, 23 Sep 2025 12:25:48 +0200 Subject: [PATCH 20/25] minor edits --- shipgen/MyGenerator.cxx | 2 -- 1 file changed, 2 deletions(-) diff --git a/shipgen/MyGenerator.cxx b/shipgen/MyGenerator.cxx index 12332d75a3..e3a4bbc67e 100644 --- a/shipgen/MyGenerator.cxx +++ b/shipgen/MyGenerator.cxx @@ -34,7 +34,6 @@ Bool_t MyGenerator::Init(const char* fileName, const int firstEvent) { fTree->SetBranchAddress("px",&px); // momentum fTree->SetBranchAddress("py",&py); fTree->SetBranchAddress("pz",&pz); - cout << "Init OK" << fn << " " << fNevents << endl; return kTRUE; } // ------------------------------------------------------------------------- @@ -68,7 +67,6 @@ Bool_t MyGenerator::ReadEvent(FairPrimaryGenerator* cpg) Double_t mass = pdgBase->GetParticle(id)->Mass(); Double_t e = TMath::Sqrt( px*px+py*py+pz*pz+ mass*mass ); cpg->AddTrack(id, px, py, pz, vx, vy, vz); - cout << vx<< " " << vy<< " " << vz< Date: Tue, 23 Sep 2025 12:27:29 +0200 Subject: [PATCH 21/25] interface improved --- macro/run_simScript.py | 125 ++++++++++++++++++++---------------- macro/update_config.py | 25 ++++++++ macro/update_from_file.py | 22 ------- sstDecouplingTools/sst.json | 17 ++--- 4 files changed, 103 insertions(+), 86 deletions(-) create mode 100644 macro/update_config.py delete mode 100644 macro/update_from_file.py diff --git a/macro/run_simScript.py b/macro/run_simScript.py index 37ad2c276a..9fd090aefc 100755 --- a/macro/run_simScript.py +++ b/macro/run_simScript.py @@ -1,7 +1,4 @@ #!/usr/bin/env python -# SPDX-License-Identifier: LGPL-3.0-or-later -# SPDX-FileCopyrightText: Copyright CERN for the benefit of the SHiP Collaboration - import os import sys import ROOT @@ -9,12 +6,12 @@ import shipunit as u import shipRoot_conf import rootUtils as ut -import geometry_config +from ShipGeoConfig import ConfigRegistry from argparse import ArgumentParser from array import array from backports import tdirectory634 DownScaleDiMuon = False - +from update_config import update_config # Default HNL parameters theHNLMass = 1.0*u.GeV theProductionCouplings = theDecayCouplings = None @@ -39,6 +36,21 @@ inputFile = "$EOSSHIP/eos/experiment/ship/data/Charm/Cascade-parp16-MSTP82-1-MSEL4-978Bpot.root" defaultInputFile = True +globalDesigns = { + '2023' : { + 'dy' : 6., + 'caloDesign' : 3, + 'strawDesign' : 10 + }, + '2025' : { + 'dy' : 6., + 'ds' : 8, + 'caloDesign' : 2, + 'strawDesign' : 10 + }, +} +default = '2025' + parser = ArgumentParser() group = parser.add_mutually_exclusive_group() @@ -46,8 +58,8 @@ parser.add_argument("--Pythia6", dest="pythia6", help="Use Pythia6", action="store_true") parser.add_argument("--Pythia8", dest="pythia8", help="Use Pythia8", action="store_true") parser.add_argument("--EvtGenDecayer", dest="evtgen_decayer", help="Use TEvtGenDecayer for J/psi and other quarkonium decays", action="store_true") -subparsers = parser.add_subparsers(dest="command", help="Which mode to run") # === PG subcommand === +subparsers = parser.add_subparsers(dest="command", help="Which mode to run") pg_parser = subparsers.add_parser("PG", help="Use Particle Gun") pg_parser.add_argument( @@ -87,19 +99,9 @@ "--Dy", dest="Dy", type=float, help="size of the full uniform spread of PG ypos: (Vy - Dy/2, Vy + Dy/2)" ) -# === End of PG commands === -# === Genie subcommand === -genie_parser = subparsers.add_parser("Genie", help="Genie for reading and processing neutrino interactions") -genie_parser.add_argument( - "--z_start_nu", dest="z_start_nu", default=2844.2850, type=float, - help="Genie neutrino start z start position (default=2844.2850 cm)" -) -genie_parser.add_argument( - "--z_end_nu", dest="z_end_nu", default=3180.4350, type=float, - help="Genie neutrino end z position (default=3180.4350 cm)" -) -# === End of Genie subcommand === +# === Enf of PG commands === parser.add_argument("-A", help="b: signal from b, c: from c (default), bc: from Bc, or inclusive", default='c') +parser.add_argument("--Genie", dest="genie", help="Genie for reading and processing neutrino interactions", action="store_true") parser.add_argument("--NuRadio", dest="nuradio", help="misuse GenieGenerator for neutrino radiography and geometry timing test", action="store_true") parser.add_argument("--Ntuple", dest="ntuple", help="Use ntuple as input", action="store_true") parser.add_argument("--MuonBack", dest="muonback", help="Generate events from muon background file, --Cosmics=0 for cosmic generator data", action="store_true") @@ -128,10 +130,15 @@ group.add_argument("-f", dest="inputFile", help="Input file if not default file", default=False) parser.add_argument("-g", dest="geofile", help="geofile for muon shield geometry, for experts only", default=None) parser.add_argument("-o", "--output", dest="outputDir", help="Output directory", default=".") -parser.add_argument("-Y", dest="dy", help="max height of vacuum tank", default=6.0, type=float) +parser.add_argument("-Y", dest="dy", help="max height of vacuum tank", default=globalDesigns[default]['dy']) +parser.add_argument("--caloDesign", + help="0=ECAL/HCAL TP 2=splitCal 3=ECAL/ passive HCAL", + default=globalDesigns[default]['caloDesign'], + type=int, + choices=[0,2,3]) parser.add_argument("--strawDesign", help="Tracker station frame material: 4=aluminium; 10=steel (default)", - default=10, + default=globalDesigns[default]['strawDesign'], type=int, choices=[4,10]) parser.add_argument("-F", dest="deepCopy", help="default = False: copy only stable particles to stack, except for HNL events", action="store_true") @@ -165,16 +172,25 @@ parser.add_argument("--SND", dest="SND", help="Activate SND.", action='store_true') parser.add_argument( "--SND_design", - help="Choose SND design(s) among [1,2,...] or 'all' to enable all. 1: EmulsionTarget, 2: MTC + SiliconTarget", + help="Choose SND design(s) among [1,2,...] or 'all' to enable all. 1: EmulsionTarget, 2: MTC", nargs='+', default=[2], ) parser.add_argument("--noSND", dest="SND", help="Deactivate SND. NOOP, as it currently defaults to off.", action='store_false') parser.add_argument("--target-yaml", help="Path to the yaml target config file", default=os.path.expandvars("$FAIRSHIP/geometry/target_config_Jun25.yaml")) - +# sst-decoupling flags +parser.add_argument("--decouple", dest="decouple", help="Decoupling flag", required=False, default=False, action='store_true') +parser.add_argument("--useJSON", dest="json", help="Geometry tuning from json file", required=False, default=False, action='store_true') +parser.add_argument("--prestraw", dest="PSD", help="Prestrawdetector flag", required=False, default=False, action='store_true') +parser.add_argument("--ttreegen", dest="ttreegen", help="Use events from ttree as input", required=False, action="store_true") options = parser.parse_args() +if options.json: + update_config(options, os.path.expandvars('$FAIRSHIP/sstDecouplingTools/sst.json'), 'options') + + + # Handle SND_design: allow 'all' (case-insensitive) or list of ints available_snd_designs = [1, 2] # Extend this list as new designs are added if any(str(x).lower() == 'all' for x in options.SND_design): @@ -222,7 +238,7 @@ print("cannot have HNL and SUSY or DP at the same time, abort") sys.exit(2) -if (options.command == "Genie" or options.nuradio) and defaultInputFile: +if (options.genie or options.nuradio) and defaultInputFile: inputFile = "$EOSSHIP/eos/experiment/ship/data/GenieEvents/genie-nu_mu.root" if options.mudis and defaultInputFile: print('input file required if simEngine = muonDIS') @@ -246,8 +262,10 @@ ROOT.gInterpreter.ProcessLine('fair::Logger::SetConsoleSeverity("debug1");') elif options.debug == 3: ROOT.gInterpreter.ProcessLine('fair::Logger::SetConsoleSeverity("debug2");') -ship_geo = geometry_config.create_config( +ship_geo = ConfigRegistry.loadpy( + "$FAIRSHIP/geometry/geometry_config.py", Yheight=options.dy, + CaloDesign=options.caloDesign, strawDesign=options.strawDesign, muShieldGeo=options.geofile, shieldName=options.shieldName, @@ -256,15 +274,19 @@ SND_design=options.SND_design, TARGET_YAML=options.target_yaml ) +print(ship_geo) +if options.json: + update_config(ship_geo, os.path.expandvars('$FAIRSHIP/sstDecouplingTools/sst.json'), 'geometry') + +if options.PSD: + ship_geo.PSD = options.PSD +ship_geo.decouple = options.decouple # Output file name, add dy to be able to setup geometry with ambiguities. if options.command == "PG": tag = f"PG_{options.pID}-{mcEngine}" -elif options.command == "Genie": - tag = f"Genie-{mcEngine}" - simEngine = "Genie" else: - for g in ["pythia8", "evtcalc", "pythia6", "nuradio", "ntuple", "muonback", "mudis", "fixedTarget", "cosmics"]: + for g in ["pythia8", "evtcalc", "pythia6", "genie", "nuradio", "ntuple", "muonback", "mudis", "fixedTarget", "cosmics", "ttreegen"]: if getattr(options, g): simEngine = g.capitalize() break @@ -308,6 +330,7 @@ primGen = ROOT.FairPrimaryGenerator() if options.pythia8: primGen.SetTarget(ship_geo.target.z0, 0.) + # -----Pythia8-------------------------------------- if HNL or options.RPVSUSY: P8gen = ROOT.HNLPythia8Generator() @@ -361,9 +384,9 @@ P8gen.UseExternalFile(inputFile, options.firstEvent) # Use geometry constants instead of fragile TGeo navigation P8gen.SetTargetCoordinates(ship_geo.target.z0, ship_geo.target.z0 + ship_geo.target.length) -# pion on proton 500GeV -# P8gen.SetMom(500.*u.GeV) -# P8gen.SetId(-211) + # pion on proton 500GeV + # P8gen.SetMom(500.*u.GeV) + # P8gen.SetId(-211) primGen.AddGenerator(P8gen) if options.fixedTarget: HNL = False @@ -388,6 +411,15 @@ P6gen.SetTarget("gamma/mu+","n0") # default "gamma/mu-","p+" primGen.AddGenerator(P6gen) +# -----TTree generator-------------------------------------- + +if options.ttreegen: + primGen.SetTarget(0.0, 0.0) + print(f"Opening input file for ttree generator: {inputFile}") + TtreeGen = ROOT.MyGenerator() + TtreeGen.Init(inputFile, options.firstEvent) + primGen.AddGenerator(TtreeGen) + options.nEvents = TtreeGen.GetNevents() # -----EvtCalc-------------------------------------- if options.evtcalc: primGen.SetTarget(0.0, 0.0) @@ -436,13 +468,13 @@ options.nEvents = min(options.nEvents,DISgen.GetNevents()) print('Generate ',options.nEvents,' with DIS input', ' first event',options.firstEvent) # -----Neutrino Background------------------------ -if options.command == "Genie": +if options.genie: # Genie ut.checkFileExists(inputFile) primGen.SetTarget(0., 0.) # do not interfere with GenieGenerator Geniegen = ROOT.GenieGenerator() Geniegen.Init(inputFile,options.firstEvent) - Geniegen.SetPositions(ship_geo.target.z0, options.z_start_nu, options.z_end_nu) + Geniegen.SetPositions(ship_geo.target.z0, ship_geo.tauMudet.zMudetC-5*u.m, ship_geo.TrackStation2.z) primGen.AddGenerator(Geniegen) options.nEvents = min(options.nEvents,Geniegen.GetNevents()) run.SetPythiaDecayer("DecayConfigNuAge.C") @@ -587,7 +619,7 @@ if options.print_fields: geomGeant4.printVMCFields() geomGeant4.printWeightsandFields(onlyWithField = True,\ - exclude=['DecayVolume','Tr1','Tr2','Tr3','Tr4','Veto','MuonDetector','SplitCal']) + exclude=['DecayVolume','Tr1','Tr2','Tr3','Tr4','Veto','Ecal','Hcal','MuonDetector','SplitCal']) # Plot the field example #fieldMaker.plotField(1, ROOT.TVector3(-9000.0, 6000.0, 50.0), ROOT.TVector3(-300.0, 300.0, 6.0), 'Bzx.png') #fieldMaker.plotField(2, ROOT.TVector3(-9000.0, 6000.0, 50.0), ROOT.TVector3(-400.0, 400.0, 6.0), 'Bzy.png') @@ -674,11 +706,13 @@ branches.Add(ROOT.TObjString('TTPoint')) branches.Add(ROOT.TObjString('ScoringPoint')) branches.Add(ROOT.TObjString('strawtubesPoint')) + branches.Add(ROOT.TObjString('EcalPoint')) + branches.Add(ROOT.TObjString('sEcalPointLite')) + branches.Add(ROOT.TObjString('smuonPoint')) branches.Add(ROOT.TObjString('TimeDetPoint')) branches.Add(ROOT.TObjString('MCEventHeader')) branches.Add(ROOT.TObjString('UpstreamTaggerPoint')) branches.Add(ROOT.TObjString('MTCdetPoint')) - branches.Add(ROOT.TObjString('SiliconTargetPoint')) branches.Add(ROOT.TObjString('sGeoTracks')) sTree.AutoSave() @@ -721,27 +755,6 @@ os.replace(temp_filename, outFile) print("Successfully added DISCrossSection to the output file:", outFile) -if options.command == "Genie": - # breakpoint() - # Copy Genie (gst TTree) information to the output file - f_input = ROOT.TFile.Open(inputFile, "READ") - print("check") - gst = f_input.gst - - selection_string = "(Entry$ >= "+str(options.firstEvent)+")" - if (options.firstEvent + options.nEvents) < gst.GetEntries() : - selection_string += "&&(Entry$ < "+str(options.firstEvent + options.nEvents)+")" - - # Reopen output file - f_output = ROOT.TFile.Open(outFile, "UPDATE") - - # Copy only the events used in this run - gst_copy = gst.CopyTree(selection_string) - gst_copy.Write() - - f_input.Close() - f_output.Close() - # ------------------------------------------------------------------------ import checkMagFields diff --git a/macro/update_config.py b/macro/update_config.py new file mode 100644 index 0000000000..0601addab0 --- /dev/null +++ b/macro/update_config.py @@ -0,0 +1,25 @@ +import json +from ShipGeoConfig import AttrDict + +def update_config(global_config, my_config, name): + with open(my_config, 'r') as f: + data = json.load(f) + new_values = AttrDict(data[name]) + if name == 'geometry': + stations = new_values.positions + global_config.TrackStation1.z=stations[0] + global_config.TrackStation2.z=stations[1] + global_config.TrackStation3.z=stations[2] + global_config.TrackStation4.z=stations[3] + print(f'Station Z-positions {stations}') + #global_config.strawtubes.ViewAngle = new_values.angles + #print(f"ship_geo.VA was set to {data['angles']}") + # if 'magnetic_strength' in data: + # ship_geo.magnetic_strength = data['magnetic_strength'] + # print(f'ship_geo.magnetic_strength was set to {data["magnetic_strength"]}') + + global_config.z = global_config.TrackStation2.z + 0.5 * (global_config.TrackStation3.z - global_config.TrackStation2.z) + global_config.Bfield.z = global_config.z + if name == 'options': + global_config.theSeed = new_values.theSeed + global_config.nEvents = new_values.nEvents diff --git a/macro/update_from_file.py b/macro/update_from_file.py deleted file mode 100644 index 28bad486fd..0000000000 --- a/macro/update_from_file.py +++ /dev/null @@ -1,22 +0,0 @@ -import json - -def update_from_file(ship_geo, options, file): - with open(file) as f: - data = json.load(f) - - if 'positions' in data: - stations = data['positions'] - ship_geo.TrackStation1.z=stations[0] - ship_geo.TrackStation2.z=stations[1] - ship_geo.TrackStation3.z=stations[2] - ship_geo.TrackStation4.z=stations[3] - print(f'Station Z-positions {stations}') - if 'angles' in data: - ship_geo.strawtubes.ViewAngle = data['angles'] - print(f"ship_geo.VA was set to {data['angles']}") - # if 'magnetic_strength' in data: - # ship_geo.magnetic_strength = data['magnetic_strength'] - # print(f'ship_geo.magnetic_strength was set to {data["magnetic_strength"]}') - - ship_geo.z = ship_geo.TrackStation2.z + 0.5 * (ship_geo.TrackStation3.z - ship_geo.TrackStation2.z) - ship_geo.Bfield.z = ship_geo.z diff --git a/sstDecouplingTools/sst.json b/sstDecouplingTools/sst.json index 844372d585..5602172bf5 100644 --- a/sstDecouplingTools/sst.json +++ b/sstDecouplingTools/sst.json @@ -1,10 +1,11 @@ { - "positions": [ - 2598.0, - 2698.0, - 3498.0, - 3538.0 - ], - "angles": 4.57, - "magnetic_strength": 0.45 + "geometry": { + "positions": [8407.0, 8607.0, 9307.0, 9507.0], + "angles": 4.57, + "magnetic_strength": 0.45 + }, + "options": { + "theSeed": 42, + "nEvents": 11 + } } From ea7a46b4ccac6fe5c570d86feffa786facc92cd3 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 23 Sep 2025 10:21:39 +0000 Subject: [PATCH 22/25] style: pre-commit fixes --- prestrawdetector/prestrawdetector.cxx | 242 +++++++++++++------------- prestrawdetector/prestrawdetector.h | 74 ++++---- 2 files changed, 163 insertions(+), 153 deletions(-) diff --git a/prestrawdetector/prestrawdetector.cxx b/prestrawdetector/prestrawdetector.cxx index 5036837ddb..3e848ee79b 100644 --- a/prestrawdetector/prestrawdetector.cxx +++ b/prestrawdetector/prestrawdetector.cxx @@ -1,33 +1,31 @@ #include "prestrawdetector.h" -#include "prestrawdetectorPoint.h" -#include "FairVolume.h" -#include "FairGeoVolume.h" -#include "FairGeoNode.h" -#include "FairRootManager.h" -#include "FairGeoLoader.h" +#include "FairGeoBuilder.h" #include "FairGeoInterface.h" +#include "FairGeoLoader.h" #include "FairGeoMedia.h" -#include "FairGeoBuilder.h" +#include "FairGeoNode.h" +#include "FairGeoVolume.h" +#include "FairMCEventHeader.h" +#include "FairRootManager.h" #include "FairRun.h" +#include "FairRunSim.h" #include "FairRuntimeDb.h" +#include "FairVolume.h" #include "ShipDetectorList.h" #include "ShipStack.h" -#include "FairRunSim.h" -#include "FairMCEventHeader.h" - #include "TClonesArray.h" -#include "TVirtualMC.h" -#include "TGeoManager.h" #include "TGeoBBox.h" #include "TGeoCompositeShape.h" -#include "TGeoTube.h" +#include "TGeoManager.h" #include "TGeoMaterial.h" #include "TGeoMedium.h" +#include "TGeoTube.h" #include "TMath.h" #include "TParticle.h" #include "TVector3.h" - +#include "TVirtualMC.h" +#include "prestrawdetectorPoint.h" prestrawdetector::prestrawdetector() : FairDetector("prestrawdetector", kTRUE) @@ -68,177 +66,187 @@ prestrawdetector::prestrawdetector(std::string medium) , fprestrawdetectorPointCollection(new TClonesArray("prestrawdetectorPoint")) {} -prestrawdetector::~prestrawdetector() { +prestrawdetector::~prestrawdetector() +{ if (fprestrawdetectorPointCollection) { - fprestrawdetectorPointCollection->Delete(); - delete fprestrawdetectorPointCollection; - } + fprestrawdetectorPointCollection->Delete(); + delete fprestrawdetectorPointCollection; + } } void prestrawdetector::Initialize() { - FairDetector::Initialize(); + FairDetector::Initialize(); } // ----- Private method InitMedium Int_t prestrawdetector::InitMedium(const char* name) { - static FairGeoLoader *geoLoad=FairGeoLoader::Instance(); - static FairGeoInterface *geoFace=geoLoad->getGeoInterface(); - static FairGeoMedia *media=geoFace->getMedia(); - static FairGeoBuilder *geoBuild=geoLoad->getGeoBuilder(); - - FairGeoMedium *ShipMedium=media->getMedium(name); - - if (!ShipMedium) - { - Fatal("InitMedium","Material %s not defined in media file.", name); - return -1111; - } - TGeoMedium* medium=gGeoManager->GetMedium(name); - if (medium!=NULL) - return ShipMedium->getMediumIndex(); - - return geoBuild->createMedium(ShipMedium); + static FairGeoLoader* geoLoad = FairGeoLoader::Instance(); + static FairGeoInterface* geoFace = geoLoad->getGeoInterface(); + static FairGeoMedia* media = geoFace->getMedia(); + static FairGeoBuilder* geoBuild = geoLoad->getGeoBuilder(); + + FairGeoMedium* ShipMedium = media->getMedium(name); + + if (!ShipMedium) { + Fatal("InitMedium", "Material %s not defined in media file.", name); + return -1111; + } + TGeoMedium* medium = gGeoManager->GetMedium(name); + if (medium != NULL) + return ShipMedium->getMediumIndex(); + + return geoBuild->createMedium(ShipMedium); } Bool_t prestrawdetector::ProcessHits(FairVolume* vol) { - float fEventT0 = 0.; - if (auto* run = FairRunSim::Instance()) { - if (auto* hdr = run->GetMCEventHeader()) { - fEventT0 = hdr->GetT(); - } - } - + float fEventT0 = 0.; + if (auto* run = FairRunSim::Instance()) { + if (auto* hdr = run->GetMCEventHeader()) { + fEventT0 = hdr->GetT(); + } + } + /** This method is called from the MC stepping */ - //Set parameters at entrance of volume. Reset ELoss. - if ( gMC->IsTrackEntering() ) { - fELoss = 0.; - fTime = gMC->TrackTime() * 1.0e09; - fLength = gMC->TrackLength(); - gMC->TrackPosition(fPos); - gMC->TrackMomentum(fMom); - } - // Sum energy loss for all steps in the active volume - fELoss += gMC->Edep(); - //std::cout << gMC->IsTrackExiting()<IsTrackStop()<IsTrackDisappeared(); - // Create strawtubesPoint at exit of active volume - if (gMC->IsTrackEntering() ) { - //if (fELoss == 0. ) { return kFALSE; } - //if (!gMC->IsTrackPrimary()) return kTRUE; - //if (gMC->GetStack()->GetCurrentTrackNumber() != 0) return kTRUE; - TParticle* p=gMC->GetStack()->GetCurrentTrack(); - Int_t pdgCode = p->GetPdgCode(); - fTrackID = gMC->GetStack()->GetCurrentTrackNumber(); - Int_t det_uniqueId; - gMC->CurrentVolID(det_uniqueId); - /*if (fVolumeID == det_uniqueId) { - //std::cout << pdgCode<< " same volume again ? "<< straw_uniqueId << " exit:" << gMC->IsTrackExiting() << " stop:" << gMC->IsTrackStop() << " disappeared:" << gMC->IsTrackDisappeared()<< std::endl; - return kTRUE; }*/ - fVolumeID = det_uniqueId; - // # d = |pq . u x v|/|u x v| - TLorentzVector Pos; - gMC->TrackPosition(Pos); - Double_t xmean = (fPos.X()+Pos.X())/2. ; - Double_t ymean = (fPos.Y()+Pos.Y())/2. ; - Double_t zmean = (fPos.Z()+Pos.Z())/2. ; - Double_t deltaTrackLength = gMC->TrackLength() - fLength; - AddHit(fTrackID, det_uniqueId, TVector3(xmean, ymean, zmean), - TVector3(fMom.Px(), fMom.Py(), fMom.Pz()), fTime, deltaTrackLength, - fELoss,pdgCode); - - // Increment number of strawtubes det points in TParticle - /*ShipStack* stack = (ShipStack*) gMC->GetStack(); - stack->AddPoint(kStraw);*/ - } + // Set parameters at entrance of volume. Reset ELoss. + if (gMC->IsTrackEntering()) { + fELoss = 0.; + fTime = gMC->TrackTime() * 1.0e09; + fLength = gMC->TrackLength(); + gMC->TrackPosition(fPos); + gMC->TrackMomentum(fMom); + } + // Sum energy loss for all steps in the active volume + fELoss += gMC->Edep(); + // std::cout << gMC->IsTrackExiting()<IsTrackStop()<IsTrackDisappeared(); + // Create strawtubesPoint at exit of active volume + if (gMC->IsTrackEntering()) { + // if (fELoss == 0. ) { return kFALSE; } + // if (!gMC->IsTrackPrimary()) return kTRUE; + // if (gMC->GetStack()->GetCurrentTrackNumber() != 0) return kTRUE; + TParticle* p = gMC->GetStack()->GetCurrentTrack(); + Int_t pdgCode = p->GetPdgCode(); + fTrackID = gMC->GetStack()->GetCurrentTrackNumber(); + Int_t det_uniqueId; + gMC->CurrentVolID(det_uniqueId); + /*if (fVolumeID == det_uniqueId) { + //std::cout << pdgCode<< " same volume again ? "<< straw_uniqueId << " exit:" << gMC->IsTrackExiting() << " + stop:" << gMC->IsTrackStop() << " disappeared:" << gMC->IsTrackDisappeared()<< std::endl; return kTRUE; }*/ + fVolumeID = det_uniqueId; + // # d = |pq . u x v|/|u x v| + TLorentzVector Pos; + gMC->TrackPosition(Pos); + Double_t xmean = (fPos.X() + Pos.X()) / 2.; + Double_t ymean = (fPos.Y() + Pos.Y()) / 2.; + Double_t zmean = (fPos.Z() + Pos.Z()) / 2.; + Double_t deltaTrackLength = gMC->TrackLength() - fLength; + AddHit(fTrackID, + det_uniqueId, + TVector3(xmean, ymean, zmean), + TVector3(fMom.Px(), fMom.Py(), fMom.Pz()), + fTime, + deltaTrackLength, + fELoss, + pdgCode); + + // Increment number of strawtubes det points in TParticle + /*ShipStack* stack = (ShipStack*) gMC->GetStack(); + stack->AddPoint(kStraw);*/ + } return kTRUE; } -prestrawdetectorPoint* prestrawdetector::AddHit(Int_t trackID, Int_t detID, - TVector3 pos, TVector3 mom, - Double_t time, Double_t length, - Double_t eLoss, Int_t pdgCode) +prestrawdetectorPoint* prestrawdetector::AddHit(Int_t trackID, + Int_t detID, + TVector3 pos, + TVector3 mom, + Double_t time, + Double_t length, + Double_t eLoss, + Int_t pdgCode) { - TClonesArray& clref = *fprestrawdetectorPointCollection; - Int_t size = clref.GetEntriesFast(); - return new(clref[size]) prestrawdetectorPoint(trackID, detID, pos, mom, - time, length, eLoss, pdgCode); + TClonesArray& clref = *fprestrawdetectorPointCollection; + Int_t size = clref.GetEntriesFast(); + return new (clref[size]) prestrawdetectorPoint(trackID, detID, pos, mom, time, length, eLoss, pdgCode); } void prestrawdetector::Register() { - /** This will create a branch in the output tree called - strawtubesPoint, setting the last parameter to kFALSE means: - this collection will not be written to the file, it will exist - only during the simulation. - */ - FairRootManager::Instance()->Register("prestrawdetectorPoint", "prestrawdetector", - fprestrawdetectorPointCollection, kTRUE); + /** This will create a branch in the output tree called + strawtubesPoint, setting the last parameter to kFALSE means: + this collection will not be written to the file, it will exist + only during the simulation. + */ + FairRootManager::Instance()->Register( + "prestrawdetectorPoint", "prestrawdetector", fprestrawdetectorPointCollection, kTRUE); } TClonesArray* prestrawdetector::GetCollection(Int_t iColl) const { - if (iColl == 0) { return fprestrawdetectorPointCollection; } - else { return NULL; } + if (iColl == 0) { + return fprestrawdetectorPointCollection; + } else { + return NULL; + } } void prestrawdetector::Reset() { - fprestrawdetectorPointCollection->Clear(); + fprestrawdetectorPointCollection->Clear(); } void prestrawdetector::SetZ(Double_t z0) { - fPSDz = z0; //! z-position of veto station + fPSDz = z0; //! z-position of veto station } void prestrawdetector::SetX(Double_t x) { - fPSDx = x; //! x-half-size + fPSDx = x; //! x-half-size } void prestrawdetector::SetY(Double_t y) { - fPSDy = y; //! y-half-size + fPSDy = y; //! y-half-size } void prestrawdetector::SetThickness(Double_t thickness) { - fPSD_Thickness = thickness; //! thickness + fPSD_Thickness = thickness; //! thickness } void prestrawdetector::ConstructGeometry() { - TGeoVolume *top = gGeoManager->GetTopVolume(); + TGeoVolume* top = gGeoManager->GetTopVolume(); InitMedium("air"); - TGeoMedium *air = gGeoManager->GetMedium("air"); + TGeoMedium* air = gGeoManager->GetMedium("air"); InitMedium("vacuum"); - TGeoMedium *vacuum = gGeoManager->GetMedium("vacuum"); + TGeoMedium* vacuum = gGeoManager->GetMedium("vacuum"); InitMedium("ShipSens"); - TGeoMedium *Se = gGeoManager->GetMedium("ShipSens"); + TGeoMedium* Se = gGeoManager->GetMedium("ShipSens"); InitMedium("aluminium"); - TGeoMedium *Al = gGeoManager->GetMedium("aluminium"); + TGeoMedium* Al = gGeoManager->GetMedium("aluminium"); InitMedium("mylar"); - TGeoMedium *mylar = gGeoManager->GetMedium("mylar"); + TGeoMedium* mylar = gGeoManager->GetMedium("mylar"); InitMedium("STTmix9010_2bar"); - TGeoMedium *sttmix9010_2bar = gGeoManager->GetMedium("STTmix9010_2bar"); + TGeoMedium* sttmix9010_2bar = gGeoManager->GetMedium("STTmix9010_2bar"); InitMedium("tungsten"); - TGeoMedium *tungsten = gGeoManager->GetMedium("tungsten"); - + TGeoMedium* tungsten = gGeoManager->GetMedium("tungsten"); gGeoManager->SetVisLevel(4); gGeoManager->SetTopVisible(); - Double_t eps=0.0001; + Double_t eps = 0.0001; TGeoBBox* myBox = new TGeoBBox("myBox", fPSDx, fPSDy, fPSD_Thickness); TGeoVolume* myVol = new TGeoVolume("myVol", myBox, vacuum); myVol->SetLineColor(kYellow); top->AddNode(myVol, 0, new TGeoTranslation(0, 0, fPSDz)); - AddSensitiveVolume(myVol); + AddSensitiveVolume(myVol); // Position it somewhere // top->AddNode(myVol, 1, new TGeoTranslation(0, 0, 500)); @@ -246,5 +254,5 @@ void prestrawdetector::ConstructGeometry() void prestrawdetector::EndOfEvent() { - fprestrawdetectorPointCollection->Clear(); + fprestrawdetectorPointCollection->Clear(); } diff --git a/prestrawdetector/prestrawdetector.h b/prestrawdetector/prestrawdetector.h index 718af354df..05409d2b95 100644 --- a/prestrawdetector/prestrawdetector.h +++ b/prestrawdetector/prestrawdetector.h @@ -3,9 +3,9 @@ #include "FairDetector.h" #include "FairVolume.h" -#include "TVector3.h" #include "TClonesArray.h" #include "TLorentzVector.h" +#include "TVector3.h" class prestrawdetectorPoint; class FairVolume; @@ -20,23 +20,27 @@ class prestrawdetector : public FairDetector virtual ~prestrawdetector(); /** Initialization of the detector is done here */ - virtual void Initialize(); //done + virtual void Initialize(); // done + + virtual Bool_t ProcessHits(FairVolume* v = 0); // done - virtual Bool_t ProcessHits( FairVolume* v=0);//done + prestrawdetectorPoint* AddHit(Int_t trackID, + Int_t detID, + TVector3 pos, + TVector3 mom, + Double_t time, + Double_t length, + Double_t eLoss, + Int_t pdgCode); - prestrawdetectorPoint* AddHit(Int_t trackID, Int_t detID, - TVector3 pos, TVector3 mom, - Double_t time, Double_t length, - Double_t eLoss, Int_t pdgCode); - /** Registers the produced collections in FAIRRootManager. */ - virtual void Register(); + virtual void Register(); /** Gets the produced collections */ - virtual TClonesArray* GetCollection(Int_t iColl) const ; //done + virtual TClonesArray* GetCollection(Int_t iColl) const; // done /** has to be called after each event to reset the containers */ - virtual void Reset(); + virtual void Reset(); void SetZ(Double_t z0); void SetX(Double_t x); @@ -44,39 +48,37 @@ class prestrawdetector : public FairDetector void SetThickness(Double_t thickness); void ConstructGeometry(); - - virtual void CopyClones( TClonesArray* cl1, TClonesArray* cl2 , - Int_t offset) {;} - virtual void SetSpecialPhysicsCuts() {;} - virtual void EndOfEvent(); - virtual void FinishPrimary() {;} - virtual void FinishRun() {;} - virtual void BeginPrimary() {;} - virtual void PostTrack() {;} - virtual void PreTrack() {;} - virtual void BeginEvent() {;} + virtual void CopyClones(TClonesArray* cl1, TClonesArray* cl2, Int_t offset) { ; } + virtual void SetSpecialPhysicsCuts() { ; } + virtual void EndOfEvent(); + virtual void FinishPrimary() { ; } + virtual void FinishRun() { ; } + virtual void BeginPrimary() { ; } + virtual void PostTrack() { ; } + virtual void PreTrack() { ; } + virtual void BeginEvent() { ; } private: - Double_t fPSDz; - Double_t fPSDx; //! x(half of the size) in cm - Double_t fPSDy; //! y(half of the size) in cm - Double_t fPSD_Thickness; //! thickness in cm - std::string fMedium; - Int_t fTrackID; //! track index - Int_t fVolumeID; //! volume id - TLorentzVector fPos; //! position at entrance - TLorentzVector fMom; //! momentum at entrance - Double_t fTime; //! time - Double_t fLength; //! length - Double_t fELoss; - TClonesArray* fprestrawdetectorPointCollection; // collection of hit points + Double_t fPSDz; + Double_t fPSDx; //! x(half of the size) in cm + Double_t fPSDy; //! y(half of the size) in cm + Double_t fPSD_Thickness; //! thickness in cm + std::string fMedium; + Int_t fTrackID; //! track index + Int_t fVolumeID; //! volume id + TLorentzVector fPos; //! position at entrance + TLorentzVector fMom; //! momentum at entrance + Double_t fTime; //! time + Double_t fLength; //! length + Double_t fELoss; + TClonesArray* fprestrawdetectorPointCollection; // collection of hit points prestrawdetector(const prestrawdetector&); prestrawdetector& operator=(const prestrawdetector&); Int_t InitMedium(const char* name); - ClassDef(prestrawdetector,2) + ClassDef(prestrawdetector, 2) }; #endif From 2ced8d18e3f174cbb9e070f38d24d72d48af67d5 Mon Sep 17 00:00:00 2001 From: AnalitikGnid Date: Wed, 3 Dec 2025 00:25:09 +0100 Subject: [PATCH 23/25] update --- macro/run_simScript.py | 42 ++++++++++++++++--- python/shipDigiReco.py | 2 +- shipgen/CMakeLists.txt | 2 +- shipgen/LinkDef.h | 2 +- ...{MyGenerator.cxx => PrestrawGenerator.cxx} | 20 ++++----- .../{MyGenerator.h => PrestrawGenerator.h} | 14 +++---- sstDecouplingTools/sst.json | 2 +- 7 files changed, 57 insertions(+), 27 deletions(-) rename shipgen/{MyGenerator.cxx => PrestrawGenerator.cxx} (78%) rename shipgen/{MyGenerator.h => PrestrawGenerator.h} (74%) diff --git a/macro/run_simScript.py b/macro/run_simScript.py index 9fd090aefc..5eece207dc 100755 --- a/macro/run_simScript.py +++ b/macro/run_simScript.py @@ -3,6 +3,10 @@ import sys import ROOT +import time +t0 = time.perf_counter() +c0 = time.process_time() + import shipunit as u import shipRoot_conf import rootUtils as ut @@ -307,6 +311,9 @@ # Parameter file name parFile=f"{options.outputDir}/ship.params.{tag}.root" +t1 = time.perf_counter() +c1 = time.process_time() +print(f"TIME 1 {t1-t0} and CPU {c1 - c0}") # In general, the following parts need not be touched # ======================================================================== @@ -320,12 +327,19 @@ run.SetSink(ROOT.FairRootFileSink(outFile)) # Output file run.SetUserConfig("g4Config.C") # user configuration file default g4Config.C rtdb = run.GetRuntimeDb() + +t2 = time.perf_counter() +c2 = time.process_time() +print(f"TIME 2 {t2-t1} and CPU {c2 - c1}") # -----Create geometry---------------------------------------------- # import shipMuShield_only as shipDet_conf # special use case for an attempt to convert active shielding geometry for use with FLUKA # import shipTarget_only as shipDet_conf import shipDet_conf modules = shipDet_conf.configure(run,ship_geo) +t3 = time.perf_counter() +c3 = time.process_time() +print(f"TIME 3 {t3-t2} and CPU {c3 - c2}") # -----Create PrimaryGenerator-------------------------------------- primGen = ROOT.FairPrimaryGenerator() if options.pythia8: @@ -416,7 +430,7 @@ if options.ttreegen: primGen.SetTarget(0.0, 0.0) print(f"Opening input file for ttree generator: {inputFile}") - TtreeGen = ROOT.MyGenerator() + TtreeGen = ROOT.PrestrawGenerator() TtreeGen.Init(inputFile, options.firstEvent) primGen.AddGenerator(TtreeGen) options.nEvents = TtreeGen.GetNevents() @@ -567,14 +581,21 @@ if options.evtgen_decayer: run.SetPythiaDecayer('DecayConfigTEvtGen.C') print('Using TEvtGenDecayer for J/psi and quarkonium decays with EvtGen') - +t4 = time.perf_counter() +c4 = time.process_time() +print(f"TIME 4 {t4-t3} and CPU {c4 - c3}") # -----Initialize simulation run------------------------------------ run.Init() +ta = time.perf_counter() +ca = time.process_time() +print(f"TIME a {ta-t4} and CPU {ca - c4}") if options.dryrun: # Early stop after setting up Pythia 8 sys.exit(0) gMC = ROOT.TVirtualMC.GetMC() fStack = gMC.GetStack() - +tb = time.perf_counter() +cb = time.process_time() +print(f"TIME b {tb-ta} and CPU {cb - ca}") # -----J/psi external decayer configuration handled in g4config.in------------------------------------ # VMC command /mcPhysics/setExtDecayerSelection J/psi forces external decayer usage EnergyCut = 10. * u.MeV if options.mudis else 100. * u.MeV @@ -601,7 +622,9 @@ trajFilter.SetEnergyCut(0., 400.*u.GeV) trajFilter.SetStorePrimaries(ROOT.kTRUE) trajFilter.SetStoreSecondaries(ROOT.kTRUE) - +tc = time.perf_counter() +cc = time.process_time() +print(f"TIME c {tc-tb} and CPU {cc - cb}") # The VMC sets the fields using the "/mcDet/setIsLocalMagField true" option in "gconfig/g4config.in" import geomGeant4 @@ -623,9 +646,14 @@ # Plot the field example #fieldMaker.plotField(1, ROOT.TVector3(-9000.0, 6000.0, 50.0), ROOT.TVector3(-300.0, 300.0, 6.0), 'Bzx.png') #fieldMaker.plotField(2, ROOT.TVector3(-9000.0, 6000.0, 50.0), ROOT.TVector3(-400.0, 400.0, 6.0), 'Bzy.png') - +t5 = time.perf_counter() +c5 = time.process_time() +print(f"TIME 5 {t5-tc} and CPU {c5-cc}") # -----Start run---------------------------------------------------- run.Run(options.nEvents) +t6 = time.perf_counter() +c6 = time.process_time() +print(f"TIME 6 {t6-t5} and CPU {c6 - c5}") # -----Runtime database--------------------------------------------- kParameterMerged = ROOT.kTRUE parOut = ROOT.FairParRootFileIo(kParameterMerged) @@ -640,7 +668,9 @@ import saveBasicParameters saveBasicParameters.execute(f"{options.outputDir}/geofile_full.{tag}.root",ship_geo) - +t7 = time.perf_counter() +c7 = time.process_time() +print(f"TIME 7 {t7-t6} and CPU {c7 - c6}") # checking for overlaps if options.check_overlaps: ROOT.gROOT.SetWebDisplay("off") # Workaround for https://github.com/root-project/root/issues/18881 diff --git a/python/shipDigiReco.py b/python/shipDigiReco.py index 8510d7e900..8ff8466251 100644 --- a/python/shipDigiReco.py +++ b/python/shipDigiReco.py @@ -109,7 +109,7 @@ def reconstruct(self): self.Vertexing.execute() def digitize(self): - self.sTree.t0 = self.random.Rndm()*1*u.microsecond + self.sTree.t0 = 0#self.random.Rndm()*1*u.microsecond self.header.SetEventTime( self.sTree.t0 ) self.header.SetRunId( self.sTree.MCEventHeader.GetRunID() ) self.header.SetMCEntryNumber( self.sTree.MCEventHeader.GetEventID() ) # counts from 1 diff --git a/shipgen/CMakeLists.txt b/shipgen/CMakeLists.txt index 86d607cffc..99ee18949f 100644 --- a/shipgen/CMakeLists.txt +++ b/shipgen/CMakeLists.txt @@ -78,7 +78,7 @@ set(SRCS TEvtGenDecayer.cxx BeamSmearingUtils.cxx MeanMaterialBudget.cxx - MyGenerator.cxx) + PrestrawGenerator.cxx) set(LINKDEF LinkDef.h) set(LIBRARY_NAME ShipGen) diff --git a/shipgen/LinkDef.h b/shipgen/LinkDef.h index 97cb940099..69e2ead145 100644 --- a/shipgen/LinkDef.h +++ b/shipgen/LinkDef.h @@ -13,7 +13,7 @@ #pragma link C++ class DPPythia8Generator+; #pragma link C++ class GenieGenerator+; #pragma link C++ class NtupleGenerator+; -#pragma link C++ class MyGenerator+; +#pragma link C++ class PrestrawGenerator+; #pragma link C++ class MuonBackGenerator+; #pragma link C++ class CosmicsGenerator+; #pragma link C++ class MuDISGenerator+; diff --git a/shipgen/MyGenerator.cxx b/shipgen/PrestrawGenerator.cxx similarity index 78% rename from shipgen/MyGenerator.cxx rename to shipgen/PrestrawGenerator.cxx index e3a4bbc67e..6229093bdd 100644 --- a/shipgen/MyGenerator.cxx +++ b/shipgen/PrestrawGenerator.cxx @@ -2,7 +2,7 @@ #include "TROOT.h" #include "TFile.h" #include "FairPrimaryGenerator.h" -#include "MyGenerator.h" +#include "PrestrawGenerator.h" #include "TDatabasePDG.h" // for TDatabasePDG #include "TMath.h" // for Sqrt @@ -11,18 +11,18 @@ using std::endl; // read events from ntuples produced // ----- Default constructor ------------------------------------------- -MyGenerator::MyGenerator() {} +PrestrawGenerator::PrestrawGenerator() {} // ------------------------------------------------------------------------- // ----- Default constructor ------------------------------------------- -Bool_t MyGenerator::Init(const char* fileName) { +Bool_t PrestrawGenerator::Init(const char* fileName) { return Init(fileName, 0); } // ----- Default constructor ------------------------------------------- -Bool_t MyGenerator::Init(const char* fileName, const int firstEvent) { - cout << "Info MyGenerator: Opening input file " << fileName << endl; +Bool_t PrestrawGenerator::Init(const char* fileName, const int firstEvent) { + cout << "Info PrestrawGenerator: Opening input file " << fileName << endl; fInputFile = new TFile(fileName); if (fInputFile->IsZombie()) { - cout << "-E MyGenerator: Error opening the Signal file" << fileName << endl; + cout << "-E PrestrawGenerator: Error opening the Signal file" << fileName << endl; } fTree = (TTree *)fInputFile->Get("mytree"); fNevents = fTree->GetEntries(); @@ -40,9 +40,9 @@ Bool_t MyGenerator::Init(const char* fileName, const int firstEvent) { // ----- Destructor ---------------------------------------------------- -MyGenerator::~MyGenerator() +PrestrawGenerator::~PrestrawGenerator() { - // cout << "destroy My" << endl; + // cout << "destroy prestraw" << endl; fInputFile->Close(); fInputFile->Delete(); delete fInputFile; @@ -50,7 +50,7 @@ MyGenerator::~MyGenerator() // ------------------------------------------------------------------------- // ----- Passing the event --------------------------------------------- -Bool_t MyGenerator::ReadEvent(FairPrimaryGenerator* cpg) +Bool_t PrestrawGenerator::ReadEvent(FairPrimaryGenerator* cpg) { cout < Date: Wed, 3 Dec 2025 04:33:37 +0000 Subject: [PATCH 24/25] style: pre-commit fixes --- macro/run_simScript.py | 2 +- macro/update_config.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/macro/run_simScript.py b/macro/run_simScript.py index 5eece207dc..66f890a2a8 100755 --- a/macro/run_simScript.py +++ b/macro/run_simScript.py @@ -3,7 +3,7 @@ import sys import ROOT -import time +import time t0 = time.perf_counter() c0 = time.process_time() diff --git a/macro/update_config.py b/macro/update_config.py index 0601addab0..25aa940b8a 100644 --- a/macro/update_config.py +++ b/macro/update_config.py @@ -2,7 +2,7 @@ from ShipGeoConfig import AttrDict def update_config(global_config, my_config, name): - with open(my_config, 'r') as f: + with open(my_config) as f: data = json.load(f) new_values = AttrDict(data[name]) if name == 'geometry': From 91ff7fdf74e2a07e77b06482ebf995d41e92314a Mon Sep 17 00:00:00 2001 From: AnalitikGnid Date: Wed, 3 Dec 2025 05:34:54 +0100 Subject: [PATCH 25/25] bug fixes --- macro/ShipReco.py | 6 +- macro/run_simScript.py | 141 ++++++++++++++++++----------------------- python/shipDet_conf.py | 26 +++++++- python/shipDigiReco.py | 14 ++-- 4 files changed, 94 insertions(+), 93 deletions(-) diff --git a/macro/ShipReco.py b/macro/ShipReco.py index 89a4ad7e15..e985d79368 100644 --- a/macro/ShipReco.py +++ b/macro/ShipReco.py @@ -44,6 +44,7 @@ def mem_monitor(): , required=False, choices=['FH','AR','TemplateMatching'], default='') parser.add_argument("-dy", dest="dy", help="Max height of tank", required=False, default=None,type=int) parser.add_argument("--Debug", dest="Debug", help="Switch on debugging", required=False, action="store_true") +parser.add_argument("--digiOnly", dest="digionly", help="only digitization", default=False, required=False, action="store_true") options = parser.parse_args() vertexing = not options.noVertexing @@ -131,8 +132,9 @@ def mem_monitor(): print('event ', global_variables.iEvent) rc = SHiP.sTree.GetEvent(global_variables.iEvent) SHiP.digitize() - SHiP.reconstruct() - SHiP.recoTree.Fill() + if not options.digionly: + SHiP.reconstruct() + SHiP.recoTree.Fill() # memory monitoring # mem_monitor() # end loop over events diff --git a/macro/run_simScript.py b/macro/run_simScript.py index 5eece207dc..a424eed865 100755 --- a/macro/run_simScript.py +++ b/macro/run_simScript.py @@ -1,20 +1,20 @@ #!/usr/bin/env python +# SPDX-License-Identifier: LGPL-3.0-or-later +# SPDX-FileCopyrightText: Copyright CERN for the benefit of the SHiP Collaboration + import os import sys import ROOT -import time -t0 = time.perf_counter() -c0 = time.process_time() - import shipunit as u import shipRoot_conf import rootUtils as ut -from ShipGeoConfig import ConfigRegistry +import geometry_config from argparse import ArgumentParser from array import array from backports import tdirectory634 DownScaleDiMuon = False + from update_config import update_config # Default HNL parameters theHNLMass = 1.0*u.GeV @@ -40,21 +40,6 @@ inputFile = "$EOSSHIP/eos/experiment/ship/data/Charm/Cascade-parp16-MSTP82-1-MSEL4-978Bpot.root" defaultInputFile = True -globalDesigns = { - '2023' : { - 'dy' : 6., - 'caloDesign' : 3, - 'strawDesign' : 10 - }, - '2025' : { - 'dy' : 6., - 'ds' : 8, - 'caloDesign' : 2, - 'strawDesign' : 10 - }, -} -default = '2025' - parser = ArgumentParser() group = parser.add_mutually_exclusive_group() @@ -62,8 +47,8 @@ parser.add_argument("--Pythia6", dest="pythia6", help="Use Pythia6", action="store_true") parser.add_argument("--Pythia8", dest="pythia8", help="Use Pythia8", action="store_true") parser.add_argument("--EvtGenDecayer", dest="evtgen_decayer", help="Use TEvtGenDecayer for J/psi and other quarkonium decays", action="store_true") -# === PG subcommand === subparsers = parser.add_subparsers(dest="command", help="Which mode to run") +# === PG subcommand === pg_parser = subparsers.add_parser("PG", help="Use Particle Gun") pg_parser.add_argument( @@ -103,9 +88,19 @@ "--Dy", dest="Dy", type=float, help="size of the full uniform spread of PG ypos: (Vy - Dy/2, Vy + Dy/2)" ) -# === Enf of PG commands === +# === End of PG commands === +# === Genie subcommand === +genie_parser = subparsers.add_parser("Genie", help="Genie for reading and processing neutrino interactions") +genie_parser.add_argument( + "--z_start_nu", dest="z_start_nu", default=2844.2850, type=float, + help="Genie neutrino start z start position (default=2844.2850 cm)" +) +genie_parser.add_argument( + "--z_end_nu", dest="z_end_nu", default=3180.4350, type=float, + help="Genie neutrino end z position (default=3180.4350 cm)" +) +# === End of Genie subcommand === parser.add_argument("-A", help="b: signal from b, c: from c (default), bc: from Bc, or inclusive", default='c') -parser.add_argument("--Genie", dest="genie", help="Genie for reading and processing neutrino interactions", action="store_true") parser.add_argument("--NuRadio", dest="nuradio", help="misuse GenieGenerator for neutrino radiography and geometry timing test", action="store_true") parser.add_argument("--Ntuple", dest="ntuple", help="Use ntuple as input", action="store_true") parser.add_argument("--MuonBack", dest="muonback", help="Generate events from muon background file, --Cosmics=0 for cosmic generator data", action="store_true") @@ -134,15 +129,10 @@ group.add_argument("-f", dest="inputFile", help="Input file if not default file", default=False) parser.add_argument("-g", dest="geofile", help="geofile for muon shield geometry, for experts only", default=None) parser.add_argument("-o", "--output", dest="outputDir", help="Output directory", default=".") -parser.add_argument("-Y", dest="dy", help="max height of vacuum tank", default=globalDesigns[default]['dy']) -parser.add_argument("--caloDesign", - help="0=ECAL/HCAL TP 2=splitCal 3=ECAL/ passive HCAL", - default=globalDesigns[default]['caloDesign'], - type=int, - choices=[0,2,3]) +parser.add_argument("-Y", dest="dy", help="max height of vacuum tank", default=6.0, type=float) parser.add_argument("--strawDesign", help="Tracker station frame material: 4=aluminium; 10=steel (default)", - default=globalDesigns[default]['strawDesign'], + default=10, type=int, choices=[4,10]) parser.add_argument("-F", dest="deepCopy", help="default = False: copy only stable particles to stack, except for HNL events", action="store_true") @@ -176,7 +166,7 @@ parser.add_argument("--SND", dest="SND", help="Activate SND.", action='store_true') parser.add_argument( "--SND_design", - help="Choose SND design(s) among [1,2,...] or 'all' to enable all. 1: EmulsionTarget, 2: MTC", + help="Choose SND design(s) among [1,2,...] or 'all' to enable all. 1: EmulsionTarget, 2: MTC + SiliconTarget", nargs='+', default=[2], ) @@ -192,9 +182,6 @@ options = parser.parse_args() if options.json: update_config(options, os.path.expandvars('$FAIRSHIP/sstDecouplingTools/sst.json'), 'options') - - - # Handle SND_design: allow 'all' (case-insensitive) or list of ints available_snd_designs = [1, 2] # Extend this list as new designs are added if any(str(x).lower() == 'all' for x in options.SND_design): @@ -242,7 +229,7 @@ print("cannot have HNL and SUSY or DP at the same time, abort") sys.exit(2) -if (options.genie or options.nuradio) and defaultInputFile: +if (options.command == "Genie" or options.nuradio) and defaultInputFile: inputFile = "$EOSSHIP/eos/experiment/ship/data/GenieEvents/genie-nu_mu.root" if options.mudis and defaultInputFile: print('input file required if simEngine = muonDIS') @@ -266,10 +253,8 @@ ROOT.gInterpreter.ProcessLine('fair::Logger::SetConsoleSeverity("debug1");') elif options.debug == 3: ROOT.gInterpreter.ProcessLine('fair::Logger::SetConsoleSeverity("debug2");') -ship_geo = ConfigRegistry.loadpy( - "$FAIRSHIP/geometry/geometry_config.py", +ship_geo = geometry_config.create_config( Yheight=options.dy, - CaloDesign=options.caloDesign, strawDesign=options.strawDesign, muShieldGeo=options.geofile, shieldName=options.shieldName, @@ -289,8 +274,11 @@ # Output file name, add dy to be able to setup geometry with ambiguities. if options.command == "PG": tag = f"PG_{options.pID}-{mcEngine}" +elif options.command == "Genie": + tag = f"Genie-{mcEngine}" + simEngine = "Genie" else: - for g in ["pythia8", "evtcalc", "pythia6", "genie", "nuradio", "ntuple", "muonback", "mudis", "fixedTarget", "cosmics", "ttreegen"]: + for g in ["pythia8", "evtcalc", "pythia6", "nuradio", "ntuple", "muonback", "mudis", "fixedTarget", "cosmics", "ttreegen"]: if getattr(options, g): simEngine = g.capitalize() break @@ -311,9 +299,6 @@ # Parameter file name parFile=f"{options.outputDir}/ship.params.{tag}.root" -t1 = time.perf_counter() -c1 = time.process_time() -print(f"TIME 1 {t1-t0} and CPU {c1 - c0}") # In general, the following parts need not be touched # ======================================================================== @@ -327,24 +312,16 @@ run.SetSink(ROOT.FairRootFileSink(outFile)) # Output file run.SetUserConfig("g4Config.C") # user configuration file default g4Config.C rtdb = run.GetRuntimeDb() - -t2 = time.perf_counter() -c2 = time.process_time() -print(f"TIME 2 {t2-t1} and CPU {c2 - c1}") # -----Create geometry---------------------------------------------- # import shipMuShield_only as shipDet_conf # special use case for an attempt to convert active shielding geometry for use with FLUKA # import shipTarget_only as shipDet_conf import shipDet_conf modules = shipDet_conf.configure(run,ship_geo) -t3 = time.perf_counter() -c3 = time.process_time() -print(f"TIME 3 {t3-t2} and CPU {c3 - c2}") # -----Create PrimaryGenerator-------------------------------------- primGen = ROOT.FairPrimaryGenerator() if options.pythia8: primGen.SetTarget(ship_geo.target.z0, 0.) - # -----Pythia8-------------------------------------- if HNL or options.RPVSUSY: P8gen = ROOT.HNLPythia8Generator() @@ -398,9 +375,9 @@ P8gen.UseExternalFile(inputFile, options.firstEvent) # Use geometry constants instead of fragile TGeo navigation P8gen.SetTargetCoordinates(ship_geo.target.z0, ship_geo.target.z0 + ship_geo.target.length) - # pion on proton 500GeV - # P8gen.SetMom(500.*u.GeV) - # P8gen.SetId(-211) +# pion on proton 500GeV +# P8gen.SetMom(500.*u.GeV) +# P8gen.SetId(-211) primGen.AddGenerator(P8gen) if options.fixedTarget: HNL = False @@ -424,7 +401,6 @@ P6gen.SetMom(50.*u.GeV) P6gen.SetTarget("gamma/mu+","n0") # default "gamma/mu-","p+" primGen.AddGenerator(P6gen) - # -----TTree generator-------------------------------------- if options.ttreegen: @@ -482,13 +458,13 @@ options.nEvents = min(options.nEvents,DISgen.GetNevents()) print('Generate ',options.nEvents,' with DIS input', ' first event',options.firstEvent) # -----Neutrino Background------------------------ -if options.genie: +if options.command == "Genie": # Genie ut.checkFileExists(inputFile) primGen.SetTarget(0., 0.) # do not interfere with GenieGenerator Geniegen = ROOT.GenieGenerator() Geniegen.Init(inputFile,options.firstEvent) - Geniegen.SetPositions(ship_geo.target.z0, ship_geo.tauMudet.zMudetC-5*u.m, ship_geo.TrackStation2.z) + Geniegen.SetPositions(ship_geo.target.z0, options.z_start_nu, options.z_end_nu) primGen.AddGenerator(Geniegen) options.nEvents = min(options.nEvents,Geniegen.GetNevents()) run.SetPythiaDecayer("DecayConfigNuAge.C") @@ -581,21 +557,14 @@ if options.evtgen_decayer: run.SetPythiaDecayer('DecayConfigTEvtGen.C') print('Using TEvtGenDecayer for J/psi and quarkonium decays with EvtGen') -t4 = time.perf_counter() -c4 = time.process_time() -print(f"TIME 4 {t4-t3} and CPU {c4 - c3}") + # -----Initialize simulation run------------------------------------ run.Init() -ta = time.perf_counter() -ca = time.process_time() -print(f"TIME a {ta-t4} and CPU {ca - c4}") if options.dryrun: # Early stop after setting up Pythia 8 sys.exit(0) gMC = ROOT.TVirtualMC.GetMC() fStack = gMC.GetStack() -tb = time.perf_counter() -cb = time.process_time() -print(f"TIME b {tb-ta} and CPU {cb - ca}") + # -----J/psi external decayer configuration handled in g4config.in------------------------------------ # VMC command /mcPhysics/setExtDecayerSelection J/psi forces external decayer usage EnergyCut = 10. * u.MeV if options.mudis else 100. * u.MeV @@ -622,9 +591,7 @@ trajFilter.SetEnergyCut(0., 400.*u.GeV) trajFilter.SetStorePrimaries(ROOT.kTRUE) trajFilter.SetStoreSecondaries(ROOT.kTRUE) -tc = time.perf_counter() -cc = time.process_time() -print(f"TIME c {tc-tb} and CPU {cc - cb}") + # The VMC sets the fields using the "/mcDet/setIsLocalMagField true" option in "gconfig/g4config.in" import geomGeant4 @@ -642,18 +609,13 @@ if options.print_fields: geomGeant4.printVMCFields() geomGeant4.printWeightsandFields(onlyWithField = True,\ - exclude=['DecayVolume','Tr1','Tr2','Tr3','Tr4','Veto','Ecal','Hcal','MuonDetector','SplitCal']) + exclude=['DecayVolume','Tr1','Tr2','Tr3','Tr4','Veto','MuonDetector','SplitCal']) # Plot the field example #fieldMaker.plotField(1, ROOT.TVector3(-9000.0, 6000.0, 50.0), ROOT.TVector3(-300.0, 300.0, 6.0), 'Bzx.png') #fieldMaker.plotField(2, ROOT.TVector3(-9000.0, 6000.0, 50.0), ROOT.TVector3(-400.0, 400.0, 6.0), 'Bzy.png') -t5 = time.perf_counter() -c5 = time.process_time() -print(f"TIME 5 {t5-tc} and CPU {c5-cc}") + # -----Start run---------------------------------------------------- run.Run(options.nEvents) -t6 = time.perf_counter() -c6 = time.process_time() -print(f"TIME 6 {t6-t5} and CPU {c6 - c5}") # -----Runtime database--------------------------------------------- kParameterMerged = ROOT.kTRUE parOut = ROOT.FairParRootFileIo(kParameterMerged) @@ -668,9 +630,7 @@ import saveBasicParameters saveBasicParameters.execute(f"{options.outputDir}/geofile_full.{tag}.root",ship_geo) -t7 = time.perf_counter() -c7 = time.process_time() -print(f"TIME 7 {t7-t6} and CPU {c7 - c6}") + # checking for overlaps if options.check_overlaps: ROOT.gROOT.SetWebDisplay("off") # Workaround for https://github.com/root-project/root/issues/18881 @@ -736,13 +696,11 @@ branches.Add(ROOT.TObjString('TTPoint')) branches.Add(ROOT.TObjString('ScoringPoint')) branches.Add(ROOT.TObjString('strawtubesPoint')) - branches.Add(ROOT.TObjString('EcalPoint')) - branches.Add(ROOT.TObjString('sEcalPointLite')) - branches.Add(ROOT.TObjString('smuonPoint')) branches.Add(ROOT.TObjString('TimeDetPoint')) branches.Add(ROOT.TObjString('MCEventHeader')) branches.Add(ROOT.TObjString('UpstreamTaggerPoint')) branches.Add(ROOT.TObjString('MTCdetPoint')) + branches.Add(ROOT.TObjString('SiliconTargetPoint')) branches.Add(ROOT.TObjString('sGeoTracks')) sTree.AutoSave() @@ -785,6 +743,27 @@ os.replace(temp_filename, outFile) print("Successfully added DISCrossSection to the output file:", outFile) +if options.command == "Genie": + # breakpoint() + # Copy Genie (gst TTree) information to the output file + f_input = ROOT.TFile.Open(inputFile, "READ") + print("check") + gst = f_input.gst + + selection_string = "(Entry$ >= "+str(options.firstEvent)+")" + if (options.firstEvent + options.nEvents) < gst.GetEntries() : + selection_string += "&&(Entry$ < "+str(options.firstEvent + options.nEvents)+")" + + # Reopen output file + f_output = ROOT.TFile.Open(outFile, "UPDATE") + + # Copy only the events used in this run + gst_copy = gst.CopyTree(selection_string) + gst_copy.Write() + + f_input.Close() + f_output.Close() + # ------------------------------------------------------------------------ import checkMagFields diff --git a/python/shipDet_conf.py b/python/shipDet_conf.py index 8738e84f59..6c2b5493b2 100644 --- a/python/shipDet_conf.py +++ b/python/shipDet_conf.py @@ -257,12 +257,22 @@ def configure_strawtubes(yaml_file, ship_geo): ship_geo.strawtubesDigi.v_drift, ship_geo.strawtubesDigi.sigma_spatial, ) - Prestrawdetector = ROOT.prestrawdetector('Prestrawdetector', True) - Prestrawdetector.SetZ(ship_geo.psd) - detectorList.append(Prestrawdetector) detectorList.append(strawtubes) +def configure_prestrawdetector(yaml_file, ship_geo): + with open(yaml_file) as file: + config = yaml.safe_load(file) + + ship_geo.prestrawdetector_geo = AttrDict(config['PSD']) + + Prestrawdetector = ROOT.prestrawdetector(ship_geo.prestrawdetector_geo.material) + Prestrawdetector.SetZ(ship_geo.prestrawdetector_geo.z_position) + Prestrawdetector.SetX(ship_geo.prestrawdetector_geo.width) + Prestrawdetector.SetY(ship_geo.prestrawdetector_geo.height) + Prestrawdetector.SetThickness(ship_geo.prestrawdetector_geo.wall_thickness) + detectorList.append(Prestrawdetector) + def configure(run, ship_geo): # ---- for backward compatibility ---- if not hasattr(ship_geo, "DecayVolumeMedium"): @@ -395,6 +405,13 @@ def configure(run, ship_geo): ship_geo, ) + if hasattr(ship_geo, "PSD"): + if ship_geo.PSD: + configure_prestrawdetector( + os.path.join(os.environ["FAIRSHIP"], "geometry", "prestrawdetector_config.yaml"), + ship_geo, + ) + if ship_geo.EcalOption == 2: # splitCal with pointing information SplitCal = ROOT.splitcal("SplitCal", ROOT.kTRUE) x = ship_geo.SplitCal @@ -484,6 +501,9 @@ def configure(run, ship_geo): run.SetField(fMagField) exclusionList = [] + if ship_geo.decouple: + exclusionList = ["Muon","Ecal","Hcal","TargetTrackers","NuTauTarget","HighPrecisionTrackers",\ + "Veto","MuonShield","TargetStation","NuTauMudet","EmuMagnet", "TimeDet", "UpstreamTagger"] # exclusionList = ["strawtubes","TargetTrackers","NuTauTarget",\ # "SiliconTarget","Veto","Magnet","MuonShield","TargetStation", "TimeDet", "UpstreamTagger"] diff --git a/python/shipDigiReco.py b/python/shipDigiReco.py index 8ff8466251..6aee74ea51 100644 --- a/python/shipDigiReco.py +++ b/python/shipDigiReco.py @@ -114,15 +114,15 @@ def digitize(self): self.header.SetRunId( self.sTree.MCEventHeader.GetRunID() ) self.header.SetMCEntryNumber( self.sTree.MCEventHeader.GetEventID() ) # counts from 1 self.eventHeader.Fill() - self.digiSBT.process() + #self.digiSBT.process() self.strawtubes.process() - self.timeDetector.process() - self.upstreamTaggerDetector.process() + #self.timeDetector.process() + #self.upstreamTaggerDetector.process() # adding digitization of SND/MTC - if self.sTree.GetBranch("MTCDetPoint"): - self.digiMTC.process() - if self.sTree.GetBranch("splitcalPoint"): - self.splitcalDetector.process() + #if self.sTree.GetBranch("MTCDetPoint"): + #self.digiMTC.process() + #if self.sTree.GetBranch("splitcalPoint"): + #self.splitcalDetector.process() def findTracks(self): hitPosLists = {}