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/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 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 37ad2c276a..a424eed865 100755 --- a/macro/run_simScript.py +++ b/macro/run_simScript.py @@ -15,6 +15,7 @@ from backports import tdirectory634 DownScaleDiMuon = False +from update_config import update_config # Default HNL parameters theHNLMass = 1.0*u.GeV theProductionCouplings = theDecayCouplings = None @@ -172,9 +173,15 @@ 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): @@ -256,6 +263,13 @@ 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": @@ -264,7 +278,7 @@ 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", "nuradio", "ntuple", "muonback", "mudis", "fixedTarget", "cosmics", "ttreegen"]: if getattr(options, g): simEngine = g.capitalize() break @@ -387,7 +401,15 @@ P6gen.SetMom(50.*u.GeV) 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.PrestrawGenerator() + TtreeGen.Init(inputFile, options.firstEvent) + primGen.AddGenerator(TtreeGen) + options.nEvents = TtreeGen.GetNevents() # -----EvtCalc-------------------------------------- if options.evtcalc: primGen.SetTarget(0.0, 0.0) diff --git a/macro/update_config.py b/macro/update_config.py new file mode 100644 index 0000000000..25aa940b8a --- /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) 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/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..3e848ee79b --- /dev/null +++ b/prestrawdetector/prestrawdetector.cxx @@ -0,0 +1,258 @@ +#include "prestrawdetector.h" + +#include "FairGeoBuilder.h" +#include "FairGeoInterface.h" +#include "FairGeoLoader.h" +#include "FairGeoMedia.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 "TClonesArray.h" +#include "TGeoBBox.h" +#include "TGeoCompositeShape.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) + , 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(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; + } +} + +void prestrawdetector::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); +} + +Bool_t prestrawdetector::ProcessHits(FairVolume* vol) +{ + 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);*/ + } + 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(); + 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); +} + +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::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(); + InitMedium("air"); + TGeoMedium* air = gGeoManager->GetMedium("air"); + InitMedium("vacuum"); + TGeoMedium* vacuum = gGeoManager->GetMedium("vacuum"); + + 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", 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); + + // 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..05409d2b95 --- /dev/null +++ b/prestrawdetector/prestrawdetector.h @@ -0,0 +1,84 @@ +#ifndef PRESTRAWDETECOR_H +#define PRESTRAWDETECOR_H + +#include "FairDetector.h" +#include "FairVolume.h" +#include "TClonesArray.h" +#include "TLorentzVector.h" +#include "TVector3.h" + +class prestrawdetectorPoint; +class FairVolume; +class TClonesArray; + +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 + + 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 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() { ; } + + 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 + + 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..383452c409 --- /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..88ba5ff6b7 --- /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 diff --git a/python/shipDet_conf.py b/python/shipDet_conf.py index 160c95d59a..6c2b5493b2 100644 --- a/python/shipDet_conf.py +++ b/python/shipDet_conf.py @@ -257,10 +257,22 @@ def configure_strawtubes(yaml_file, ship_geo): ship_geo.strawtubesDigi.v_drift, ship_geo.strawtubesDigi.sigma_spatial, ) - 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"): @@ -393,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 @@ -482,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 8510d7e900..6aee74ea51 100644 --- a/python/shipDigiReco.py +++ b/python/shipDigiReco.py @@ -109,20 +109,20 @@ 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 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 = {} diff --git a/shipgen/CMakeLists.txt b/shipgen/CMakeLists.txt index 63a6a7e6d5..99ee18949f 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 + PrestrawGenerator.cxx) set(LINKDEF LinkDef.h) set(LIBRARY_NAME ShipGen) 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) diff --git a/shipgen/LinkDef.h b/shipgen/LinkDef.h index ed34b9f276..69e2ead145 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 PrestrawGenerator+; #pragma link C++ class MuonBackGenerator+; #pragma link C++ class CosmicsGenerator+; #pragma link C++ class MuDISGenerator+; diff --git a/shipgen/PrestrawGenerator.cxx b/shipgen/PrestrawGenerator.cxx new file mode 100644 index 0000000000..6229093bdd --- /dev/null +++ b/shipgen/PrestrawGenerator.cxx @@ -0,0 +1,78 @@ +#include +#include "TROOT.h" +#include "TFile.h" +#include "FairPrimaryGenerator.h" +#include "PrestrawGenerator.h" +#include "TDatabasePDG.h" // for TDatabasePDG +#include "TMath.h" // for Sqrt + +using std::cout; +using std::endl; +// read events from ntuples produced + +// ----- Default constructor ------------------------------------------- +PrestrawGenerator::PrestrawGenerator() {} +// ------------------------------------------------------------------------- +// ----- Default constructor ------------------------------------------- +Bool_t PrestrawGenerator::Init(const char* fileName) { + return Init(fileName, 0); +} +// ----- Default constructor ------------------------------------------- +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 PrestrawGenerator: 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); + return kTRUE; +} +// ------------------------------------------------------------------------- + + +// ----- Destructor ---------------------------------------------------- +PrestrawGenerator::~PrestrawGenerator() +{ + // cout << "destroy prestraw" << endl; + fInputFile->Close(); + fInputFile->Delete(); + delete fInputFile; +} +// ------------------------------------------------------------------------- + +// ----- Passing the event --------------------------------------------- +Bool_t PrestrawGenerator::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*px+py*py+pz*pz+ mass*mass ); + cpg->AddTrack(id, px, py, pz, vx, vy, vz); + + return kTRUE; +} + +// ------------------------------------------------------------------------- +Int_t PrestrawGenerator::GetNevents() +{ + return fNevents; +} diff --git a/shipgen/PrestrawGenerator.h b/shipgen/PrestrawGenerator.h new file mode 100644 index 0000000000..d9ce2092fe --- /dev/null +++ b/shipgen/PrestrawGenerator.h @@ -0,0 +1,40 @@ +#ifndef PRESTRAWGENERATOR_H +#define PRESTRAWGENERATOR_H 1 + +#include "TROOT.h" +#include "FairGenerator.h" +#include "TTree.h" // for TTree +#include "FairLogger.h" // for FairLogger, MESSAGE_ORIGIN + +class FairPrimaryGenerator; + +class PrestrawGenerator : public FairGenerator +{ + public: + + /** default constructor **/ + PrestrawGenerator(); + + /** destructor **/ + virtual ~PrestrawGenerator(); + + /** public method ReadEvent **/ + Bool_t ReadEvent(FairPrimaryGenerator*); + virtual Bool_t Init(const char*, int); //! + virtual Bool_t Init(const char*); //! + Int_t GetNevents(); + private: + + protected: + + Double_t px, py, pz, vx, vy, vz; + TFile* fInputFile; + TTree* fTree; + FairLogger* fLogger; //! don't make it persistent, magic ROOT command + int fNevents; + int fn; + Int_t id; + ClassDef(PrestrawGenerator,1); +}; + +#endif /* !PRESTRAWGENERATOR_H */ diff --git a/sstDecouplingTools/README.md b/sstDecouplingTools/README.md new file mode 100644 index 0000000000..314a5524c1 --- /dev/null +++ b/sstDecouplingTools/README.md @@ -0,0 +1,30 @@ +### Geometry tuning + +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` + +### SST decoupling + +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 +``` + +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 +``` 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..f796b3cb6c --- /dev/null +++ b/sstDecouplingTools/produce_ttree.C @@ -0,0 +1,86 @@ +#include "FairRootManager.h" +#include "ShipMCTrack.h" +#include "TLorentzVector.h" + +#include +#include +#include +#include +#include + +// Include class definition for strawtubesPoint if available +#include "prestrawdetectorPoint.h" +#include "strawtubesPoint.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..2926f3f0d5 --- /dev/null +++ b/sstDecouplingTools/sst.csv @@ -0,0 +1,9 @@ +key,value +options.theSeed,1234 +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 +ship_geo.psd,2588.0 diff --git a/sstDecouplingTools/sst.json b/sstDecouplingTools/sst.json new file mode 100644 index 0000000000..1f6a73f878 --- /dev/null +++ b/sstDecouplingTools/sst.json @@ -0,0 +1,11 @@ +{ + "geometry": { + "positions": [8407.0, 8607.0, 9307.0, 9507.0], + "angles": 4.57, + "magnetic_strength": 0.45 + }, + "options": { + "theSeed": 42, + "nEvents": 100 + } +}