From 5a580f6ab2f7aeb645c5ae2705a0b161b8f7d9a9 Mon Sep 17 00:00:00 2001 From: anupama-reghunath Date: Wed, 21 Aug 2024 23:13:19 +0200 Subject: [PATCH 01/10] Timing Implementation of MuDIS Added array for softTracks for MuDIS --- macro/run_simScript.py | 81 ++++++++++--------- shipgen/MuDISGenerator.cxx | 154 ++++++++++++++++++++++++++----------- shipgen/MuDISGenerator.h | 1 + 3 files changed, 157 insertions(+), 79 deletions(-) diff --git a/macro/run_simScript.py b/macro/run_simScript.py index fb11e32304..6e5eb4bed8 100755 --- a/macro/run_simScript.py +++ b/macro/run_simScript.py @@ -1,4 +1,6 @@ #!/usr/bin/env python +from __future__ import print_function +from __future__ import division import os import sys import ROOT @@ -31,7 +33,7 @@ MCTracksWithEnergyCutOnly = True # copy particles above a certain kin energy cut MCTracksWithHitsOrEnergyCut = False # or of above, factor 2 file size increase compared to MCTracksWithEnergyCutOnly -charmonly = False # option to be set with -A to enable only charm decays, charm x-sec measurement +charmonly = False # option to be set with -A to enable only charm decays, charm x-sec measurement HNL = True inputFile = "/eos/experiment/ship/data/Charm/Cascade-parp16-MSTP82-1-MSEL4-978Bpot.root" @@ -71,6 +73,7 @@ default = '2023' inactivateMuonProcesses = False # provisionally for making studies of various muon background sources +inactivateMuonNuclearProcess = False #Only True for muonDIS parser = ArgumentParser() group = parser.add_mutually_exclusive_group() @@ -170,7 +173,7 @@ if options.A =='b': inputFile = "/eos/experiment/ship/data/Beauty/Cascade-run0-19-parp16-MSTP82-1-MSEL5-5338Bpot.root" if options.A.lower() == 'charmonly': charmonly = True - HNL = False + HNL = False if options.A not in ['b','c','bc','meson','pbrem','qcd']: inclusive = True if options.MM: motherMode=options.MM @@ -197,11 +200,11 @@ #sanity check -if (HNL and options.RPVSUSY) or (HNL and options.DarkPhoton) or (options.DarkPhoton and options.RPVSUSY): +if (HNL and options.RPVSUSY) or (HNL and options.DarkPhoton) or (options.DarkPhoton and options.RPVSUSY): print("cannot have HNL and SUSY or DP at the same time, abort") sys.exit(2) -if (simEngine == "Genie" or simEngine == "nuRadiography") and defaultInputFile: +if (simEngine == "Genie" or simEngine == "nuRadiography") and defaultInputFile: inputFile = "/eos/experiment/ship/data/GenieEvents/genie-nu_mu.root" # "/eos/experiment/ship/data/GenieEvents/genie-nu_mu_bar.root" if simEngine == "muonDIS" and defaultInputFile: @@ -244,11 +247,12 @@ if charmonly: tag = simEngine+"CharmOnly-"+mcEngine if options.eventDisplay: tag = tag+'_D' if options.dv > 4 : tag = 'conical.'+tag +elif dy: tag = str(options.dy)+'.'+tag if not os.path.exists(options.outputDir): os.makedirs(options.outputDir) outFile = "%s/ship.%s.root" % (options.outputDir, tag) -# rm older files !!! +# rm older files !!! for x in os.listdir(options.outputDir): if not x.find(tag)<0: os.system("rm %s/%s" % (options.outputDir, x) ) # Parameter file name @@ -265,8 +269,8 @@ run = ROOT.FairRunSim() run.SetName(mcEngine) # Transport engine run.SetSink(ROOT.FairRootFileSink(outFile)) # Output file -run.SetUserConfig("g4Config.C") # user configuration file default g4Config.C -rtdb = run.GetRuntimeDb() +run.SetUserConfig("g4Config.C") # user configuration file default g4Config.C +rtdb = run.GetRuntimeDb() # -----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 @@ -276,7 +280,7 @@ # -----Create PrimaryGenerator-------------------------------------- primGen = ROOT.FairPrimaryGenerator() if simEngine == "Pythia8": - primGen.SetTarget(ship_geo.target.z0, 0.) + primGen.SetTarget(ship_geo.target.z0, 0.) # -----Pythia8-------------------------------------- if HNL or options.RPVSUSY: P8gen = ROOT.HNLPythia8Generator() @@ -299,7 +303,7 @@ pythia8_conf.configurerpvsusy(P8gen,options.theMass,[theCouplings[0],theCouplings[1]], theCouplings[2],options.RPVSUSYbench,inclusive,options.deepCopy) P8gen.SetParameters("ProcessLevel:all = off") - if inputFile: + if inputFile: ut.checkFileExists(inputFile) # read from external file P8gen.UseExternalFile(inputFile, options.firstEvent) @@ -312,7 +316,7 @@ import pythia8darkphoton_conf passDPconf = pythia8darkphoton_conf.configure(P8gen,options.theMass,options.theDPepsilon,inclusive, motherMode, options.deepCopy) if (passDPconf!=1): sys.exit() - if HNL or options.RPVSUSY or options.DarkPhoton: + if HNL or options.RPVSUSY or options.DarkPhoton: P8gen.SetSmearBeam(1*u.cm) # finite beam size P8gen.SetLmin((ship_geo.Chamber1.z - ship_geo.chambers.Tub1length) - ship_geo.target.z0 ) P8gen.SetLmax(ship_geo.TrackStation1.z - ship_geo.target.z0 ) @@ -320,7 +324,7 @@ primGen.SetTarget(0., 0.) #vertex is setted in pythia8Generator ut.checkFileExists(inputFile) if ship_geo.Box.gausbeam: - primGen.SetBeam(0.,0., 0.5, 0.5) #more central beam, for hits in downstream detectors + primGen.SetBeam(0.,0., 0.5, 0.5) #more central beam, for hits in downstream detectors primGen.SmearGausVertexXY(True) #sigma = x else: primGen.SetBeam(0.,0., ship_geo.Box.TX-1., ship_geo.Box.TY-1.) #Uniform distribution in x/y on the target (0.5 cm of margin at both sides) @@ -332,7 +336,6 @@ # P8gen.SetMom(500.*u.GeV) # P8gen.SetId(-211) primGen.AddGenerator(P8gen) - if simEngine == "FixedTarget": P8gen = ROOT.FixedTargetGenerator() P8gen.SetTarget("volTarget_1",0.,0.) @@ -343,7 +346,7 @@ primGen.AddGenerator(P8gen) if simEngine == "Pythia6": # set muon interaction close to decay volume - primGen.SetTarget(ship_geo.target.z0+ship_geo.muShield.length, 0.) + primGen.SetTarget(ship_geo.target.z0+ship_geo.muShield.length, 0.) # -----Pythia6------------------------- test = ROOT.TPythia6() # don't know any other way of forcing to load lib P6gen = ROOT.tPythia6Generator() @@ -366,7 +369,7 @@ ) # -----Particle Gun----------------------- -if simEngine == "PG": +if simEngine == "PG": myPgun = ROOT.FairBoxGenerator(options.pID,1) myPgun.SetPRange(options.Estart,options.Eend) myPgun.SetPhiRange(0, 360) # // Azimuth angle range [degree] @@ -376,8 +379,9 @@ # -----muon DIS Background------------------------ if simEngine == "muonDIS": ut.checkFileExists(inputFile) - primGen.SetTarget(0., 0.) + primGen.SetTarget(0., 0.) DISgen = ROOT.MuDISGenerator() + # from nu_tau detector to tracking station 2 # mu_start, mu_end = ship_geo.tauMudet.zMudetC,ship_geo.TrackStation2.z # @@ -388,7 +392,8 @@ DISgen.Init(inputFile,options.firstEvent) primGen.AddGenerator(DISgen) options.nEvents = min(options.nEvents,DISgen.GetNevents()) - inactivateMuonProcesses = True # avoid unwanted hadronic events of "incoming" muon flying backward + inactivateMuonProcesses = False + inactivateMuonNuclearProcess = True #avoid double counting print('Generate ',options.nEvents,' with DIS input', ' first event',options.firstEvent) # -----neutrino interactions from nuage------------------------ if simEngine == "Nuage": @@ -407,7 +412,7 @@ nZcells = ntt -1 startx = -ship_geo.NuTauTarget.xdim/2. + nXcells*ship_geo.NuTauTarget.BrX endx = -ship_geo.NuTauTarget.xdim/2. + (nXcells+1)*ship_geo.NuTauTarget.BrX - starty = -ship_geo.NuTauTarget.ydim/2. + nYcells*ship_geo.NuTauTarget.BrY + starty = -ship_geo.NuTauTarget.ydim/2. + nYcells*ship_geo.NuTauTarget.BrY endy = - ship_geo.NuTauTarget.ydim/2. + (nYcells+1)*ship_geo.NuTauTarget.BrY startz = ship_geo.EmuMagnet.zC - ship_geo.NuTauTarget.zdim/2. + ntt *ship_geo.NuTauTT.TTZ + nZcells * ship_geo.NuTauTarget.CellW endz = ship_geo.EmuMagnet.zC - ship_geo.NuTauTarget.zdim/2. + ntt *ship_geo.NuTauTT.TTZ + nZcells * ship_geo.NuTauTarget.CellW + ship_geo.NuTauTarget.BrZ @@ -425,7 +430,7 @@ ut.checkFileExists(inputFile) primGen.SetTarget(0., 0.) # do not interfere with GenieGenerator Geniegen = ROOT.GenieGenerator() - Geniegen.Init(inputFile,options.firstEvent) + Geniegen.Init(inputFile,options.firstEvent) 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()) @@ -435,7 +440,7 @@ ut.checkFileExists(inputFile) primGen.SetTarget(0., 0.) # do not interfere with GenieGenerator Geniegen = ROOT.GenieGenerator() - Geniegen.Init(inputFile,options.firstEvent) + Geniegen.Init(inputFile,options.firstEvent) # Geniegen.SetPositions(ship_geo.target.z0, ship_geo.target.z0, ship_geo.MuonStation3.z) Geniegen.SetPositions(ship_geo.target.z0, ship_geo.tauMudet.zMudetC, ship_geo.MuonStation3.z) Geniegen.NuOnly() @@ -446,7 +451,7 @@ pdg.AddParticle('W','Ion', 1.71350e+02, True, 0., 74, 'XXX', 1000741840) # run.SetPythiaDecayer('DecayConfigPy8.C') - # this requires writing a C macro, would have been easier to do directly in python! + # this requires writing a C macro, would have been easier to do directly in python! # for i in [431,421,411,-431,-421,-411]: # ROOT.gMC.SetUserDecay(i) # Force the decay to be done w/external decayer if simEngine == "Ntuple": @@ -460,10 +465,10 @@ print('Process ',options.nEvents,' from input file') # if simEngine == "MuonBack": -# reading muon tracks from previous Pythia8/Geant4 simulation with charm replaced by cascade production +# reading muon tracks from previous Pythia8/Geant4 simulation with charm replaced by cascade production fileType = ut.checkFileExists(inputFile) if fileType == 'tree': - # 2018 background production + # 2018 background production primGen.SetTarget(ship_geo.target.z0+70.845*u.m,0.) else: primGen.SetTarget(ship_geo.target.z0+50*u.m,0.) @@ -485,20 +490,20 @@ options.nEvents = min(options.nEvents,MuonBackgen.GetNevents()) MCTracksWithHitsOnly = True # otherwise, output file becomes too big print('Process ',options.nEvents,' from input file, with Phi random=',options.phiRandom, ' with MCTracksWithHitsOnly',MCTracksWithHitsOnly) - if options.followMuon : + if options.followMuon : options.fastMuon = True modules['Veto'].SetFollowMuon() if options.fastMuon : modules['Veto'].SetFastMuon() # optional, boost gamma2muon conversion - # ROOT.kShipMuonsCrossSectionFactor = 100. + # ROOT.kShipMuonsCrossSectionFactor = 100. # if simEngine == "Cosmics": primGen.SetTarget(0., 0.) Cosmicsgen = ROOT.CosmicsGenerator() import CMBG_conf CMBG_conf.configure(Cosmicsgen, ship_geo) - if not Cosmicsgen.Init(Opt_high): + if not Cosmicsgen.Init(Opt_high): print("initialization of cosmic background generator failed ",Opt_high) sys.exit(0) Cosmicsgen.n_EVENTS = options.nEvents @@ -524,17 +529,17 @@ elif MCTracksWithEnergyCutOnly: fStack.SetMinPoints(-1) fStack.SetEnergyCut(100.*u.MeV) -elif MCTracksWithHitsOrEnergyCut: +elif MCTracksWithHitsOrEnergyCut: fStack.SetMinPoints(1) fStack.SetEnergyCut(100.*u.MeV) -elif options.deepCopy: +elif options.deepCopy: fStack.SetMinPoints(0) fStack.SetEnergyCut(0.*u.MeV) if options.eventDisplay: # Set cuts for storing the trajectories, can only be done after initialization of run (?!) trajFilter = ROOT.FairTrajFilter.Instance() - trajFilter.SetStepSizeCut(1*u.mm) + trajFilter.SetStepSizeCut(1*u.mm); trajFilter.SetVertexCut(-20*u.m, -20*u.m,ship_geo.target.z0-1*u.m, 20*u.m, 20*u.m, 200.*u.m) trajFilter.SetMomentumCutP(0.1*u.GeV) trajFilter.SetEnergyCut(0., 400.*u.GeV) @@ -572,6 +577,10 @@ gProcessTable = ROOT.G4ProcessTable.GetProcessTable() procmu = gProcessTable.FindProcess(ROOT.G4String('muIoni'),ROOT.G4String('mu+')) procmu.SetVerboseLevel(2) +if inactivateMuonNuclearProcess: + ROOT.gROOT.ProcessLine('#include "Geant4/G4ProcessTable.hh"') + mygMC = ROOT.TGeant4.GetMC() + mygMC.ProcessGeantCommand("/process/inactivate muonNuclear") # -----Start run---------------------------------------------------- run.Run(options.nEvents) # -----Runtime database--------------------------------------------- @@ -595,22 +604,22 @@ fGeo.CheckOverlaps(0.1) # 1 micron takes 5minutes fGeo.PrintOverlaps() # check subsystems in more detail - for x in fGeo.GetTopNode().GetNodes(): + for x in fGeo.GetTopNode().GetNodes(): x.CheckOverlaps(0.0001) fGeo.PrintOverlaps() # -----Finish------------------------------------------------------- timer.Stop() rtime = timer.RealTime() ctime = timer.CpuTime() -print(' ') -print("Macro finished succesfully.") -if "P8gen" in globals() : +print(' ') +print("Macro finished succesfully.") +if "P8gen" in globals() : if (HNL): print("number of retries, events without HNL ",P8gen.nrOfRetries()) - elif (options.DarkPhoton): + elif (options.DarkPhoton): print("number of retries, events without Dark Photons ",P8gen.nrOfRetries()) print("total number of dark photons (including multiple meson decays per single collision) ",P8gen.nrOfDP()) -print("Output file is ", outFile) +print("Output file is ", outFile) print("Parameter file is ",parFile) print("Real time ",rtime, " s, CPU time ",ctime,"s") @@ -632,11 +641,11 @@ nEvents = 0 pointContainers = [] for x in sTree.GetListOfBranches(): - name = x.GetName() + name = x.GetName() if not name.find('Point')<0: pointContainers.append('sTree.'+name+'.GetEntries()') # makes use of convention that all sensitive detectors fill XXXPoint containers for n in range(t.GetEntries()): rc = t.GetEvent(n) - empty = True + empty = True for x in pointContainers: if eval(x)>0: empty = False if not empty: diff --git a/shipgen/MuDISGenerator.cxx b/shipgen/MuDISGenerator.cxx index 2d955faeef..e473859c89 100644 --- a/shipgen/MuDISGenerator.cxx +++ b/shipgen/MuDISGenerator.cxx @@ -15,11 +15,12 @@ using std::cout; using std::endl; -// read events from ntuples produced with MuDIS -// http://MuDIS.hepforge.org/manuals/MuDIS_PhysicsAndUserManual_20130615.pdf + +const Double_t c_light = 29.9792458;//speed of light in cm/ns + // MuDIS momentum GeV // Vertex in SI units, assume this means m -// important to read back number of events to give to FairRoot + // ----- Default constructor ------------------------------------------- MuDISGenerator::MuDISGenerator() {} @@ -28,31 +29,55 @@ MuDISGenerator::MuDISGenerator() {} Bool_t MuDISGenerator::Init(const char* fileName) { return Init(fileName, 0); } -// ----- Default constructor ------------------------------------------- Bool_t MuDISGenerator::Init(const char* fileName, const int firstEvent) { - LOGF(info, "Opening input file %s", fileName); - - iMuon = 0; - dPart = 0; - if (0 == strncmp("/eos",fileName,4) ) { - TString tmp = gSystem->Getenv("EOSSHIP"); - tmp+=fileName; - fInputFile = TFile::Open(tmp); - LOGF(info, "Open external file on eos: %s", tmp.Data()); - }else{ - fInputFile = new TFile(fileName); - } - if (fInputFile->IsZombie() or !fInputFile) { - LOG(FATAL) << "Error opening input file"; - return kFALSE; } - fTree = (TTree *)fInputFile->Get("DIS"); - fNevents = fTree->GetEntries(); - fn = firstEvent; - fTree->SetBranchAddress("InMuon",&iMuon); // incoming muon - fTree->SetBranchAddress("Particles",&dPart); - // cout << "muon DIS Generator number of events "<< fNevents << endl; - return kTRUE; + LOGF(info, "Opening input file %s", fileName); + + iMuon = 0; + dPart = 0; + dPartSoft = 0; + + if (0 == strncmp("/eos",fileName,4)) { + TString tmp = gSystem->Getenv("EOSSHIP"); + tmp += fileName; + fInputFile = TFile::Open(tmp); + LOGF(info, "Open external file on eos: %s", tmp.Data()); + } else { + fInputFile = new TFile(fileName); + } + + if (!fInputFile || fInputFile->IsZombie()) { + LOG(FATAL) << "Error opening input file"; + return kFALSE; + } + + fTree = (TTree *)fInputFile->Get("DIS"); + if (!fTree) { + LOG(FATAL) << "Error: Tree 'DIS' not found in file"; + return kFALSE; + } + + // Print all branches of the TTree + //fTree->Print(); + + // Print out available branches + //TObjArray* branches = fTree->GetListOfBranches(); + //for (Int_t i = 0; i < branches->GetEntries(); i++) { + //TBranch* branch = (TBranch*)branches->At(i); + //cout << "Branch name: " << branch->GetName() << ", Class name: " << branch->GetClassName() << endl; + //} + + fNevents = fTree->GetEntries(); + fn = firstEvent; + + // Set branch addresses + fTree->SetBranchAddress("InMuon", &iMuon); // incoming muon + fTree->SetBranchAddress("DISParticles", &dPart); + fTree->SetBranchAddress("SoftParticles", &dPartSoft); // Soft interaction particles + + cout << "MuDISGenerator: Initialization successful." << endl; + return kTRUE; } + Double_t MuDISGenerator::MeanMaterialBudget(const Double_t *start, const Double_t *end, Double_t *mparam) { // @@ -76,6 +101,7 @@ Double_t MuDISGenerator::MeanMaterialBudget(const Double_t *start, const Double_ // Corrections and improvements by // Andrea Dainese, Andrea.Dainese@lnl.infn.it, // Andrei Gheata, Andrei.Gheata@cern.ch + // Anupama Reghunath, anupama.reghunath@cern.ch // mparam[0]=0; mparam[1]=1; mparam[2] =0; mparam[3] =0; @@ -216,25 +242,31 @@ Bool_t MuDISGenerator::ReadEvent(FairPrimaryGenerator* cpg) if (fn==fNevents) {LOG(WARNING) << "End of input file. Rewind.";} fTree->GetEntry(fn%fNevents); fn++; - if (fn%100==0) { + if (fn%10==0) { cout << "Info MuDISGenerator: MuDIS event-nr "<< fn << endl; } + + int nf = dPart->GetEntries(); + int ns = dPartSoft->GetEntries(); //cout << "*********************************************************" << endl; - //cout << "muon DIS Generator debug " << iMuon->GetEntries()<< " "<< iMuon->AddrAt(0) << " nf "<< nf << " fn=" << fn <GetEntries()<< " "<< iMuon->AddrAt(0) << " nf "<< nf << " fn=" << fn <(iMuon->AddrAt(0)); - //cout << "muon DIS Generator in muon " << int(mu[0][0])<< endl; + //cout << "muon DIS Generator in muon " << int(mu[0][0])<< endl; Double_t x = mu[0][5]*100.; // come in m -> cm Double_t y = mu[0][6]*100.; // come in m -> cm Double_t z = mu[0][7]*100.; // come in m -> cm Double_t w = mu[0][8]; + Double_t t_muon = mu[0][11]; // in ns + // calculate start/end positions along this muon, and amout of material in between + Double_t txmu=mu[0][1]/mu[0][3]; Double_t tymu=mu[0][2]/mu[0][3]; start[0]=x-(z-start[2])*txmu; @@ -257,6 +289,8 @@ Bool_t MuDISGenerator::ReadEvent(FairPrimaryGenerator* cpg) Int_t count=0; //cout << "Info MuDISGenerator Start prob2int while loop, bparam= " << bparam << ", " << bparam*1.e8 <Uniform(0.,1.)) { zmu=gRandom->Uniform(start[2],end[2]); xmu=x-(z-zmu)*txmu; @@ -274,24 +308,57 @@ Bool_t MuDISGenerator::ReadEvent(FairPrimaryGenerator* cpg) //cout << "Info MuDISGenerator: prob2int " << prob2int << ", " << count << endl; //cout << "MuDIS: put position " << xmu << ", " << ymu << ", " << zmu << endl; - //modify weight, by multiplying with average densiy along track?? - cpg->AddTrack(int(mu[0][0]),mu[0][1],mu[0][2],mu[0][3],xmu,ymu,zmu,-1,false,mu[0][4],0.,w); - // in order to have a simulation of possible veto detector hits, let Geant4 follow the muon backward - // due to dE/dx, as soon as the muon travers thick material, this approximation will become bad. - // a proper implementation however would be to have this better integrated in Geant4, follow the muon, call DIS event, continue - cpg->AddTrack(int(mu[0][0]),-mu[0][1],-mu[0][2],-mu[0][3],xmu,ymu,zmu,0,true,mu[0][4],0.,w); -// outgoing particles, [did,dpx,dpy,dpz,E], put density along trajectory as weight, g/cm^2 - w=mparam[0]*mparam[4]; + + double total_mom = TMath::Sqrt(TMath::Power(mu[0][1], 2) + + TMath::Power(mu[0][2], 2) + + TMath::Power(mu[0][3], 2)); //in GeV + + double distance = TMath::Sqrt(TMath::Power(x - xmu, 2) + + TMath::Power(y - ymu, 2) + + TMath::Power(z - zmu, 2)); // in cm + + + double mass = 0.10565999895334244; //muon mass in GeV + double v = c_light*total_mom/TMath::Sqrt(TMath::Power(total_mom,2)+ + TMath::Power(mass,2)); + double t_rmu=distance/v; + + // Adjust time based on the relative positions + if (zmu < z) { + t_rmu = -t_rmu; + } + + double t_DIS=(t_muon+t_rmu)/1e9; //time taken in seconds to reach [xmu,ymu,zmu] + + cpg->AddTrack(int(mu[0][0]),mu[0][1],mu[0][2],mu[0][3],xmu,ymu,zmu,-1,false,mu[0][4],t_DIS,w);//set time of the track start as t_muon from the input file + + // outgoing particles, [did,dpx,dpy,dpz,E], put density along trajectory as weight, g/cm^2 + w=mparam[0]*mparam[4];//modify weight, by multiplying with average densiy along track?? + for(int i=0; i(dPart->AddrAt(i)); - //cout << "muon DIS Generator out part " << int(Part[0][0]) << endl; - //cout << "muon DIS Generator out part mom " << Part[0][1]<<" " << Part[0][2] <<" " << Part[0][3] << " " << Part[0][4] << endl; - //cout << "muon DIS Generator out part pos " << mu[0][5]<<" " << mu[0][6] <<" " << mu[0][7] << endl; - //cout << "muon DIS Generator out part w " << mu[0][8] << endl; - cpg->AddTrack(int(Part[0][0]),Part[0][1],Part[0][2],Part[0][3],xmu,ymu,zmu,0,true,Part[0][4],0.,w); + //cout << "muon DIS Generator out part " << int(Part[0][0]) << endl; + //cout << "muon DIS Generator out part mom " << Part[0][1]<<" " << Part[0][2] <<" " << Part[0][3] << " " << Part[0][4] << endl; + //cout << "muon DIS Generator out part pos " << mu[0][5]<<" " << mu[0][6] <<" " << mu[0][7] << endl; + //cout << "muon DIS Generator out part w " << mu[0][8] << endl; + cpg->AddTrack(int(Part[0][0]),Part[0][1],Part[0][2],Part[0][3],xmu,ymu,zmu,0,true,Part[0][4],t_DIS,w); //it takes times in seconds. //cout << "muon DIS part added "<(dPartSoft->AddrAt(i)); + if ( SoftPart[0][7] > zmu ) { //if the softinteractions are after the DIS point, don't save the soft interaction. + continue; + } + double t_soft=SoftPart[0][8]/1e9; + cpg->AddTrack(int(SoftPart[0][0]), SoftPart[0][1], SoftPart[0][2], SoftPart[0][3], SoftPart[0][5], SoftPart[0][6], SoftPart[0][7], 0, true, SoftPart[0][4], t_soft, w); + + } + return kTRUE; } // ------------------------------------------------------------------------- @@ -299,3 +366,4 @@ Int_t MuDISGenerator::GetNevents() { return fNevents; } + diff --git a/shipgen/MuDISGenerator.h b/shipgen/MuDISGenerator.h index a723e5da58..ec8ae66136 100644 --- a/shipgen/MuDISGenerator.h +++ b/shipgen/MuDISGenerator.h @@ -42,6 +42,7 @@ class MuDISGenerator : public FairGenerator Double_t startZ,endZ; TClonesArray* iMuon ; TClonesArray* dPart ; + TClonesArray* dPartSoft ; FairLogger* fLogger; //! don't make it persistent, magic ROOT command TFile* fInputFile; TTree* fTree; From dceb8858bebe6262caa2671736a16fe292750bac Mon Sep 17 00:00:00 2001 From: anupama-reghunath Date: Thu, 29 Aug 2024 15:30:53 +0200 Subject: [PATCH 02/10] edit run_simScript to inactivate Muon Nuclear for MuDIS --- macro/run_simScript.py | 90 ++++++++++++++++++++---------------------- 1 file changed, 42 insertions(+), 48 deletions(-) diff --git a/macro/run_simScript.py b/macro/run_simScript.py index 6e5eb4bed8..308aa2485f 100755 --- a/macro/run_simScript.py +++ b/macro/run_simScript.py @@ -1,6 +1,4 @@ #!/usr/bin/env python -from __future__ import print_function -from __future__ import division import os import sys import ROOT @@ -33,7 +31,7 @@ MCTracksWithEnergyCutOnly = True # copy particles above a certain kin energy cut MCTracksWithHitsOrEnergyCut = False # or of above, factor 2 file size increase compared to MCTracksWithEnergyCutOnly -charmonly = False # option to be set with -A to enable only charm decays, charm x-sec measurement +charmonly = False # option to be set with -A to enable only charm decays, charm x-sec measurement HNL = True inputFile = "/eos/experiment/ship/data/Charm/Cascade-parp16-MSTP82-1-MSEL4-978Bpot.root" @@ -72,8 +70,7 @@ } default = '2023' -inactivateMuonProcesses = False # provisionally for making studies of various muon background sources -inactivateMuonNuclearProcess = False #Only True for muonDIS +inactivateMuonProcesses = False # provisionally for making studies of various muon background sources parser = ArgumentParser() group = parser.add_mutually_exclusive_group() @@ -153,6 +150,7 @@ const="vacuums", default="helium" ) + parser.add_argument("--SND", dest="SND", help="Activate SND.", action='store_true') parser.add_argument("--noSND", dest="SND", help="Deactivate SND. NOOP, as it currently defaults to off.", action='store_false') @@ -173,7 +171,7 @@ if options.A =='b': inputFile = "/eos/experiment/ship/data/Beauty/Cascade-run0-19-parp16-MSTP82-1-MSEL5-5338Bpot.root" if options.A.lower() == 'charmonly': charmonly = True - HNL = False + HNL = False if options.A not in ['b','c','bc','meson','pbrem','qcd']: inclusive = True if options.MM: motherMode=options.MM @@ -200,11 +198,11 @@ #sanity check -if (HNL and options.RPVSUSY) or (HNL and options.DarkPhoton) or (options.DarkPhoton and options.RPVSUSY): +if (HNL and options.RPVSUSY) or (HNL and options.DarkPhoton) or (options.DarkPhoton and options.RPVSUSY): print("cannot have HNL and SUSY or DP at the same time, abort") sys.exit(2) -if (simEngine == "Genie" or simEngine == "nuRadiography") and defaultInputFile: +if (simEngine == "Genie" or simEngine == "nuRadiography") and defaultInputFile: inputFile = "/eos/experiment/ship/data/GenieEvents/genie-nu_mu.root" # "/eos/experiment/ship/data/GenieEvents/genie-nu_mu_bar.root" if simEngine == "muonDIS" and defaultInputFile: @@ -247,12 +245,11 @@ if charmonly: tag = simEngine+"CharmOnly-"+mcEngine if options.eventDisplay: tag = tag+'_D' if options.dv > 4 : tag = 'conical.'+tag -elif dy: tag = str(options.dy)+'.'+tag if not os.path.exists(options.outputDir): os.makedirs(options.outputDir) outFile = "%s/ship.%s.root" % (options.outputDir, tag) -# rm older files !!! +# rm older files !!! for x in os.listdir(options.outputDir): if not x.find(tag)<0: os.system("rm %s/%s" % (options.outputDir, x) ) # Parameter file name @@ -269,18 +266,17 @@ run = ROOT.FairRunSim() run.SetName(mcEngine) # Transport engine run.SetSink(ROOT.FairRootFileSink(outFile)) # Output file -run.SetUserConfig("g4Config.C") # user configuration file default g4Config.C -rtdb = run.GetRuntimeDb() +run.SetUserConfig("g4Config.C") # user configuration file default g4Config.C +rtdb = run.GetRuntimeDb() # -----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) # -----Create PrimaryGenerator-------------------------------------- primGen = ROOT.FairPrimaryGenerator() if simEngine == "Pythia8": - primGen.SetTarget(ship_geo.target.z0, 0.) + primGen.SetTarget(ship_geo.target.z0, 0.) # -----Pythia8-------------------------------------- if HNL or options.RPVSUSY: P8gen = ROOT.HNLPythia8Generator() @@ -303,7 +299,7 @@ pythia8_conf.configurerpvsusy(P8gen,options.theMass,[theCouplings[0],theCouplings[1]], theCouplings[2],options.RPVSUSYbench,inclusive,options.deepCopy) P8gen.SetParameters("ProcessLevel:all = off") - if inputFile: + if inputFile: ut.checkFileExists(inputFile) # read from external file P8gen.UseExternalFile(inputFile, options.firstEvent) @@ -316,7 +312,7 @@ import pythia8darkphoton_conf passDPconf = pythia8darkphoton_conf.configure(P8gen,options.theMass,options.theDPepsilon,inclusive, motherMode, options.deepCopy) if (passDPconf!=1): sys.exit() - if HNL or options.RPVSUSY or options.DarkPhoton: + if HNL or options.RPVSUSY or options.DarkPhoton: P8gen.SetSmearBeam(1*u.cm) # finite beam size P8gen.SetLmin((ship_geo.Chamber1.z - ship_geo.chambers.Tub1length) - ship_geo.target.z0 ) P8gen.SetLmax(ship_geo.TrackStation1.z - ship_geo.target.z0 ) @@ -324,7 +320,7 @@ primGen.SetTarget(0., 0.) #vertex is setted in pythia8Generator ut.checkFileExists(inputFile) if ship_geo.Box.gausbeam: - primGen.SetBeam(0.,0., 0.5, 0.5) #more central beam, for hits in downstream detectors + primGen.SetBeam(0.,0., 0.5, 0.5) #more central beam, for hits in downstream detectors primGen.SmearGausVertexXY(True) #sigma = x else: primGen.SetBeam(0.,0., ship_geo.Box.TX-1., ship_geo.Box.TY-1.) #Uniform distribution in x/y on the target (0.5 cm of margin at both sides) @@ -346,7 +342,7 @@ primGen.AddGenerator(P8gen) if simEngine == "Pythia6": # set muon interaction close to decay volume - primGen.SetTarget(ship_geo.target.z0+ship_geo.muShield.length, 0.) + primGen.SetTarget(ship_geo.target.z0+ship_geo.muShield.length, 0.) # -----Pythia6------------------------- test = ROOT.TPythia6() # don't know any other way of forcing to load lib P6gen = ROOT.tPythia6Generator() @@ -369,7 +365,7 @@ ) # -----Particle Gun----------------------- -if simEngine == "PG": +if simEngine == "PG": myPgun = ROOT.FairBoxGenerator(options.pID,1) myPgun.SetPRange(options.Estart,options.Eend) myPgun.SetPhiRange(0, 360) # // Azimuth angle range [degree] @@ -379,9 +375,8 @@ # -----muon DIS Background------------------------ if simEngine == "muonDIS": ut.checkFileExists(inputFile) - primGen.SetTarget(0., 0.) + primGen.SetTarget(0., 0.) DISgen = ROOT.MuDISGenerator() - # from nu_tau detector to tracking station 2 # mu_start, mu_end = ship_geo.tauMudet.zMudetC,ship_geo.TrackStation2.z # @@ -392,8 +387,6 @@ DISgen.Init(inputFile,options.firstEvent) primGen.AddGenerator(DISgen) options.nEvents = min(options.nEvents,DISgen.GetNevents()) - inactivateMuonProcesses = False - inactivateMuonNuclearProcess = True #avoid double counting print('Generate ',options.nEvents,' with DIS input', ' first event',options.firstEvent) # -----neutrino interactions from nuage------------------------ if simEngine == "Nuage": @@ -412,7 +405,7 @@ nZcells = ntt -1 startx = -ship_geo.NuTauTarget.xdim/2. + nXcells*ship_geo.NuTauTarget.BrX endx = -ship_geo.NuTauTarget.xdim/2. + (nXcells+1)*ship_geo.NuTauTarget.BrX - starty = -ship_geo.NuTauTarget.ydim/2. + nYcells*ship_geo.NuTauTarget.BrY + starty = -ship_geo.NuTauTarget.ydim/2. + nYcells*ship_geo.NuTauTarget.BrY endy = - ship_geo.NuTauTarget.ydim/2. + (nYcells+1)*ship_geo.NuTauTarget.BrY startz = ship_geo.EmuMagnet.zC - ship_geo.NuTauTarget.zdim/2. + ntt *ship_geo.NuTauTT.TTZ + nZcells * ship_geo.NuTauTarget.CellW endz = ship_geo.EmuMagnet.zC - ship_geo.NuTauTarget.zdim/2. + ntt *ship_geo.NuTauTT.TTZ + nZcells * ship_geo.NuTauTarget.CellW + ship_geo.NuTauTarget.BrZ @@ -430,7 +423,7 @@ ut.checkFileExists(inputFile) primGen.SetTarget(0., 0.) # do not interfere with GenieGenerator Geniegen = ROOT.GenieGenerator() - Geniegen.Init(inputFile,options.firstEvent) + Geniegen.Init(inputFile,options.firstEvent) 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()) @@ -440,7 +433,7 @@ ut.checkFileExists(inputFile) primGen.SetTarget(0., 0.) # do not interfere with GenieGenerator Geniegen = ROOT.GenieGenerator() - Geniegen.Init(inputFile,options.firstEvent) + Geniegen.Init(inputFile,options.firstEvent) # Geniegen.SetPositions(ship_geo.target.z0, ship_geo.target.z0, ship_geo.MuonStation3.z) Geniegen.SetPositions(ship_geo.target.z0, ship_geo.tauMudet.zMudetC, ship_geo.MuonStation3.z) Geniegen.NuOnly() @@ -451,7 +444,7 @@ pdg.AddParticle('W','Ion', 1.71350e+02, True, 0., 74, 'XXX', 1000741840) # run.SetPythiaDecayer('DecayConfigPy8.C') - # this requires writing a C macro, would have been easier to do directly in python! + # this requires writing a C macro, would have been easier to do directly in python! # for i in [431,421,411,-431,-421,-411]: # ROOT.gMC.SetUserDecay(i) # Force the decay to be done w/external decayer if simEngine == "Ntuple": @@ -465,10 +458,10 @@ print('Process ',options.nEvents,' from input file') # if simEngine == "MuonBack": -# reading muon tracks from previous Pythia8/Geant4 simulation with charm replaced by cascade production +# reading muon tracks from previous Pythia8/Geant4 simulation with charm replaced by cascade production fileType = ut.checkFileExists(inputFile) if fileType == 'tree': - # 2018 background production + # 2018 background production primGen.SetTarget(ship_geo.target.z0+70.845*u.m,0.) else: primGen.SetTarget(ship_geo.target.z0+50*u.m,0.) @@ -490,20 +483,20 @@ options.nEvents = min(options.nEvents,MuonBackgen.GetNevents()) MCTracksWithHitsOnly = True # otherwise, output file becomes too big print('Process ',options.nEvents,' from input file, with Phi random=',options.phiRandom, ' with MCTracksWithHitsOnly',MCTracksWithHitsOnly) - if options.followMuon : + if options.followMuon : options.fastMuon = True modules['Veto'].SetFollowMuon() if options.fastMuon : modules['Veto'].SetFastMuon() # optional, boost gamma2muon conversion - # ROOT.kShipMuonsCrossSectionFactor = 100. + # ROOT.kShipMuonsCrossSectionFactor = 100. # if simEngine == "Cosmics": primGen.SetTarget(0., 0.) Cosmicsgen = ROOT.CosmicsGenerator() import CMBG_conf CMBG_conf.configure(Cosmicsgen, ship_geo) - if not Cosmicsgen.Init(Opt_high): + if not Cosmicsgen.Init(Opt_high): print("initialization of cosmic background generator failed ",Opt_high) sys.exit(0) Cosmicsgen.n_EVENTS = options.nEvents @@ -529,17 +522,17 @@ elif MCTracksWithEnergyCutOnly: fStack.SetMinPoints(-1) fStack.SetEnergyCut(100.*u.MeV) -elif MCTracksWithHitsOrEnergyCut: +elif MCTracksWithHitsOrEnergyCut: fStack.SetMinPoints(1) fStack.SetEnergyCut(100.*u.MeV) -elif options.deepCopy: +elif options.deepCopy: fStack.SetMinPoints(0) fStack.SetEnergyCut(0.*u.MeV) if options.eventDisplay: # Set cuts for storing the trajectories, can only be done after initialization of run (?!) trajFilter = ROOT.FairTrajFilter.Instance() - trajFilter.SetStepSizeCut(1*u.mm); + trajFilter.SetStepSizeCut(1*u.mm) trajFilter.SetVertexCut(-20*u.m, -20*u.m,ship_geo.target.z0-1*u.m, 20*u.m, 20*u.m, 200.*u.m) trajFilter.SetMomentumCutP(0.1*u.GeV) trajFilter.SetEnergyCut(0., 400.*u.GeV) @@ -565,7 +558,7 @@ #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') -if inactivateMuonProcesses : +if inactivateMuonProcesses: ROOT.gROOT.ProcessLine('#include "Geant4/G4ProcessTable.hh"') mygMC = ROOT.TGeant4.GetMC() mygMC.ProcessGeantCommand("/process/inactivate muPairProd") @@ -577,10 +570,11 @@ gProcessTable = ROOT.G4ProcessTable.GetProcessTable() procmu = gProcessTable.FindProcess(ROOT.G4String('muIoni'),ROOT.G4String('mu+')) procmu.SetVerboseLevel(2) -if inactivateMuonNuclearProcess: - ROOT.gROOT.ProcessLine('#include "Geant4/G4ProcessTable.hh"') - mygMC = ROOT.TGeant4.GetMC() - mygMC.ProcessGeantCommand("/process/inactivate muonNuclear") + +if simEngine == "muonDIS": + ROOT.gROOT.ProcessLine('#include "Geant4/G4ProcessTable.hh"') + mygMC = ROOT.TGeant4.GetMC() + mygMC.ProcessGeantCommand("/process/inactivate muonNuclear") # -----Start run---------------------------------------------------- run.Run(options.nEvents) # -----Runtime database--------------------------------------------- @@ -604,22 +598,22 @@ fGeo.CheckOverlaps(0.1) # 1 micron takes 5minutes fGeo.PrintOverlaps() # check subsystems in more detail - for x in fGeo.GetTopNode().GetNodes(): + for x in fGeo.GetTopNode().GetNodes(): x.CheckOverlaps(0.0001) fGeo.PrintOverlaps() # -----Finish------------------------------------------------------- timer.Stop() rtime = timer.RealTime() ctime = timer.CpuTime() -print(' ') -print("Macro finished succesfully.") -if "P8gen" in globals() : +print(' ') +print("Macro finished succesfully.") +if "P8gen" in globals() : if (HNL): print("number of retries, events without HNL ",P8gen.nrOfRetries()) - elif (options.DarkPhoton): + elif (options.DarkPhoton): print("number of retries, events without Dark Photons ",P8gen.nrOfRetries()) print("total number of dark photons (including multiple meson decays per single collision) ",P8gen.nrOfDP()) -print("Output file is ", outFile) +print("Output file is ", outFile) print("Parameter file is ",parFile) print("Real time ",rtime, " s, CPU time ",ctime,"s") @@ -641,11 +635,11 @@ nEvents = 0 pointContainers = [] for x in sTree.GetListOfBranches(): - name = x.GetName() + name = x.GetName() if not name.find('Point')<0: pointContainers.append('sTree.'+name+'.GetEntries()') # makes use of convention that all sensitive detectors fill XXXPoint containers for n in range(t.GetEntries()): rc = t.GetEvent(n) - empty = True + empty = True for x in pointContainers: if eval(x)>0: empty = False if not empty: From 08e00635b756f069ccb2e396cf191b89dc3f082d Mon Sep 17 00:00:00 2001 From: anupama-reghunath Date: Sat, 7 Sep 2024 12:39:48 +0200 Subject: [PATCH 03/10] Added MuonDIS directory for updated scripts --- macro/run_simScript.py | 19 --- muonDIS/add_muonresponse.py | 88 ++++++++++++++ muonDIS/makeMuonDIS.py | 226 +++++++++++++++++++++++++++++++++++ muonDIS/make_nTuple_SBT.py | 202 +++++++++++++++++++++++++++++++ muonDIS/make_nTuple_Tr.py | 231 ++++++++++++++++++++++++++++++++++++ shipgen/MuDISGenerator.cxx | 78 ++++++------ shipgen/MuDISGenerator.h | 4 + 7 files changed, 784 insertions(+), 64 deletions(-) create mode 100644 muonDIS/add_muonresponse.py create mode 100644 muonDIS/makeMuonDIS.py create mode 100644 muonDIS/make_nTuple_SBT.py create mode 100644 muonDIS/make_nTuple_Tr.py diff --git a/macro/run_simScript.py b/macro/run_simScript.py index 308aa2485f..020483d9cb 100755 --- a/macro/run_simScript.py +++ b/macro/run_simScript.py @@ -70,8 +70,6 @@ } default = '2023' -inactivateMuonProcesses = False # provisionally for making studies of various muon background sources - parser = ArgumentParser() group = parser.add_mutually_exclusive_group() parser.add_argument("--evtcalc", help="Use EventCalc", action="store_true") @@ -558,23 +556,6 @@ #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') -if inactivateMuonProcesses: - ROOT.gROOT.ProcessLine('#include "Geant4/G4ProcessTable.hh"') - mygMC = ROOT.TGeant4.GetMC() - mygMC.ProcessGeantCommand("/process/inactivate muPairProd") - mygMC.ProcessGeantCommand("/process/inactivate muBrems") - mygMC.ProcessGeantCommand("/process/inactivate muIoni") - mygMC.ProcessGeantCommand("/process/inactivate muonNuclear") - mygMC.ProcessGeantCommand("/particle/select mu+") - mygMC.ProcessGeantCommand("/particle/process/dump") - gProcessTable = ROOT.G4ProcessTable.GetProcessTable() - procmu = gProcessTable.FindProcess(ROOT.G4String('muIoni'),ROOT.G4String('mu+')) - procmu.SetVerboseLevel(2) - -if simEngine == "muonDIS": - ROOT.gROOT.ProcessLine('#include "Geant4/G4ProcessTable.hh"') - mygMC = ROOT.TGeant4.GetMC() - mygMC.ProcessGeantCommand("/process/inactivate muonNuclear") # -----Start run---------------------------------------------------- run.Run(options.nEvents) # -----Runtime database--------------------------------------------- diff --git a/muonDIS/add_muonresponse.py b/muonDIS/add_muonresponse.py new file mode 100644 index 0000000000..8e98ff932d --- /dev/null +++ b/muonDIS/add_muonresponse.py @@ -0,0 +1,88 @@ +#!/usr/bin/env python3 + +"""Script to add the incoming muon's MC hits in the SBT to the simulation file.""" + +import argparse +import logging + +import ROOT as r + +parser = argparse.ArgumentParser(description=__doc__) +parser.add_argument( + "-f", + "--inputfile", + required=True, + help="Simulation file produced from DIS, WARNING: File will be modified.", +) +parser.add_argument( + "-m", + "--muonfile", + required=True, + help="Muon file used for DIS production (muonDis_<>.root)", +) +args = parser.parse_args() + +# Open the existing ROOT file in update mode +logging.warning(f"{args.inputfile} opened for modification") +inputFile = r.TFile.Open(args.inputfile, "UPDATE") +try: + input_tree = inputFile.cbmsim +except Exception as e: + print(f"Error: {e}") + inputFile.Close() + exit(1) + +# Open the external file with additional vetoPoints +muonFile = r.TFile.Open(args.muonfile, "READ") +try: + muon_tree = muonFile.DIS +except Exception as e: + print(f"Error: {e}") + muonFile.Close() + exit(1) + +inputFile.cd() +output_tree = input_tree.CloneTree( + 0 +) # Clone the structure of the existing tree, but do not copy the entries + +# Access the vetoPoint branch in the original and external trees +originalVetoPoints = r.TClonesArray("vetoPoint") +muonVetoPoints = r.TClonesArray("vetoPoint") +combinedVetoPoints = r.TClonesArray("vetoPoint") + +input_tree.SetBranchAddress("vetoPoint", originalVetoPoints) +muon_tree.SetBranchAddress("muon_vetoPoints", muonVetoPoints) +output_tree.SetBranchAddress("vetoPoint", combinedVetoPoints) + + +for input_event, muon_event in zip(input_tree, muon_tree): + interaction_point = r.TVector3() + input_event.MCTrack[0].GetStartVertex(interaction_point) + + combinedVetoPoints.Clear() + + index = 0 + + for hit in originalVetoPoints: + if combinedVetoPoints.GetSize() == index: + combinedVetoPoints.Expand(index + 1) + combinedVetoPoints[index] = hit + index += 1 + + for hit in muonVetoPoints: + if hit.GetZ() < interaction_point.Z(): + if combinedVetoPoints.GetSize() == index: + combinedVetoPoints.Expand(index + 1) + combinedVetoPoints[index] = hit + index += 1 + + output_tree.Fill() + +# Write the updated tree back to the same file +inputFile.cd() # Ensure we are in the correct directory of the ROOT file +output_tree.Write("cbmsim", r.TObject.kOverwrite) +inputFile.Close() +muonFile.Close() + +print(f"Updated file saved as {args.inputfile}") diff --git a/muonDIS/makeMuonDIS.py b/muonDIS/makeMuonDIS.py new file mode 100644 index 0000000000..69fa271d23 --- /dev/null +++ b/muonDIS/makeMuonDIS.py @@ -0,0 +1,226 @@ +#!/usr/bin/env python +"""Script to generate DIS events for muons in Pythia6, and sve to a ROOT file (along with the original muon's soft interactions).""" + +import argparse +import logging +import time +from array import array + +import ROOT as r + +logging.basicConfig(level=logging.INFO) +PDG = r.TDatabasePDG.Instance() + + +def rotate(px, py, pz, theta, phi): + """Rotate the daughter particle momentum to align with respect to the muon's momentum.""" + momentum = r.TVector3(px, py, pz) + + rotation = r.TRotation() + rotation.RotateY(theta) # Rotate around the Y-axis + rotation.RotateZ(phi) # Rotate around the Z-axis + + # Apply the rotation to the momentum vector + rotated_momentum = rotation * momentum + + return rotated_momentum.X(), rotated_momentum.Y(), rotated_momentum.Z() + + +def makeMuonDIS(): + """Generate DIS events.""" + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("-f", "--inputFile", help="Input file to use", required=True) + parser.add_argument( + "-nJob", "--nJob", type=int, help="Process ID, gives muon index", required=True + ) + parser.add_argument( + "-nPerJobs", + "--nPerJobs", + type=int, + help="The number of muons per file", + required=True, + ) + parser.add_argument( + "-nDISPerMuon", + "--nDIS", + type=int, + help="Number of DIS per muon to generate", + required=True, + ) + + args = parser.parse_args() + + logging.info(f"Opening input file: {args.inputFile}") + muonFile = r.TFile.Open(args.inputFile, "read") + + try: + muon_tree = muonFile.MuonAndSoftInteractions + except Exception as e: + print(f"Error: {e}") + muonFile.Close() + exit(1) + + nPerJob = args.nPerJobs # number of muons handled by the python script + nStart = args.nPerJobs * args.nJob + logging.debug(f"Total entries in the tree: {muon_tree.GetEntries()}") + nEnd = min(muon_tree.GetEntries(), nStart + nPerJob) + + logging.info(f"Creating output file: muonDis_{args.nJob}.root") + + outputFile = r.TFile.Open(f"muonDis_{args.nJob}.root", "recreate") + output_tree = r.TTree("DIS", "muon DIS") + + iMuon = r.TClonesArray("TVectorD") + output_tree.Branch("InMuon", iMuon, 32000, -1) + + dPartDIS = r.TClonesArray("TVectorD") + output_tree.Branch("DISParticles", dPartDIS, 32000, -1) + + dPartSoft = r.TClonesArray("TVectorD") + output_tree.Branch("SoftParticles", dPartSoft, 32000, -1) + + muon_vetoPoints = r.TClonesArray("vetoPoint") + output_tree.Branch("muon_vetoPoints", muon_vetoPoints, 32000, -1) + + myPythia = r.TPythia6() + myPythia.SetMSEL(2) + myPythia.SetPARP(2, 2) + for kf in [211, 321, 130, 310, 3112, 3122, 3222, 3312, 3322, 3334]: + kc = myPythia.Pycomp(kf) + myPythia.SetMDCY(kc, 1, 0) + + seed = int(time.time() % 900000000) + myPythia.SetMRPY(1, seed) + mutype = {-13: "gamma/mu+", 13: "gamma/mu-"} + + myPythia.SetMSTU(11, 11) + logging.info(f"Processing events from {nStart} to {nEnd}...") + + nMade = 0 + + for k in range(nStart, nEnd): + muon_tree.GetEvent(k) + + imuondata = muon_tree.imuondata + + pid = imuondata[0] + px = imuondata[1] + py = imuondata[2] + pz = imuondata[3] + x = imuondata[4] + y = imuondata[5] + z = imuondata[6] + w = imuondata[7] + time_muon = imuondata[8] + + p = r.TMath.Sqrt(px**2 + py**2 + pz**2) + mass = PDG.GetParticle(abs(int(pid))).Mass() + E = r.TMath.Sqrt(mass**2 + p**2) + + theta = r.TMath.ACos(pz / p) + phi = r.TMath.ATan2(py, px) + + logging.debug( + f"Muon: PID={pid}, px={px}, py={py}, pz={pz}, E={E}, x={x}, y={y}, z={z}, w={w}" + ) + + isProton = 1 + xsec = 0 + + mu = array("d", [pid, px, py, pz, E, x, y, z, w, isProton, xsec, time_muon]) + muPart = r.TVectorD(12, mu) + myPythia.Initialize("FIXT", mutype[pid], "p+", p) # target = "p+" + myPythia.Pylist(1) + + for a in range(args.nDIS): + if a == args.nDIS // 2: + myPythia.Initialize("FIXT", mutype[pid], "n0", p) # target = "n0" + isProton = 0 + logging.info("Switching to neutron interaction") + + dPartDIS.Clear() + iMuon.Clear() + muPart[9] = isProton + iMuon[0] = muPart + myPythia.GenerateEvent() + myPythia.Pyedit(1) + logging.info(f"Event {a} generated, number of particles: {myPythia.GetN()}") + + for itrk in range(1, myPythia.GetN() + 1): + xsec = myPythia.GetPARI(1) + muPart[10] = xsec + did = myPythia.GetK(itrk, 2) + dpx, dpy, dpz = rotate( + myPythia.GetP(itrk, 1), + myPythia.GetP(itrk, 2), + myPythia.GetP(itrk, 3), + theta, + phi, + ) + psq = dpx**2 + dpy**2 + dpz**2 + masssq = PDG.GetParticle(did).Mass() ** 2 + E = r.TMath.Sqrt(masssq + psq) + m = array("d", [did, dpx, dpy, dpz, E]) + part = r.TVectorD(5, m) + nPart = dPartDIS.GetEntries() + if dPartDIS.GetSize() == nPart: + dPartDIS.Expand(nPart + 10) + # dPartDIS.ConstructedAt(nPart).Use(part) #to be adapted later + dPartDIS[nPart] = part + if itrk == 1: + with open(f"sigmadata_{args.nJob}.txt", "a") as fcross: + fcross.write(f"{xsec}\n") + + dPartSoft.Clear() + + for softTrack in muon_tree.tracks: + did = softTrack.GetPdgCode() + dpx = softTrack.GetPx() + dpy = softTrack.GetPy() + dpz = softTrack.GetPz() + + psq = dpx**2 + dpy**2 + dpz**2 + masssq = PDG.GetParticle(did).Mass() ** 2 + E = r.TMath.Sqrt(masssq + psq) + + softx = softTrack.GetStartX() + softy = softTrack.GetStartY() + softz = softTrack.GetStartZ() + time_ = softTrack.GetStartT() + + m = array("d", [did, dpx, dpy, dpz, E, softx, softy, softz, time_]) + + part = r.TVectorD(9, m) + nPart = dPartSoft.GetEntries() + if dPartSoft.GetSize() == nPart: + dPartSoft.Expand(nPart + 10) + # dPartSoft.ConstructedAt(nPart).Use(part) #to be adapted later + dPartSoft[nPart] = part + + muon_vetoPoints.Clear() + + index = 0 + for hit in muon_tree.muon_vetoPoints: + if muon_vetoPoints.GetSize() == index: + muon_vetoPoints.Expand(index + 1) + hit.SetTrackID(0) # Set TrackID to match for muon ID for new simulation + muon_vetoPoints[index] = hit + index += 1 + + output_tree.Fill() + + nMade += 1 + + if nMade % 10 == 0: + logging.info(f"Muons processed: {nMade}") + + outputFile.cd() + output_tree.Write() + myPythia.SetMSTU(11, 6) + logging.info( + f"Created DIS for saved in muonDis_{args.nJob}.root, muon event index from {nStart} - {nEnd} , nDISPerMuon {args.nDIS}" + ) + + +if __name__ == "__main__": + makeMuonDIS() diff --git a/muonDIS/make_nTuple_SBT.py b/muonDIS/make_nTuple_SBT.py new file mode 100644 index 0000000000..86e6f7317f --- /dev/null +++ b/muonDIS/make_nTuple_SBT.py @@ -0,0 +1,202 @@ +#!/usr/bin/env python +"""Script to collect muons hitting the SBT (including soft interaction products) to a ROOT file.""" + +import argparse +import logging +import os + +import ROOT as r +import shipunit as u + +# Argument parser setup +parser = argparse.ArgumentParser(description=__doc__) +parser.add_argument("-test", dest="testing_code", help="Run Test", action="store_true") +parser.add_argument( + "-p", dest="path", help="path to muon background files", required=False +) +parser.add_argument( + "-o", + "--outputfile", + default="muonsProduction_wsoft_SBT.root", + help="custom outputfile name", +) +args = parser.parse_args() + +# Histogram setup +hnumSegPermmuon = r.TH1I( + "hnumSegPermmuon", "Numbers of fired segments per muon", 200, 0.0, 200 +) +hPmuon = r.TH1F("hPmuon", "The momentum of the muons hitting the SBT", 400, 0.0, 400) + +# Initialize variables +ev = 0 +processed_events = set() + +if args.testing_code: + print( + "test code, output file name overwritten as: muonsProduction_wsoft_SBT_test.root" + ) + args.outputfile = "muonsProduction_wsoft_SBT_test.root" + selectedmuons = "SelectedMuonsSBT_test.txt" +else: + selectedmuons = "SelectedMuonsSBT.txt" + +if args.path: + path = args.path +else: + path = "/eos/experiment/ship/simulation/bkg/MuonBack_2024helium/sc_v6_10_spills" + +fsel = open(selectedmuons, "w") + +output_file = r.TFile.Open(args.outputfile, "recreate") +output_tree = r.TTree( + "MuonAndSoftInteractions", "Muon information and soft interaction tracks" +) + +imuondata = r.TVectorD(9) # 9 values: pid, px, py, pz, x, y, z, weight,time_of_hit +output_tree.Branch("imuondata", imuondata) + +track_array = r.TObjArray() +output_tree.Branch("tracks", track_array) + +muon_vetoPoints = r.TClonesArray("vetoPoint") +output_tree.Branch("muon_vetoPoints", muon_vetoPoints) + +logging.basicConfig(level=logging.INFO) + +for inputFolder in os.listdir(path): + for subFolder in os.listdir(os.path.join(path, inputFolder)): + if not os.path.isdir(os.path.join(path, inputFolder, subFolder)): + continue + + if args.testing_code and ev >= 100000: + break + + logging.info(f"Processing folder: {inputFolder}/{subFolder}") + f = r.TFile.Open( + os.path.join( + path, inputFolder, subFolder, "ship.conical.MuonBack-TGeant4.root" + ), + "read", + ) + + try: + tree = f.cbmsim + except Exception as e: + print(f"Error :{e}") + f.Close() + continue + + for event in tree: + ev += 1 + numHitsPermuon = 0 + + # saving soft tracks + track_array.Clear() + + muon_id = None + for itrk in range(event.MCTrack.GetEntries()): + # loops through MCTracks to find the incoming Muon's track id. + if abs(event.MCTrack[itrk].GetPdgCode()) == 13: + muon_id = itrk + break + + for track in event.MCTrack: + if track.GetMotherId() == muon_id and ( + not track.GetProcName().Data() == "Muon nuclear interaction" + ): + track_array.Add(track) + + index = 0 + muon_vetoPoints.Clear() + + # saving the incoming muon's veto response + for hit in event.vetoPoint: + detID = hit.GetDetectorID() + pid = hit.PdgCode() + if 1000 < detID < 999999 and abs(pid) == 13: + if muon_vetoPoints.GetSize() == index: + muon_vetoPoints.Expand(index + 1) + muon_vetoPoints[index] = hit + index += 1 + + for hit in event.vetoPoint: + detID = hit.GetDetectorID() + pid = hit.PdgCode() + trackID = hit.GetTrackID() + + if abs(pid) == 13: + numHitsPermuon += 1 + if ev not in processed_events: + processed_events.add(ev) + P = r.TMath.Sqrt( + hit.GetPx() ** 2 + hit.GetPy() ** 2 + hit.GetPz() ** 2 + ) + weight = event.MCTrack[trackID].GetWeight() + hPmuon.Fill(P, weight) + + if P > 3 / u.GeV: + imuondata[0] = float(pid) + imuondata[1] = float(hit.GetPx() / u.GeV) + imuondata[2] = float(hit.GetPy() / u.GeV) + imuondata[3] = float(hit.GetPz() / u.GeV) + imuondata[4] = float(hit.GetX() / u.m) + imuondata[5] = float(hit.GetY() / u.m) + imuondata[6] = float(hit.GetZ() / u.m) + imuondata[7] = float(weight) + imuondata[8] = float(hit.GetTime()) + + output_tree.Fill() + + fsel.write( + f"{pid} {hit.GetPx() / u.GeV} {hit.GetPy() / u.GeV} {hit.GetPz() / u.GeV} {hit.GetX() / u.m} {hit.GetY() / u.m} {hit.GetZ() / u.m} {weight}\n" + ) + + if numHitsPermuon != 0: + hnumSegPermmuon.Fill(numHitsPermuon) + + f.Close() + + +output_file.cd() +output_tree.Write() +hnumSegPermmuon.Write() +hPmuon.Write() +output_file.Close() +fsel.close() + +print( + "------------------------------------------------------file saved, reading", + args.outputfile, + " now----------------------------------------------------------------", +) + +with r.TFile.Open(args.outputfile, "read") as file: + try: + tree = file.MuonAndSoftInteractions + except Exception as e: + print(f"Error: {e}") + exit(1) + + print(f"Processing tree: {tree.GetName()}") + print(f"Total number of entries: {tree.GetEntries()}") + + for event in tree: + imuondata = event.imuondata + pid = imuondata[0] + px = imuondata[1] + py = imuondata[2] + pz = imuondata[3] + x = imuondata[4] + y = imuondata[5] + z = imuondata[6] + weight = imuondata[7] + time_hit = imuondata[8] + num_tracks = len(event.tracks) + num_muonhits = len(event.muon_vetoPoints) + + print( + f"Muon PID: {pid}, x: {x}, y: {y}, z: {z}, t_muon: {time_hit}, " + f"Number of soft tracks in this event: {num_tracks}, " + f"Number of SBT hits from the muon: {num_muonhits}" + ) diff --git a/muonDIS/make_nTuple_Tr.py b/muonDIS/make_nTuple_Tr.py new file mode 100644 index 0000000000..3e04d4f59c --- /dev/null +++ b/muonDIS/make_nTuple_Tr.py @@ -0,0 +1,231 @@ +#!/usr/bin/env python +"""Script to collect muons hitting the 1st Tracking Station (including soft interaction products) to a ROOT file.""" + +import argparse +import logging +import os + +import ROOT as r +import shipunit as u + +# Histogram +hPmuon = r.TH1F( + "hPmuon", "Momentum of muons hitting the SBT;P [GeV];Entries", 400, 0.0, 400 +) + +# Argument parser +parser = argparse.ArgumentParser(description=__doc__) +parser.add_argument( + "-o", + "--outputfile", + default="muonsProduction_wsoft_Tr.root", + help="custom outputfile name", +) +parser.add_argument( + "-p", dest="path", help="path to muon background files", required=False +) +parser.add_argument("-test", dest="testing_code", help="Run Test", action="store_true") +args = parser.parse_args() + +# Variables +ev = 0 +evs = [0, 0, 0] +processed_events = set() + + +if args.testing_code: + print( + "test code, output file name overwritten as: muonsProduction_wsoft_Tr_test.root" + ) + args.outputfile = "muonsProduction_wsoft_Tr_test.root" + selectedmuons = "SelectedMuonsTr_test.txt" +else: + selectedmuons = "SelectedMuonsTr.txt" + +# Output file and tree setup +output_file = r.TFile.Open(args.outputfile, "recreate") +output_tree = r.TTree( + "MuonAndSoftInteractions", "Muon information and soft interaction tracks" +) + +imuondata = r.TVectorD(9) # 9 values: pid, px, py, pz, x, y, z, weight,time_of_hit +output_tree.Branch("imuondata", imuondata) + +track_array = r.TObjArray() +output_tree.Branch("tracks", track_array) + +muon_vetoPoints = r.TClonesArray("vetoPoint") +output_tree.Branch("muon_vetoPoints", muon_vetoPoints) + +if args.path: + path = args.path +else: + path = "/eos/experiment/ship/simulation/bkg/MuonBack_2024helium/sc_v6_10_spills" + +fsel = open(selectedmuons, "w") + +logging.basicConfig(level=logging.INFO) + +for inputFolder in os.listdir(path): + for subFolder in os.listdir(os.path.join(path, inputFolder)): + if not os.path.isdir(os.path.join(path, inputFolder, subFolder)): + continue + + if args.testing_code and evs[2] >= 2: + break + + logging.info(f"Processing folder: {inputFolder}/{subFolder}") + f = r.TFile.Open( + os.path.join( + path, inputFolder, subFolder, "ship.conical.MuonBack-TGeant4.root" + ), + "read", + ) + + try: + tree = f.cbmsim + except Exception as e: + print(f"Error :{e}") + f.Close() + continue + + for event in tree: + ev += 1 + + # saving soft tracks + track_array.Clear() + + muon_id = None + for itrk in range(event.MCTrack.GetEntries()): + # loops through MCTracks to find the incoming Muon's track id. + if abs(event.MCTrack[itrk].GetPdgCode()) == 13: + muon_id = itrk + break + + for track in event.MCTrack: + if track.GetMotherId() == muon_id and ( + not track.GetProcName().Data() == "Muon nuclear interaction" + ): + track_array.Add(track) + + # saving the incoming muon's veto response + index = 0 + muon_vetoPoints.Clear() + for hit in event.vetoPoint: + detID = hit.GetDetectorID() + pid = hit.PdgCode() + if 1000 < detID < 999999 and abs(pid) == 13: + if muon_vetoPoints.GetSize() == index: + muon_vetoPoints.Expand(index + 1) + muon_vetoPoints[index] = hit + index += 1 + + strawDic, trackDic = [], [] + + # Process straw tube hits + for strawHit in event.strawtubesPoint: + P = r.TMath.Sqrt( + strawHit.GetPx() ** 2 + + strawHit.GetPy() ** 2 + + strawHit.GetPz() ** 2 + ) + if P > 3: + detIDmuonS = ( + strawHit.GetDetectorID() // 10000000 + ) # will need to update if strawtubes detID is changed + if abs(strawHit.PdgCode()) == 13 and detIDmuonS == 1: + # print(detIDmuonS,strawHit.GetTrackID()) + if strawHit.GetTrackID() not in strawDic: + strawDic.append(strawHit.GetTrackID()) + evs[0] += 1 + if evs[0] % 100 == 0: + print(f"evs0: {evs[0]}") + + # Process vetoPoint hits + for hit in event.vetoPoint: + detID = hit.GetDetectorID() + pid = hit.PdgCode() + trackID = hit.GetTrackID() + if 1000 < detID < 999999 and abs(pid) == 13: + if trackID not in trackDic: + trackDic.append(trackID) + evs[1] += 1 + if evs[1] % 100 == 0: + print(f"evs1: {evs[1]}") + + # Check for muons hitting Tr but not SBT and only save them + for m in strawDic: + if m not in trackDic: + for Hit in event.strawtubesPoint: + if m == Hit.GetTrackID() and ev not in processed_events: + processed_events.add(ev) + evs[2] += 1 + print(f"evs2: {evs[2]}") + P = r.TMath.Sqrt( + Hit.GetPx() ** 2 + Hit.GetPy() ** 2 + Hit.GetPz() ** 2 + ) + weight = tree.MCTrack[m].GetWeight() + hPmuon.Fill(P, weight) + + # Fill imuondata + imuondata[0] = float(Hit.PdgCode()) + imuondata[1] = float(Hit.GetPx() / u.GeV) + imuondata[2] = float(Hit.GetPy() / u.GeV) + imuondata[3] = float(Hit.GetPz() / u.GeV) + imuondata[4] = float(Hit.GetX() / u.m) + imuondata[5] = float(Hit.GetY() / u.m) + imuondata[6] = float(Hit.GetZ() / u.m) + imuondata[7] = float(weight) + imuondata[8] = float(hit.GetTime()) + + output_tree.Fill() + + fsel.write( + f"{Hit.PdgCode()} {Hit.GetPx() / u.GeV} {Hit.GetPy() / u.GeV} " + f"{Hit.GetPz() / u.GeV} {Hit.GetX() / u.m} {Hit.GetY() / u.m} " + f"{Hit.GetZ() / u.m} {weight}\n" + ) + +output_file.cd() +output_tree.Write() +hPmuon.Write() +output_file.Close() +print("Finish:", evs) +fsel.close() + + +print( + "------------------------------------------------------file saved, reading", + args.outputfile, + " now----------------------------------------------------------------", +) + +with r.TFile.Open(args.outputfile, "read") as file: + try: + tree = file.MuonAndSoftInteractions + except Exception as e: + print(f"Error: {e}") + exit(1) + + print(f"Processing tree: {tree.GetName()}") + print(f"Total number of entries: {tree.GetEntries()}") + + for event in tree: + imuondata = event.imuondata + pid = imuondata[0] + px = imuondata[1] + py = imuondata[2] + pz = imuondata[3] + x = imuondata[4] + y = imuondata[5] + z = imuondata[6] + weight = imuondata[7] + time_hit = imuondata[8] + num_tracks = len(event.tracks) + num_muonhits = len(event.muon_vetoPoints) + + print( + f"Muon PID: {pid}, x: {x}, y: {y}, z: {z}, t_muon: {time_hit}, " + f"Number of soft tracks in this event: {num_tracks}, " + f"Number of SBT hits from the muon: {num_muonhits}" + ) diff --git a/shipgen/MuDISGenerator.cxx b/shipgen/MuDISGenerator.cxx index e473859c89..c9b5387074 100644 --- a/shipgen/MuDISGenerator.cxx +++ b/shipgen/MuDISGenerator.cxx @@ -12,11 +12,12 @@ #include "TGeoEltu.h" #include "TVectorD.h" #include "TGeoCompositeShape.h" +#include "FairLogger.h" // for FairLogger, etc + -using std::cout; -using std::endl; const Double_t c_light = 29.9792458;//speed of light in cm/ns +const Double_t muon_mass = 0.10565999895334244; //muon mass in GeV // MuDIS momentum GeV // Vertex in SI units, assume this means m @@ -56,16 +57,7 @@ Bool_t MuDISGenerator::Init(const char* fileName, const int firstEvent) { return kFALSE; } - // Print all branches of the TTree - //fTree->Print(); - - // Print out available branches - //TObjArray* branches = fTree->GetListOfBranches(); - //for (Int_t i = 0; i < branches->GetEntries(); i++) { - //TBranch* branch = (TBranch*)branches->At(i); - //cout << "Branch name: " << branch->GetName() << ", Class name: " << branch->GetClassName() << endl; - //} - + fNevents = fTree->GetEntries(); fn = firstEvent; @@ -74,7 +66,7 @@ Bool_t MuDISGenerator::Init(const char* fileName, const int firstEvent) { fTree->SetBranchAddress("DISParticles", &dPart); fTree->SetBranchAddress("SoftParticles", &dPartSoft); // Soft interaction particles - cout << "MuDISGenerator: Initialization successful." << endl; + LOG(INFO) << "MuDISGenerator: Initialization successful."; return kTRUE; } @@ -243,14 +235,14 @@ Bool_t MuDISGenerator::ReadEvent(FairPrimaryGenerator* cpg) fTree->GetEntry(fn%fNevents); fn++; if (fn%10==0) { - cout << "Info MuDISGenerator: MuDIS event-nr "<< fn << endl; + LOG(INFO) << "Info MuDISGenerator: MuDIS event-nr "<< fn; } int nf = dPart->GetEntries(); int ns = dPartSoft->GetEntries(); - //cout << "*********************************************************" << endl; - //cout << "muon DIS Generator debug " << iMuon->GetEntries()<< " "<< iMuon->AddrAt(0) << " nf "<< nf << " fn=" << fn <GetEntries()<< " "<< iMuon->AddrAt(0) << " nf "<< nf << " fn=" << fn; //some start/end positions in z (f.i. emulsion to Tracker 1) Double_t start[3]={0.,0.,startZ}; @@ -258,7 +250,7 @@ Bool_t MuDISGenerator::ReadEvent(FairPrimaryGenerator* cpg) // incoming muon array('d',[pid,px,py,pz,E,x,y,z,w,t]) TVectorD* mu = dynamic_cast(iMuon->AddrAt(0)); - //cout << "muon DIS Generator in muon " << int(mu[0][0])<< endl; + //LOG(DEBUG) << "muon DIS Generator in muon " << int(mu[0][0]); Double_t x = mu[0][5]*100.; // come in m -> cm Double_t y = mu[0][6]*100.; // come in m -> cm Double_t z = mu[0][7]*100.; // come in m -> cm @@ -273,11 +265,11 @@ Bool_t MuDISGenerator::ReadEvent(FairPrimaryGenerator* cpg) start[1]=y-(z-start[2])*tymu; end[0]=x-(z-end[2])*txmu; end[1]=y-(z-end[2])*tymu; - //cout << "MuDIS: mu xyz position " << x << ", " << y << ", " << z << endl; - //cout << "MuDIS: mu pxyz position " << mu[0][1] << ", " << mu[0][2] << ", " << mu[0][3] << endl; - //cout << "MuDIS: mu weight " << w << endl; - //cout << "MuDIS: start position " << start[0] << ", " << start[1] << ", " << start[2] << endl; - //cout << "MuDIS: end position " << end[0] << ", " << end[1] << ", " << end[2] << endl; + //LOG(INFO) << "MuDIS: mu xyz position " << x << ", " << y << ", " << z; + //LOG(INFO) << "MuDIS: mu pxyz position " << mu[0][1] << ", " << mu[0][2] << ", " << mu[0][3]; + //LOG(INFO) << "MuDIS: mu weight " << w; + //LOG(INFO) << "MuDIS: start position " << start[0] << ", " << start[1] << ", " << start[2]; + //LOG(INFO) << "MuDIS: end position " << end[0] << ", " << end[1] << ", " << end[2]; Double_t bparam; Double_t mparam[8]; bparam=MeanMaterialBudget(start, end, mparam); @@ -287,8 +279,8 @@ Bool_t MuDISGenerator::ReadEvent(FairPrimaryGenerator* cpg) Double_t ymu; Double_t zmu; Int_t count=0; - //cout << "Info MuDISGenerator Start prob2int while loop, bparam= " << bparam << ", " << bparam*1.e8 <Uniform(0.,1.)) { @@ -299,36 +291,34 @@ Bool_t MuDISGenerator::ReadEvent(FairPrimaryGenerator* cpg) TGeoNode *node = gGeoManager->FindNode(xmu,ymu,zmu); TGeoMaterial *mat = 0; if (node && !gGeoManager->IsOutside()) mat = node->GetVolume()->GetMaterial(); - //cout << "Info MuDISGenerator: mat " << count << ", " << mat->GetName() << ", " << mat->GetDensity() << endl; + //LOG(DEBUG) << "Info MuDISGenerator: mat " << count << ", " << mat->GetName() << ", " << mat->GetDensity(); //density relative to Prob largest density along this trajectory, i.e. use rho(Pt) prob2int= mat->GetDensity()/mparam[7]; - if (prob2int>1.) cout << "***WARNING*** MuDISGenerator: prob2int > Maximum density????" << prob2int << " maxrho:" << mparam[7] << " material: " << mat->GetName() << endl; + if (prob2int>1.) LOG(WARNING) << "***WARNING*** MuDISGenerator: prob2int > Maximum density????" << prob2int << " maxrho:" << mparam[7] << " material: " << mat->GetName(); count+=1; } - //cout << "Info MuDISGenerator: prob2int " << prob2int << ", " << count << endl; - - //cout << "MuDIS: put position " << xmu << ", " << ymu << ", " << zmu << endl; + + //LOG(DEBUG) << "Info MuDISGenerator: prob2int " << prob2int << ", " << count; + //LOG(DEBUG) << "MuDIS: put position " << xmu << ", " << ymu << ", " << zmu ; - double total_mom = TMath::Sqrt(TMath::Power(mu[0][1], 2) + + Double_t total_mom = TMath::Sqrt(TMath::Power(mu[0][1], 2) + TMath::Power(mu[0][2], 2) + TMath::Power(mu[0][3], 2)); //in GeV - double distance = TMath::Sqrt(TMath::Power(x - xmu, 2) + + Double_t distance = TMath::Sqrt(TMath::Power(x - xmu, 2) + TMath::Power(y - ymu, 2) + TMath::Power(z - zmu, 2)); // in cm - - double mass = 0.10565999895334244; //muon mass in GeV - double v = c_light*total_mom/TMath::Sqrt(TMath::Power(total_mom,2)+ - TMath::Power(mass,2)); - double t_rmu=distance/v; + Double_t v = c_light*total_mom/TMath::Sqrt(TMath::Power(total_mom,2)+ + TMath::Power(muon_mass,2)); + Double_t t_rmu=distance/v; // Adjust time based on the relative positions if (zmu < z) { t_rmu = -t_rmu; } - double t_DIS=(t_muon+t_rmu)/1e9; //time taken in seconds to reach [xmu,ymu,zmu] + Double_t t_DIS=(t_muon+t_rmu)/1e9; //time taken in seconds to reach [xmu,ymu,zmu] cpg->AddTrack(int(mu[0][0]),mu[0][1],mu[0][2],mu[0][3],xmu,ymu,zmu,-1,false,mu[0][4],t_DIS,w);//set time of the track start as t_muon from the input file @@ -338,13 +328,11 @@ Bool_t MuDISGenerator::ReadEvent(FairPrimaryGenerator* cpg) for(int i=0; i(dPart->AddrAt(i)); - //cout << "muon DIS Generator out part " << int(Part[0][0]) << endl; - //cout << "muon DIS Generator out part mom " << Part[0][1]<<" " << Part[0][2] <<" " << Part[0][3] << " " << Part[0][4] << endl; - //cout << "muon DIS Generator out part pos " << mu[0][5]<<" " << mu[0][6] <<" " << mu[0][7] << endl; - //cout << "muon DIS Generator out part w " << mu[0][8] << endl; - cpg->AddTrack(int(Part[0][0]),Part[0][1],Part[0][2],Part[0][3],xmu,ymu,zmu,0,true,Part[0][4],t_DIS,w); //it takes times in seconds. - //cout << "muon DIS part added "<AddTrack(int(Part[0][0]),Part[0][1],Part[0][2],Part[0][3],xmu,ymu,zmu,0,true,Part[0][4],t_DIS,w); } //Soft interaction tracks @@ -354,7 +342,7 @@ Bool_t MuDISGenerator::ReadEvent(FairPrimaryGenerator* cpg) if ( SoftPart[0][7] > zmu ) { //if the softinteractions are after the DIS point, don't save the soft interaction. continue; } - double t_soft=SoftPart[0][8]/1e9; + Double_t t_soft=SoftPart[0][8]/1e9;//time in seconds cpg->AddTrack(int(SoftPart[0][0]), SoftPart[0][1], SoftPart[0][2], SoftPart[0][3], SoftPart[0][5], SoftPart[0][6], SoftPart[0][7], 0, true, SoftPart[0][4], t_soft, w); } diff --git a/shipgen/MuDISGenerator.h b/shipgen/MuDISGenerator.h index ec8ae66136..8e5f8cf01a 100644 --- a/shipgen/MuDISGenerator.h +++ b/shipgen/MuDISGenerator.h @@ -49,6 +49,10 @@ class MuDISGenerator : public FairGenerator int fNevents; int fn; bool fFirst; + + const Double_t c_light = 29.9792458;//speed of light in cm/ns + const Double_t muon_mass = 0.10565999895334244; //muon mass in GeV + ClassDef(MuDISGenerator,1); }; #endif /* !PNDMuGENERATOR_H */ From 908e32a812080c483ae48a7ccbd147d4fcadd67a Mon Sep 17 00:00:00 2001 From: anupama-reghunath Date: Wed, 20 Nov 2024 14:38:15 +0100 Subject: [PATCH 04/10] final fixes --- muonDIS/add_muonresponse.py | 145 +++++---- muonDIS/makeMuonDIS.py | 12 +- shipgen/MuDISGenerator.cxx | 572 +++++++++++++++++++----------------- 3 files changed, 397 insertions(+), 332 deletions(-) diff --git a/muonDIS/add_muonresponse.py b/muonDIS/add_muonresponse.py index 8e98ff932d..cf52191a36 100644 --- a/muonDIS/add_muonresponse.py +++ b/muonDIS/add_muonresponse.py @@ -1,9 +1,9 @@ -#!/usr/bin/env python3 - +#!/usr/bin/env python """Script to add the incoming muon's MC hits in the SBT to the simulation file.""" import argparse import logging +import os import ROOT as r @@ -22,67 +22,100 @@ ) args = parser.parse_args() -# Open the existing ROOT file in update mode -logging.warning(f"{args.inputfile} opened for modification") -inputFile = r.TFile.Open(args.inputfile, "UPDATE") -try: - input_tree = inputFile.cbmsim -except Exception as e: - print(f"Error: {e}") - inputFile.Close() - exit(1) - -# Open the external file with additional vetoPoints -muonFile = r.TFile.Open(args.muonfile, "READ") -try: - muon_tree = muonFile.DIS -except Exception as e: - print(f"Error: {e}") - muonFile.Close() - exit(1) - -inputFile.cd() -output_tree = input_tree.CloneTree( - 0 -) # Clone the structure of the existing tree, but do not copy the entries - -# Access the vetoPoint branch in the original and external trees -originalVetoPoints = r.TClonesArray("vetoPoint") -muonVetoPoints = r.TClonesArray("vetoPoint") -combinedVetoPoints = r.TClonesArray("vetoPoint") - -input_tree.SetBranchAddress("vetoPoint", originalVetoPoints) -muon_tree.SetBranchAddress("muon_vetoPoints", muonVetoPoints) -output_tree.SetBranchAddress("vetoPoint", combinedVetoPoints) - - -for input_event, muon_event in zip(input_tree, muon_tree): - interaction_point = r.TVector3() - input_event.MCTrack[0].GetStartVertex(interaction_point) - combinedVetoPoints.Clear() - - index = 0 +def read_file(): + """Inspecting the produced file for successfully added muon veto points.""" + inputFile = r.TFile.Open(args.inputfile, "READ") + input_tree = inputFile.cbmsim + input_tree.GetEntry(0) - for hit in originalVetoPoints: - if combinedVetoPoints.GetSize() == index: - combinedVetoPoints.Expand(index + 1) - combinedVetoPoints[index] = hit - index += 1 + muons = False + for hit in input_tree.vetoPoint: + if hit.GetTrackID() == 0: + muons = True + inputFile.Close() - for hit in muonVetoPoints: - if hit.GetZ() < interaction_point.Z(): + if muons: + print("File is already updated with incoming muon info. Nothing to do.") + return False + else: + print( + "Incoming muon's vetopoint info missing in file, proceeding with modification" + ) + return True + + +def modify_file(): + """Modify the input file with muon info.""" + if args.inputfile[0:4] == "/eos": + args.inputfile = os.environ["EOSSHIP"] + args.inputfile + + logging.warning(f"{args.inputfile} opened for modification") + + inputFile = r.TFile.Open(args.inputfile, "UPDATE") + try: + input_tree = inputFile.cbmsim + except Exception as e: + print(f"Error: {e}") + inputFile.Close() + exit(1) + + # Open the external file with additional vetoPoints + muonFile = r.TFile.Open(args.muonfile, "READ") + try: + muon_tree = muonFile.DIS + except Exception as e: + print(f"Error: {e}") + muonFile.Close() + exit(1) + + inputFile.cd() + output_tree = input_tree.CloneTree( + 0 + ) # Clone the structure of the existing tree, but do not copy the entries + + # Access the vetoPoint branch in the original and external trees + originalVetoPoints = r.TClonesArray("vetoPoint") + muonVetoPoints = r.TClonesArray("vetoPoint") + combinedVetoPoints = r.TClonesArray("vetoPoint") + + input_tree.SetBranchAddress("vetoPoint", originalVetoPoints) + muon_tree.SetBranchAddress("muon_vetoPoints", muonVetoPoints) + output_tree.SetBranchAddress("vetoPoint", combinedVetoPoints) + + for input_event, muon_event in zip(input_tree, muon_tree): + interaction_point = r.TVector3() + input_event.MCTrack[0].GetStartVertex(interaction_point) + + combinedVetoPoints.Clear() + + index = 0 + + for hit in originalVetoPoints: if combinedVetoPoints.GetSize() == index: combinedVetoPoints.Expand(index + 1) combinedVetoPoints[index] = hit index += 1 - output_tree.Fill() + for hit in muonVetoPoints: + if hit.GetZ() < interaction_point.Z(): + if combinedVetoPoints.GetSize() == index: + combinedVetoPoints.Expand(index + 1) + combinedVetoPoints[index] = hit + index += 1 + + output_tree.Fill() + + # Write the updated tree back to the same file + inputFile.cd() + output_tree.Write("cbmsim", r.TObject.kOverwrite) + inputFile.Close() + muonFile.Close() + + print(f"File updated with incoming muon info as {args.inputfile}") + -# Write the updated tree back to the same file -inputFile.cd() # Ensure we are in the correct directory of the ROOT file -output_tree.Write("cbmsim", r.TObject.kOverwrite) -inputFile.Close() -muonFile.Close() +add_muons = read_file() -print(f"Updated file saved as {args.inputfile}") +if add_muons: + modify_file() diff --git a/muonDIS/makeMuonDIS.py b/muonDIS/makeMuonDIS.py index 69fa271d23..f659678417 100644 --- a/muonDIS/makeMuonDIS.py +++ b/muonDIS/makeMuonDIS.py @@ -94,7 +94,7 @@ def makeMuonDIS(): mutype = {-13: "gamma/mu+", 13: "gamma/mu-"} myPythia.SetMSTU(11, 11) - logging.info(f"Processing events from {nStart} to {nEnd}...") + logging.info(f"Processing events from {nStart} to {nEnd-1}...") nMade = 0 @@ -121,7 +121,7 @@ def makeMuonDIS(): phi = r.TMath.ATan2(py, px) logging.debug( - f"Muon: PID={pid}, px={px}, py={py}, pz={pz}, E={E}, x={x}, y={y}, z={z}, w={w}" + f"\n\tMuon index:{k} \n\tPID = {pid}, weight = {w} \n\tpx = {px}, py = {py}, pz = {pz}, E = {E},\n\tx = {x}, y = {y}, z = {z}\n" ) isProton = 1 @@ -136,7 +136,7 @@ def makeMuonDIS(): if a == args.nDIS // 2: myPythia.Initialize("FIXT", mutype[pid], "n0", p) # target = "n0" isProton = 0 - logging.info("Switching to neutron interaction") + logging.debug("Switching to neutron interaction") dPartDIS.Clear() iMuon.Clear() @@ -144,7 +144,9 @@ def makeMuonDIS(): iMuon[0] = muPart myPythia.GenerateEvent() myPythia.Pyedit(1) - logging.info(f"Event {a} generated, number of particles: {myPythia.GetN()}") + logging.debug( + f"DIS Event {a} generated, number of particles: {myPythia.GetN()}" + ) for itrk in range(1, myPythia.GetN() + 1): xsec = myPythia.GetPARI(1) @@ -218,7 +220,7 @@ def makeMuonDIS(): output_tree.Write() myPythia.SetMSTU(11, 6) logging.info( - f"Created DIS for saved in muonDis_{args.nJob}.root, muon event index from {nStart} - {nEnd} , nDISPerMuon {args.nDIS}" + f"DIS generated for muons (index {nStart} - {nEnd-1}) , output saved in muonDis_{args.nJob}.root, nDISPerMuon = {args.nDIS}" ) diff --git a/shipgen/MuDISGenerator.cxx b/shipgen/MuDISGenerator.cxx index c9b5387074..f64dc20167 100644 --- a/shipgen/MuDISGenerator.cxx +++ b/shipgen/MuDISGenerator.cxx @@ -1,357 +1,387 @@ -#include -#include "TSystem.h" -#include "TROOT.h" -#include "TMath.h" -#include "TFile.h" -#include "TRandom.h" -#include "FairPrimaryGenerator.h" #include "MuDISGenerator.h" -#include "TGeoVolume.h" -#include "TGeoNode.h" -#include "TGeoManager.h" + +#include "FairLogger.h" +#include "FairPrimaryGenerator.h" +#include "TFile.h" +#include "TGeoCompositeShape.h" #include "TGeoEltu.h" +#include "TGeoManager.h" +#include "TGeoNode.h" +#include "TGeoVolume.h" +#include "TMath.h" +#include "TROOT.h" +#include "TRandom.h" +#include "TSystem.h" #include "TVectorD.h" -#include "TGeoCompositeShape.h" -#include "FairLogger.h" // for FairLogger, etc - - -const Double_t c_light = 29.9792458;//speed of light in cm/ns -const Double_t muon_mass = 0.10565999895334244; //muon mass in GeV +#include // MuDIS momentum GeV // Vertex in SI units, assume this means m - // ----- Default constructor ------------------------------------------- MuDISGenerator::MuDISGenerator() {} // ------------------------------------------------------------------------- // ----- Default constructor ------------------------------------------- -Bool_t MuDISGenerator::Init(const char* fileName) { - return Init(fileName, 0); +Bool_t MuDISGenerator::Init(const char* fileName) +{ + return Init(fileName, 0); } -Bool_t MuDISGenerator::Init(const char* fileName, const int firstEvent) { +Bool_t MuDISGenerator::Init(const char* fileName, const int firstEvent) +{ LOGF(info, "Opening input file %s", fileName); iMuon = 0; - dPart = 0; - dPartSoft = 0; + dPart = 0; + dPartSoft = 0; - if (0 == strncmp("/eos",fileName,4)) { + if (0 == strncmp("/eos", fileName, 4)) { TString tmp = gSystem->Getenv("EOSSHIP"); tmp += fileName; - fInputFile = TFile::Open(tmp); + fInputFile = TFile::Open(tmp, "read"); LOGF(info, "Open external file on eos: %s", tmp.Data()); } else { - fInputFile = new TFile(fileName); + fInputFile = TFile::Open(fileName, "read"); } - if (!fInputFile || fInputFile->IsZombie()) { - LOG(FATAL) << "Error opening input file"; - return kFALSE; - } - - fTree = (TTree *)fInputFile->Get("DIS"); + fTree = dynamic_cast(fInputFile->Get("DIS")); if (!fTree) { LOG(FATAL) << "Error: Tree 'DIS' not found in file"; return kFALSE; } - fNevents = fTree->GetEntries(); fn = firstEvent; // Set branch addresses - fTree->SetBranchAddress("InMuon", &iMuon); // incoming muon + fTree->SetBranchAddress("InMuon", &iMuon); // incoming muon fTree->SetBranchAddress("DISParticles", &dPart); - fTree->SetBranchAddress("SoftParticles", &dPartSoft); // Soft interaction particles + fTree->SetBranchAddress("SoftParticles", &dPartSoft); // Soft interaction particles LOG(INFO) << "MuDISGenerator: Initialization successful."; return kTRUE; } -Double_t MuDISGenerator::MeanMaterialBudget(const Double_t *start, const Double_t *end, Double_t *mparam) +Double_t MuDISGenerator::MeanMaterialBudget(const Double_t* start, const Double_t* end, Double_t* mparam) { - // - // Calculate mean material budget and material properties between - // the points "start" and "end". - // - // "mparam" - parameters used for the energy and multiple scattering - // corrections: - // - // mparam[0] - mean density: sum(x_i*rho_i)/sum(x_i) [g/cm3] - // mparam[1] - equivalent rad length fraction: sum(x_i/X0_i) [adimensional] - // mparam[2] - mean A: sum(x_i*A_i)/sum(x_i) [adimensional] - // mparam[3] - mean Z: sum(x_i*Z_i)/sum(x_i) [adimensional] - // mparam[4] - length: sum(x_i) [cm] - // mparam[5] - Z/A mean: sum(x_i*Z_i/A_i)/sum(x_i) [adimensional] - // mparam[6] - number of boundary crosses - // mparam[7] - maximum density encountered (g/cm^3) - // - // Origin: Marian Ivanov, Marian.Ivanov@cern.ch - // - // Corrections and improvements by - // Andrea Dainese, Andrea.Dainese@lnl.infn.it, - // Andrei Gheata, Andrei.Gheata@cern.ch - // Anupama Reghunath, anupama.reghunath@cern.ch - // - - mparam[0]=0; mparam[1]=1; mparam[2] =0; mparam[3] =0; - mparam[4]=0; mparam[5]=0; mparam[6]=0; mparam[7]=0; - // - Double_t bparam[6]; // total parameters - Double_t lparam[6]; // local parameters - - for (Int_t i=0;i<6;i++) bparam[i]=0; - - if (!gGeoManager) { - //AliFatalClass("No TGeo\n"); - return 0.; - } - // - Double_t length; - Double_t dir[3]; - length = TMath::Sqrt((end[0]-start[0])*(end[0]-start[0])+ - (end[1]-start[1])*(end[1]-start[1])+ - (end[2]-start[2])*(end[2]-start[2])); - mparam[4]=length; - if (lengthInitTrack(start, dir); - if (!startnode) { - //AliErrorClass(Form("start point out of geometry: x %f, y %f, z %f", - // start[0],start[1],start[2])); - return 0.0; - } - TGeoMaterial *material = startnode->GetVolume()->GetMedium()->GetMaterial(); - lparam[0] = material->GetDensity(); - if (lparam[0] > mparam[7]) mparam[7]=lparam[0]; - lparam[1] = material->GetRadLen(); - lparam[2] = material->GetA(); - lparam[3] = material->GetZ(); - lparam[4] = length; - lparam[5] = lparam[3]/lparam[2]; - if (material->IsMixture()) { - TGeoMixture * mixture = (TGeoMixture*)material; - lparam[5] =0; - Double_t sum =0; - for (Int_t iel=0;ielGetNelements();iel++){ - sum += mixture->GetWmixt()[iel]; - lparam[5]+= mixture->GetZmixt()[iel]*mixture->GetWmixt()[iel]/mixture->GetAmixt()[iel]; + // + // Calculate mean material budget and material properties between + // the points "start" and "end". + // + // "mparam" - parameters used for the energy and multiple scattering + // corrections: + // + // mparam[0] - mean density: sum(x_i*rho_i)/sum(x_i) [g/cm3] + // mparam[1] - equivalent rad length fraction: sum(x_i/X0_i) [adimensional] + // mparam[2] - mean A: sum(x_i*A_i)/sum(x_i) [adimensional] + // mparam[3] - mean Z: sum(x_i*Z_i)/sum(x_i) [adimensional] + // mparam[4] - length: sum(x_i) [cm] + // mparam[5] - Z/A mean: sum(x_i*Z_i/A_i)/sum(x_i) [adimensional] + // mparam[6] - number of boundary crosses + // mparam[7] - maximum density encountered (g/cm^3) + // + // Origin: Marian Ivanov, Marian.Ivanov@cern.ch + // + // Corrections and improvements by + // Andrea Dainese, Andrea.Dainese@lnl.infn.it, + // Andrei Gheata, Andrei.Gheata@cern.ch + // Anupama Reghunath, anupama.reghunath@cern.ch + // + + mparam[0] = 0; + mparam[1] = 1; + mparam[2] = 0; + mparam[3] = 0; + mparam[4] = 0; + mparam[5] = 0; + mparam[6] = 0; + mparam[7] = 0; + // + Double_t bparam[6]; // total parameters + Double_t lparam[6]; // local parameters + + for (Int_t i = 0; i < 6; i++) + bparam[i] = 0; + + if (!gGeoManager) { + // AliFatalClass("No TGeo\n"); + return 0.; } - lparam[5]/=sum; - } - - // Locate next boundary within length without computing safety. - // Propagate either with length (if no boundary found) or just cross boundary - gGeoManager->FindNextBoundaryAndStep(length, kFALSE); - Double_t step = 0.0; // Step made - Double_t snext = gGeoManager->GetStep(); - // If no boundary within proposed length, return current density - if (!gGeoManager->IsOnBoundary()) { - mparam[0] = lparam[0]; - mparam[1] = lparam[4]/lparam[1]; - mparam[2] = lparam[2]; - mparam[3] = lparam[3]; - mparam[4] = lparam[4]; - return lparam[0]; - } - // Try to cross the boundary and see what is next - Int_t nzero = 0; - while (length>TGeoShape::Tolerance()) { - currentnode = gGeoManager->GetCurrentNode(); - if (snext<2.*TGeoShape::Tolerance()) nzero++; - else nzero = 0; - if (nzero>3) { - // This means navigation has problems on one boundary - // Try to cross by making a small step - //AliErrorClass("Cannot cross boundary\n"); - mparam[0] = bparam[0]/step; - mparam[1] = bparam[1]; - mparam[2] = bparam[2]/step; - mparam[3] = bparam[3]/step; - mparam[5] = bparam[5]/step; - mparam[4] = step; - mparam[0] = 0.; // if crash of navigation take mean density 0 - mparam[1] = 1000000; // and infinite rad length - return bparam[0]/step; + // + Double_t length; + Double_t dir[3]; + length = TMath::Sqrt((end[0] - start[0]) * (end[0] - start[0]) + (end[1] - start[1]) * (end[1] - start[1]) + + (end[2] - start[2]) * (end[2] - start[2])); + mparam[4] = length; + if (length < TGeoShape::Tolerance()) + return 0.0; + Double_t invlen = 1. / length; + dir[0] = (end[0] - start[0]) * invlen; + dir[1] = (end[1] - start[1]) * invlen; + dir[2] = (end[2] - start[2]) * invlen; + + // Initialize start point and direction + TGeoNode* currentnode = 0; + TGeoNode* startnode = gGeoManager->InitTrack(start, dir); + if (!startnode) { + // AliErrorClass(Form("start point out of geometry: x %f, y %f, z %f", + // start[0],start[1],start[2])); + return 0.0; } - mparam[6]+=1.; - step += snext; - bparam[1] += snext/lparam[1]; - bparam[2] += snext*lparam[2]; - bparam[3] += snext*lparam[3]; - bparam[5] += snext*lparam[5]; - bparam[0] += snext*lparam[0]; - - if (snext>=length) break; - if (!currentnode) break; - length -= snext; - material = currentnode->GetVolume()->GetMedium()->GetMaterial(); + TGeoMaterial* material = startnode->GetVolume()->GetMedium()->GetMaterial(); lparam[0] = material->GetDensity(); - if (lparam[0] > mparam[7]) mparam[7]=lparam[0]; - lparam[1] = material->GetRadLen(); - lparam[2] = material->GetA(); - lparam[3] = material->GetZ(); - lparam[5] = lparam[3]/lparam[2]; + if (lparam[0] > mparam[7]) + mparam[7] = lparam[0]; + lparam[1] = material->GetRadLen(); + lparam[2] = material->GetA(); + lparam[3] = material->GetZ(); + lparam[4] = length; + lparam[5] = lparam[3] / lparam[2]; if (material->IsMixture()) { - TGeoMixture * mixture = (TGeoMixture*)material; - lparam[5]=0; - Double_t sum =0; - for (Int_t iel=0;ielGetNelements();iel++){ - sum+= mixture->GetWmixt()[iel]; - lparam[5]+= mixture->GetZmixt()[iel]*mixture->GetWmixt()[iel]/mixture->GetAmixt()[iel]; - } - lparam[5]/=sum; + TGeoMixture* mixture = (TGeoMixture*)material; + lparam[5] = 0; + Double_t sum = 0; + for (Int_t iel = 0; iel < mixture->GetNelements(); iel++) { + sum += mixture->GetWmixt()[iel]; + lparam[5] += mixture->GetZmixt()[iel] * mixture->GetWmixt()[iel] / mixture->GetAmixt()[iel]; + } + lparam[5] /= sum; } + + // Locate next boundary within length without computing safety. + // Propagate either with length (if no boundary found) or just cross boundary gGeoManager->FindNextBoundaryAndStep(length, kFALSE); - snext = gGeoManager->GetStep(); - } - mparam[0] = bparam[0]/step; - mparam[1] = bparam[1]; - mparam[2] = bparam[2]/step; - mparam[3] = bparam[3]/step; - mparam[5] = bparam[5]/step; - return bparam[0]/step; + Double_t step = 0.0; // Step made + Double_t snext = gGeoManager->GetStep(); + // If no boundary within proposed length, return current density + if (!gGeoManager->IsOnBoundary()) { + mparam[0] = lparam[0]; + mparam[1] = lparam[4] / lparam[1]; + mparam[2] = lparam[2]; + mparam[3] = lparam[3]; + mparam[4] = lparam[4]; + return lparam[0]; + } + // Try to cross the boundary and see what is next + Int_t nzero = 0; + while (length > TGeoShape::Tolerance()) { + currentnode = gGeoManager->GetCurrentNode(); + if (snext < 2. * TGeoShape::Tolerance()) + nzero++; + else + nzero = 0; + if (nzero > 3) { + // This means navigation has problems on one boundary + // Try to cross by making a small step + // AliErrorClass("Cannot cross boundary\n"); + mparam[0] = bparam[0] / step; + mparam[1] = bparam[1]; + mparam[2] = bparam[2] / step; + mparam[3] = bparam[3] / step; + mparam[5] = bparam[5] / step; + mparam[4] = step; + mparam[0] = 0.; // if crash of navigation take mean density 0 + mparam[1] = 1000000; // and infinite rad length + return bparam[0] / step; + } + mparam[6] += 1.; + step += snext; + bparam[1] += snext / lparam[1]; + bparam[2] += snext * lparam[2]; + bparam[3] += snext * lparam[3]; + bparam[5] += snext * lparam[5]; + bparam[0] += snext * lparam[0]; + + if (snext >= length) + break; + if (!currentnode) + break; + length -= snext; + material = currentnode->GetVolume()->GetMedium()->GetMaterial(); + lparam[0] = material->GetDensity(); + if (lparam[0] > mparam[7]) + mparam[7] = lparam[0]; + lparam[1] = material->GetRadLen(); + lparam[2] = material->GetA(); + lparam[3] = material->GetZ(); + lparam[5] = lparam[3] / lparam[2]; + if (material->IsMixture()) { + TGeoMixture* mixture = (TGeoMixture*)material; + lparam[5] = 0; + Double_t sum = 0; + for (Int_t iel = 0; iel < mixture->GetNelements(); iel++) { + sum += mixture->GetWmixt()[iel]; + lparam[5] += mixture->GetZmixt()[iel] * mixture->GetWmixt()[iel] / mixture->GetAmixt()[iel]; + } + lparam[5] /= sum; + } + gGeoManager->FindNextBoundaryAndStep(length, kFALSE); + snext = gGeoManager->GetStep(); + } + mparam[0] = bparam[0] / step; + mparam[1] = bparam[1]; + mparam[2] = bparam[2] / step; + mparam[3] = bparam[3] / step; + mparam[5] = bparam[5] / step; + return bparam[0] / step; } // ----- Destructor ---------------------------------------------------- MuDISGenerator::~MuDISGenerator() { - fInputFile->Close(); - fInputFile->Delete(); - delete fInputFile; + fInputFile->Close(); + fInputFile->Delete(); + delete fInputFile; } // ----- Passing the event --------------------------------------------- Bool_t MuDISGenerator::ReadEvent(FairPrimaryGenerator* cpg) { - if (fn==fNevents) {LOG(WARNING) << "End of input file. Rewind.";} - fTree->GetEntry(fn%fNevents); + if (fn == fNevents) { + LOG(WARNING) << "End of input file. Rewind."; + } + fTree->GetEntry(fn % fNevents); fn++; - if (fn%10==0) { - LOG(INFO) << "Info MuDISGenerator: MuDIS event-nr "<< fn; + if (fn % 10 == 0) { + LOG(INFO) << "Info MuDISGenerator: MuDIS event-nr " << fn; } - int nf = dPart->GetEntries(); - int ns = dPartSoft->GetEntries(); - //LOG(DEBUG) << "*********************************************************" ; - //LOG(DEBUG) << "muon DIS Generator debug " << iMuon->GetEntries()<< " "<< iMuon->AddrAt(0) << " nf "<< nf << " fn=" << fn; + int ns = dPartSoft->GetEntries(); + LOG(DEBUG) << "*********************************************************"; + LOG(DEBUG) << "muon DIS Generator debug " << iMuon->GetEntries() << " " << iMuon->AddrAt(0) << " nf " << nf + << " fn=" << fn; - //some start/end positions in z (f.i. emulsion to Tracker 1) - Double_t start[3]={0.,0.,startZ}; - Double_t end[3]={0.,0.,endZ}; + // some start/end positions in z (f.i. emulsion to Tracker 1) + Double_t start[3] = {0., 0., startZ}; + Double_t end[3] = {0., 0., endZ}; -// incoming muon array('d',[pid,px,py,pz,E,x,y,z,w,t]) + // incoming muon array('d',[pid,px,py,pz,E,x,y,z,w,t]) TVectorD* mu = dynamic_cast(iMuon->AddrAt(0)); - //LOG(DEBUG) << "muon DIS Generator in muon " << int(mu[0][0]); - Double_t x = mu[0][5]*100.; // come in m -> cm - Double_t y = mu[0][6]*100.; // come in m -> cm - Double_t z = mu[0][7]*100.; // come in m -> cm + LOG(DEBUG) << "muon DIS Generator in muon " << int(mu[0][0]); + Double_t x = mu[0][5] * 100.; // come in m -> cm + Double_t y = mu[0][6] * 100.; // come in m -> cm + Double_t z = mu[0][7] * 100.; // come in m -> cm Double_t w = mu[0][8]; - Double_t t_muon = mu[0][11]; // in ns - -// calculate start/end positions along this muon, and amout of material in between - - Double_t txmu=mu[0][1]/mu[0][3]; - Double_t tymu=mu[0][2]/mu[0][3]; - start[0]=x-(z-start[2])*txmu; - start[1]=y-(z-start[2])*tymu; - end[0]=x-(z-end[2])*txmu; - end[1]=y-(z-end[2])*tymu; - //LOG(INFO) << "MuDIS: mu xyz position " << x << ", " << y << ", " << z; - //LOG(INFO) << "MuDIS: mu pxyz position " << mu[0][1] << ", " << mu[0][2] << ", " << mu[0][3]; - //LOG(INFO) << "MuDIS: mu weight " << w; - //LOG(INFO) << "MuDIS: start position " << start[0] << ", " << start[1] << ", " << start[2]; - //LOG(INFO) << "MuDIS: end position " << end[0] << ", " << end[1] << ", " << end[2]; + Double_t t_muon = mu[0][11]; // in ns + + // calculate start/end positions along this muon, and amout of material in between + + Double_t txmu = mu[0][1] / mu[0][3]; + Double_t tymu = mu[0][2] / mu[0][3]; + start[0] = x - (z - start[2]) * txmu; + start[1] = y - (z - start[2]) * tymu; + end[0] = x - (z - end[2]) * txmu; + end[1] = y - (z - end[2]) * tymu; + LOG(DEBUG) << "MuDIS: mu xyz position " << x << ", " << y << ", " << z; + LOG(DEBUG) << "MuDIS: mu pxyz position " << mu[0][1] << ", " << mu[0][2] << ", " << mu[0][3]; + LOG(DEBUG) << "MuDIS: mu weight " << w; + LOG(DEBUG) << "MuDIS: start position " << start[0] << ", " << start[1] << ", " << start[2]; + LOG(DEBUG) << "MuDIS: end position " << end[0] << ", " << end[1] << ", " << end[2]; Double_t bparam; Double_t mparam[8]; - bparam=MeanMaterialBudget(start, end, mparam); - //loop over trajectory between start and end to pick an interaction point + bparam = MeanMaterialBudget(start, end, mparam); + // loop over trajectory between start and end to pick an interaction point Double_t prob2int = 0.; Double_t xmu; Double_t ymu; Double_t zmu; - Int_t count=0; - //LOG(INFO) << "Info MuDISGenerator Start prob2int while loop, bparam= " << bparam << ", " << bparam*1.e8 ; - //LOG(INFO) << "Info MuDISGenerator What was maximum density, mparam[7]= " << mparam[7] << ", " << mparam[7]*1.e8; - - - while (prob2intUniform(0.,1.)) { - zmu=gRandom->Uniform(start[2],end[2]); - xmu=x-(z-zmu)*txmu; - ymu=y-(z-zmu)*tymu; - //get local material at this point - TGeoNode *node = gGeoManager->FindNode(xmu,ymu,zmu); - TGeoMaterial *mat = 0; - if (node && !gGeoManager->IsOutside()) mat = node->GetVolume()->GetMaterial(); - //LOG(DEBUG) << "Info MuDISGenerator: mat " << count << ", " << mat->GetName() << ", " << mat->GetDensity(); - //density relative to Prob largest density along this trajectory, i.e. use rho(Pt) - prob2int= mat->GetDensity()/mparam[7]; - if (prob2int>1.) LOG(WARNING) << "***WARNING*** MuDISGenerator: prob2int > Maximum density????" << prob2int << " maxrho:" << mparam[7] << " material: " << mat->GetName(); - count+=1; + Int_t count = 0; + LOG(DEBUG) << "Info MuDISGenerator Start prob2int while loop, bparam= " << bparam << ", " << bparam * 1.e8; + LOG(DEBUG) << "Info MuDISGenerator What was maximum density, mparam[7]= " << mparam[7] << ", " << mparam[7] * 1.e8; + + while (prob2int < gRandom->Uniform(0., 1.)) { + zmu = gRandom->Uniform(start[2], end[2]); + xmu = x - (z - zmu) * txmu; + ymu = y - (z - zmu) * tymu; + // get local material at this point + TGeoNode* node = gGeoManager->FindNode(xmu, ymu, zmu); + TGeoMaterial* mat = 0; + if (node && !gGeoManager->IsOutside()) + mat = node->GetVolume()->GetMaterial(); + LOG(DEBUG) << "Info MuDISGenerator: mat " << count << ", " << mat->GetName() << ", " << mat->GetDensity(); + // density relative to Prob largest density along this trajectory, i.e. use rho(Pt) + prob2int = mat->GetDensity() / mparam[7]; + if (prob2int > 1.) { + LOG(WARNING) << "MuDISGenerator: prob2int > Maximum density????"; + LOG(WARNING) << "prob2int: " << prob2int; + LOG(WARNING) << "maxrho: " << mparam[7]; + LOG(WARNING) << "material: " << mat->GetName(); + } + count += 1; } - - //LOG(DEBUG) << "Info MuDISGenerator: prob2int " << prob2int << ", " << count; - //LOG(DEBUG) << "MuDIS: put position " << xmu << ", " << ymu << ", " << zmu ; - Double_t total_mom = TMath::Sqrt(TMath::Power(mu[0][1], 2) + - TMath::Power(mu[0][2], 2) + - TMath::Power(mu[0][3], 2)); //in GeV + LOG(DEBUG) << "Info MuDISGenerator: prob2int " << prob2int << ", " << count; + LOG(DEBUG) << "MuDIS: put position " << xmu << ", " << ymu << ", " << zmu; - Double_t distance = TMath::Sqrt(TMath::Power(x - xmu, 2) + - TMath::Power(y - ymu, 2) + - TMath::Power(z - zmu, 2)); // in cm + Double_t total_mom = + TMath::Sqrt(TMath::Power(mu[0][1], 2) + TMath::Power(mu[0][2], 2) + TMath::Power(mu[0][3], 2)); // in GeV - Double_t v = c_light*total_mom/TMath::Sqrt(TMath::Power(total_mom,2)+ - TMath::Power(muon_mass,2)); - Double_t t_rmu=distance/v; + Double_t distance = + TMath::Sqrt(TMath::Power(x - xmu, 2) + TMath::Power(y - ymu, 2) + TMath::Power(z - zmu, 2)); // in cm + + Double_t v = c_light * total_mom / TMath::Sqrt(TMath::Power(total_mom, 2) + TMath::Power(muon_mass, 2)); + Double_t t_rmu = distance / v; // Adjust time based on the relative positions if (zmu < z) { t_rmu = -t_rmu; } - Double_t t_DIS=(t_muon+t_rmu)/1e9; //time taken in seconds to reach [xmu,ymu,zmu] - - cpg->AddTrack(int(mu[0][0]),mu[0][1],mu[0][2],mu[0][3],xmu,ymu,zmu,-1,false,mu[0][4],t_DIS,w);//set time of the track start as t_muon from the input file - - // outgoing particles, [did,dpx,dpy,dpz,E], put density along trajectory as weight, g/cm^2 - w=mparam[0]*mparam[4];//modify weight, by multiplying with average densiy along track?? - - for(int i=0; i(dPart->AddrAt(i)); - //LOG(DEBUG)<< "muon DIS Generator out part " << int(Part[0][0]); - //LOG(DEBUG) << "muon DIS Generator out part mom " << Part[0][1]<<" " << Part[0][2] <<" " << Part[0][3] << " " << Part[0][4]; - //LOG(DEBUG) << "muon DIS Generator out part pos " << mu[0][5]<<" " << mu[0][6] <<" " << mu[0][7]; - //LOG(DEBUG) << "muon DIS Generator out part w " << mu[0][8] ; - cpg->AddTrack(int(Part[0][0]),Part[0][1],Part[0][2],Part[0][3],xmu,ymu,zmu,0,true,Part[0][4],t_DIS,w); } - - //Soft interaction tracks - - - for (int i = 0; i < ns; i++) { - TVectorD* SoftPart = dynamic_cast(dPartSoft->AddrAt(i)); - if ( SoftPart[0][7] > zmu ) { //if the softinteractions are after the DIS point, don't save the soft interaction. - continue; - } - Double_t t_soft=SoftPart[0][8]/1e9;//time in seconds - cpg->AddTrack(int(SoftPart[0][0]), SoftPart[0][1], SoftPart[0][2], SoftPart[0][3], SoftPart[0][5], SoftPart[0][6], SoftPart[0][7], 0, true, SoftPart[0][4], t_soft, w); + Double_t t_DIS = (t_muon + t_rmu) / 1e9; // time taken in seconds to reach [xmu,ymu,zmu] + + cpg->AddTrack(int(mu[0][0]), // incoming muon track () + mu[0][1], + mu[0][2], + mu[0][3], + xmu, + ymu, + zmu, + -1, + false, // tracking disabled + mu[0][4], + t_DIS, // shift time of the incoming muon track wrt t_muon from the input file. + w); // weight associated with a spill. + + // outgoing DIS particles, [did,dpx,dpy,dpz,E], put density along trajectory as weight, g/cm^2 + + w = mparam[0] * mparam[4]; // modify weight, by multiplying with average densiy * track length + + for (auto&& particle : *dPart) { + TVectorD* Part = dynamic_cast(particle); + LOG(DEBUG) << "muon DIS Generator out part " << int((*Part)[0]); + LOG(DEBUG) << "muon DIS Generator out part mom " << (*Part)[1] << " " << (*Part)[2] << " " << (*Part)[3] << " " + << (*Part)[4]; + LOG(DEBUG) << "muon DIS Generator out part pos " << xmu << " " << ymu << "" << zmu; + LOG(DEBUG) << "muon DIS Generator out part w " << w; + cpg->AddTrack( + int((*Part)[0]), (*Part)[1], (*Part)[2], (*Part)[3], xmu, ymu, zmu, 0, true, (*Part)[4], t_DIS, w); + } + // Soft interaction tracks + for (auto&& softParticle : *dPartSoft) { + TVectorD* SoftPart = dynamic_cast(softParticle); + if ((*SoftPart)[7] > zmu) { + continue; + } // Soft interactions after the DIS point are not saved + Double_t t_soft = (*SoftPart)[8] / 1e9; // Time in seconds + cpg->AddTrack(int((*SoftPart)[0]), + (*SoftPart)[1], + (*SoftPart)[2], + (*SoftPart)[3], + (*SoftPart)[5], + (*SoftPart)[6], + (*SoftPart)[7], + 0, + true, + (*SoftPart)[4], + t_soft, + w); } - return kTRUE; + return kTRUE; } // ------------------------------------------------------------------------- Int_t MuDISGenerator::GetNevents() { - return fNevents; + return fNevents; } - From 8060675c3f422e7a9d00e7199b6cce087c569b41 Mon Sep 17 00:00:00 2001 From: anupama-reghunath Date: Wed, 20 Nov 2024 17:41:39 +0100 Subject: [PATCH 05/10] add CHANGELOG entry for MuonDIS --- CHANGELOG.md | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 373fc03bfa..66d4353110 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -16,7 +16,9 @@ it in future. ### Fixed -* Use ConstructedAt instead of remove pythonization for TClonesArray +* Use ConstructedAt instead of remove pythonization for TClonesArray +* **Corrections in MuonDIS simulation** + The DIS interactions are now time-shifted to be consistent with the original incoming muon. Additionally, tracks from soft interactions of the original muon along with the muon's veto response are preserved (in muonDis_<>.root) and included up to the DIS interaction point. To be noted that the muon veto points are manually added using add_muonresponse.py, which modifies the simulation file. This replaces the old method of "backward-travelling muon" to generate the incoming muon's veto response. All MuonDIS simulation scripts have been updated and consolidated within FairShip/muonDIS, ensuring consistency for new productions. ### Changed From 7a1a5ff85049b44913226f5b810767bd053bf489 Mon Sep 17 00:00:00 2001 From: anupama-reghunath Date: Tue, 26 Nov 2024 14:09:06 +0100 Subject: [PATCH 06/10] Avoid listdir for non-directory case --- muonDIS/make_nTuple_SBT.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/muonDIS/make_nTuple_SBT.py b/muonDIS/make_nTuple_SBT.py index 86e6f7317f..416b08b303 100644 --- a/muonDIS/make_nTuple_SBT.py +++ b/muonDIS/make_nTuple_SBT.py @@ -65,6 +65,8 @@ logging.basicConfig(level=logging.INFO) for inputFolder in os.listdir(path): + if not os.path.isdir(os.path.join(path, inputFolder)): + continue for subFolder in os.listdir(os.path.join(path, inputFolder)): if not os.path.isdir(os.path.join(path, inputFolder, subFolder)): continue From cd92c6e8d2d56d8d25cbb19f8a453952e7d4101d Mon Sep 17 00:00:00 2001 From: anupama-reghunath Date: Tue, 26 Nov 2024 14:13:31 +0100 Subject: [PATCH 07/10] fixup! Avoid listdir for non-directory case --- muonDIS/make_nTuple_Tr.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/muonDIS/make_nTuple_Tr.py b/muonDIS/make_nTuple_Tr.py index 3e04d4f59c..85e41fe126 100644 --- a/muonDIS/make_nTuple_Tr.py +++ b/muonDIS/make_nTuple_Tr.py @@ -67,6 +67,8 @@ logging.basicConfig(level=logging.INFO) for inputFolder in os.listdir(path): + if not os.path.isdir(os.path.join(path, inputFolder)): + continue for subFolder in os.listdir(os.path.join(path, inputFolder)): if not os.path.isdir(os.path.join(path, inputFolder, subFolder)): continue From c029ddfe63a82cc1848381a7225c4661173a9fb7 Mon Sep 17 00:00:00 2001 From: anupama-reghunath Date: Wed, 27 Nov 2024 14:30:34 +0100 Subject: [PATCH 08/10] fixup! final fixes --- muonDIS/add_muonresponse.py | 74 ++++++++++++++++++------------------- muonDIS/makeMuonDIS.py | 4 +- muonDIS/make_nTuple_SBT.py | 6 +-- muonDIS/make_nTuple_Tr.py | 6 +-- shipgen/MuDISGenerator.cxx | 17 +-------- 5 files changed, 43 insertions(+), 64 deletions(-) diff --git a/muonDIS/add_muonresponse.py b/muonDIS/add_muonresponse.py index cf52191a36..4b6330fc2d 100644 --- a/muonDIS/add_muonresponse.py +++ b/muonDIS/add_muonresponse.py @@ -3,7 +3,6 @@ import argparse import logging -import os import ROOT as r @@ -23,17 +22,17 @@ args = parser.parse_args() -def read_file(): +def inspect_file(inputfile): """Inspecting the produced file for successfully added muon veto points.""" - inputFile = r.TFile.Open(args.inputfile, "READ") - input_tree = inputFile.cbmsim + input_file = r.TFile.Open(inputfile, "READ") + input_tree = input_file.cbmsim input_tree.GetEntry(0) muons = False for hit in input_tree.vetoPoint: if hit.GetTrackID() == 0: muons = True - inputFile.Close() + input_file.Close() if muons: print("File is already updated with incoming muon info. Nothing to do.") @@ -45,77 +44,74 @@ def read_file(): return True -def modify_file(): - """Modify the input file with muon info.""" - if args.inputfile[0:4] == "/eos": - args.inputfile = os.environ["EOSSHIP"] + args.inputfile +def modify_file(inputfile, muonfile): + """Add information from original muon to input simulation file.""" + logging.warning(f"{inputfile} will be opened for modification") - logging.warning(f"{args.inputfile} opened for modification") - - inputFile = r.TFile.Open(args.inputfile, "UPDATE") + input_file = r.TFile.Open(inputfile, "UPDATE") try: - input_tree = inputFile.cbmsim + input_tree = input_file.cbmsim except Exception as e: print(f"Error: {e}") - inputFile.Close() + input_file.Close() exit(1) # Open the external file with additional vetoPoints - muonFile = r.TFile.Open(args.muonfile, "READ") + muon_file = r.TFile.Open(muonfile, "READ") try: - muon_tree = muonFile.DIS + muon_tree = muon_file.DIS except Exception as e: print(f"Error: {e}") - muonFile.Close() + muon_file.Close() exit(1) - inputFile.cd() + input_file.cd() output_tree = input_tree.CloneTree( 0 ) # Clone the structure of the existing tree, but do not copy the entries # Access the vetoPoint branch in the original and external trees - originalVetoPoints = r.TClonesArray("vetoPoint") - muonVetoPoints = r.TClonesArray("vetoPoint") - combinedVetoPoints = r.TClonesArray("vetoPoint") + original_vetoPoint = r.TClonesArray("vetoPoint") + muon_vetoPoint = r.TClonesArray("vetoPoint") + combined_vetoPoint = r.TClonesArray("vetoPoint") - input_tree.SetBranchAddress("vetoPoint", originalVetoPoints) - muon_tree.SetBranchAddress("muon_vetoPoints", muonVetoPoints) - output_tree.SetBranchAddress("vetoPoint", combinedVetoPoints) + input_tree.SetBranchAddress("vetoPoint", original_vetoPoint) + muon_tree.SetBranchAddress("muon_vetoPoints", muon_vetoPoint) + output_tree.SetBranchAddress("vetoPoint", combined_vetoPoint) for input_event, muon_event in zip(input_tree, muon_tree): interaction_point = r.TVector3() input_event.MCTrack[0].GetStartVertex(interaction_point) - combinedVetoPoints.Clear() + combined_vetoPoint.Clear() index = 0 - for hit in originalVetoPoints: - if combinedVetoPoints.GetSize() == index: - combinedVetoPoints.Expand(index + 1) - combinedVetoPoints[index] = hit + for hit in original_vetoPoint: + if combined_vetoPoint.GetSize() == index: + combined_vetoPoint.Expand(index + 1) + combined_vetoPoint[index] = hit # pending fix to support ROOT 6.32+ index += 1 - for hit in muonVetoPoints: + for hit in muon_vetoPoint: if hit.GetZ() < interaction_point.Z(): - if combinedVetoPoints.GetSize() == index: - combinedVetoPoints.Expand(index + 1) - combinedVetoPoints[index] = hit + if combined_vetoPoint.GetSize() == index: + combined_vetoPoint.Expand(index + 1) + combined_vetoPoint[index] = hit # pending fix to support ROOT 6.32+ index += 1 output_tree.Fill() # Write the updated tree back to the same file - inputFile.cd() + input_file.cd() output_tree.Write("cbmsim", r.TObject.kOverwrite) - inputFile.Close() - muonFile.Close() + input_file.Close() + muon_file.Close() - print(f"File updated with incoming muon info as {args.inputfile}") + print(f"File updated with incoming muon info as {inputfile}") -add_muons = read_file() +add_muons = inspect_file(args.inputfile) if add_muons: - modify_file() + modify_file(args.inputfile, args.muonfile) diff --git a/muonDIS/makeMuonDIS.py b/muonDIS/makeMuonDIS.py index f659678417..64a27b01d0 100644 --- a/muonDIS/makeMuonDIS.py +++ b/muonDIS/makeMuonDIS.py @@ -1,5 +1,5 @@ #!/usr/bin/env python -"""Script to generate DIS events for muons in Pythia6, and sve to a ROOT file (along with the original muon's soft interactions).""" +"""Script to generate DIS events for muons in Pythia6, and save them to a ROOT file (along with the original muon's soft interactions).""" import argparse import logging @@ -56,7 +56,7 @@ def makeMuonDIS(): try: muon_tree = muonFile.MuonAndSoftInteractions except Exception as e: - print(f"Error: {e}") + logging.error(e) muonFile.Close() exit(1) diff --git a/muonDIS/make_nTuple_SBT.py b/muonDIS/make_nTuple_SBT.py index 416b08b303..823fdd4f7a 100644 --- a/muonDIS/make_nTuple_SBT.py +++ b/muonDIS/make_nTuple_SBT.py @@ -10,10 +10,8 @@ # Argument parser setup parser = argparse.ArgumentParser(description=__doc__) -parser.add_argument("-test", dest="testing_code", help="Run Test", action="store_true") -parser.add_argument( - "-p", dest="path", help="path to muon background files", required=False -) +parser.add_argument("--test", dest="testing_code", help="Run Test", action="store_true") +parser.add_argument("-p", "--path", help="path to muon background files") parser.add_argument( "-o", "--outputfile", diff --git a/muonDIS/make_nTuple_Tr.py b/muonDIS/make_nTuple_Tr.py index 85e41fe126..5a834a4001 100644 --- a/muonDIS/make_nTuple_Tr.py +++ b/muonDIS/make_nTuple_Tr.py @@ -21,10 +21,8 @@ default="muonsProduction_wsoft_Tr.root", help="custom outputfile name", ) -parser.add_argument( - "-p", dest="path", help="path to muon background files", required=False -) -parser.add_argument("-test", dest="testing_code", help="Run Test", action="store_true") +parser.add_argument("-p", "--path", help="path to muon background files") +parser.add_argument("--test", dest="testing_code", help="Run Test", action="store_true") args = parser.parse_args() # Variables diff --git a/shipgen/MuDISGenerator.cxx b/shipgen/MuDISGenerator.cxx index f64dc20167..e51ab18b8d 100644 --- a/shipgen/MuDISGenerator.cxx +++ b/shipgen/MuDISGenerator.cxx @@ -35,15 +35,7 @@ Bool_t MuDISGenerator::Init(const char* fileName, const int firstEvent) dPart = 0; dPartSoft = 0; - if (0 == strncmp("/eos", fileName, 4)) { - TString tmp = gSystem->Getenv("EOSSHIP"); - tmp += fileName; - fInputFile = TFile::Open(tmp, "read"); - LOGF(info, "Open external file on eos: %s", tmp.Data()); - } else { - fInputFile = TFile::Open(fileName, "read"); - } - + fInputFile = TFile::Open(fileName, "read"); fTree = dynamic_cast(fInputFile->Get("DIS")); if (!fTree) { LOG(FATAL) << "Error: Tree 'DIS' not found in file"; @@ -103,10 +95,6 @@ Double_t MuDISGenerator::MeanMaterialBudget(const Double_t* start, const Double_ for (Int_t i = 0; i < 6; i++) bparam[i] = 0; - if (!gGeoManager) { - // AliFatalClass("No TGeo\n"); - return 0.; - } // Double_t length; Double_t dir[3]; @@ -124,8 +112,7 @@ Double_t MuDISGenerator::MeanMaterialBudget(const Double_t* start, const Double_ TGeoNode* currentnode = 0; TGeoNode* startnode = gGeoManager->InitTrack(start, dir); if (!startnode) { - // AliErrorClass(Form("start point out of geometry: x %f, y %f, z %f", - // start[0],start[1],start[2])); + LOG(ERROR) << "Start point out of geometry: x " << start[0] << ", y " << start[1] << ", z " << start[2]; return 0.0; } TGeoMaterial* material = startnode->GetVolume()->GetMedium()->GetMaterial(); From dd237ede7bf04439cc8aa51f4eccb27185726b87 Mon Sep 17 00:00:00 2001 From: anupama-reghunath Date: Wed, 27 Nov 2024 14:42:51 +0100 Subject: [PATCH 09/10] fixup! final fixes --- shipgen/MuDISGenerator.h | 92 ++++++++++++++++++++-------------------- 1 file changed, 45 insertions(+), 47 deletions(-) diff --git a/shipgen/MuDISGenerator.h b/shipgen/MuDISGenerator.h index 8e5f8cf01a..e31c2a7676 100644 --- a/shipgen/MuDISGenerator.h +++ b/shipgen/MuDISGenerator.h @@ -1,58 +1,56 @@ -#ifndef PNDMuGENERATOR_H -#define PNDMuGENERATOR_H 1 +#ifndef SHIPGEN_MUDISGENERATOR_H_ +#define SHIPGEN_MUDISGENERATOR_H_ 1 -#include "TROOT.h" #include "FairGenerator.h" -#include "TTree.h" // for TTree -#include "TF1.h" // for TF1 +#include "FairLogger.h" // for FairLogger, MESSAGE_ORIGIN #include "TClonesArray.h" +#include "TF1.h" // for TF1 +#include "TROOT.h" +#include "TTree.h" // for TTree #include "TVector3.h" -#include "FairLogger.h" // for FairLogger, MESSAGE_ORIGIN #include "vector" class FairPrimaryGenerator; class MuDISGenerator : public FairGenerator { - public: - - /** default constructor **/ - MuDISGenerator(); - - /** destructor **/ - virtual ~MuDISGenerator(); - - /** public method ReadEvent **/ - Bool_t ReadEvent(FairPrimaryGenerator*); - virtual Bool_t Init(const char*, int); //! - virtual Bool_t Init(const char*); //! - Int_t GetNevents(); - - void SetPositions(Double_t z_start, Double_t z_end) - { - startZ = z_start; - endZ = z_end; - } - - private: - Double_t MeanMaterialBudget(const Double_t *start, const Double_t *end, Double_t *mparam); - - - protected: - Double_t startZ,endZ; - TClonesArray* iMuon ; - TClonesArray* dPart ; - TClonesArray* dPartSoft ; - FairLogger* fLogger; //! don't make it persistent, magic ROOT command - TFile* fInputFile; - TTree* fTree; - int fNevents; - int fn; - bool fFirst; - - const Double_t c_light = 29.9792458;//speed of light in cm/ns - const Double_t muon_mass = 0.10565999895334244; //muon mass in GeV - - ClassDef(MuDISGenerator,1); + public: + /** default constructor **/ + MuDISGenerator(); + + /** destructor **/ + virtual ~MuDISGenerator(); + + /** public method ReadEvent **/ + Bool_t ReadEvent(FairPrimaryGenerator*); + virtual Bool_t Init(const char*, int); //! + virtual Bool_t Init(const char*); //! + Int_t GetNevents(); + + void SetPositions(Double_t z_start, Double_t z_end) + { + startZ = z_start; + endZ = z_end; + } + + private: + Double_t MeanMaterialBudget(const Double_t* start, const Double_t* end, Double_t* mparam); + + protected: + Double_t startZ, endZ; + TClonesArray* iMuon; + TClonesArray* dPart; + TClonesArray* dPartSoft; + FairLogger* fLogger; //! don't make it persistent, magic ROOT command + TFile* fInputFile; + TTree* fTree; + int fNevents; + int fn; + bool fFirst; + + const Double_t c_light = 29.9792458; // speed of light in cm/ns + const Double_t muon_mass = 0.10565999895334244; // muon mass in GeV + + ClassDef(MuDISGenerator, 1); }; -#endif /* !PNDMuGENERATOR_H */ +#endif // SHIPGEN_MUDISGENERATOR_H_ From 4e83755d195d4fc3ad0ce1d28e38db6d2790e07f Mon Sep 17 00:00:00 2001 From: anupama-reghunath Date: Wed, 11 Dec 2024 18:29:42 +0100 Subject: [PATCH 10/10] Added DIS cross section directly to the sim file Removed Condor Process variable in makeMuonDIS.py --- muonDIS/makeMuonDIS.py | 116 ++++++++++++++++++++++++++----------- shipgen/MuDISGenerator.cxx | 42 ++++++++++---- 2 files changed, 113 insertions(+), 45 deletions(-) diff --git a/muonDIS/makeMuonDIS.py b/muonDIS/makeMuonDIS.py index 64a27b01d0..2bab679906 100644 --- a/muonDIS/makeMuonDIS.py +++ b/muonDIS/makeMuonDIS.py @@ -11,6 +11,39 @@ logging.basicConfig(level=logging.INFO) PDG = r.TDatabasePDG.Instance() +parser = argparse.ArgumentParser(description=__doc__) +parser.add_argument("-f", "--inputFile", help="Input file to use", required=True) +parser.add_argument( + "-i", + "--firstEvent", + dest="first_mu_event", + help="First event of muon file to use", + required=False, + default=0, + type=int, +) +parser.add_argument( + "-n", + "--nEvents", + dest="n_events", + help="Number of muons to generate DIS for", + required=False, + default=10, + type=int, +) +parser.add_argument( + "-nDISPerMuon", + "--nDIS", + help="Number of DIS per muon to generate", + required=False, + default=1000, + type=int, +) + +args = parser.parse_args() +n_events = args.n_events +first_mu_event = args.first_mu_event + def rotate(px, py, pz, theta, phi): """Rotate the daughter particle momentum to align with respect to the muon's momentum.""" @@ -26,29 +59,32 @@ def rotate(px, py, pz, theta, phi): return rotated_momentum.X(), rotated_momentum.Y(), rotated_momentum.Z() +def update_file(filename, final_xsec): + """Update the DIS cross section of the muon to the converged value from Pythia.""" + file = r.TFile.Open(filename, "update") + + original_tree = file.DIS + updated_tree = original_tree.CloneTree(0) + + iMuon = r.TClonesArray("TVectorD") + updated_tree.Branch("InMuon", iMuon, 32000, -1) + + for i, event in enumerate(original_tree): + in_muon = event.InMuon + mu = in_muon[0] + mu[10] = final_xsec[int(first_mu_event + i / args.nDIS)] + iMuon[0] = mu + updated_tree.Fill() + + file.cd() + updated_tree.Write("DIS", r.TObject.kOverwrite) + file.Close() + logging.info("Muon DIS successfully updated.") + + def makeMuonDIS(): """Generate DIS events.""" - parser = argparse.ArgumentParser(description=__doc__) - parser.add_argument("-f", "--inputFile", help="Input file to use", required=True) - parser.add_argument( - "-nJob", "--nJob", type=int, help="Process ID, gives muon index", required=True - ) - parser.add_argument( - "-nPerJobs", - "--nPerJobs", - type=int, - help="The number of muons per file", - required=True, - ) - parser.add_argument( - "-nDISPerMuon", - "--nDIS", - type=int, - help="Number of DIS per muon to generate", - required=True, - ) - - args = parser.parse_args() + final_xsec = {} logging.info(f"Opening input file: {args.inputFile}") muonFile = r.TFile.Open(args.inputFile, "read") @@ -60,14 +96,12 @@ def makeMuonDIS(): muonFile.Close() exit(1) - nPerJob = args.nPerJobs # number of muons handled by the python script - nStart = args.nPerJobs * args.nJob logging.debug(f"Total entries in the tree: {muon_tree.GetEntries()}") - nEnd = min(muon_tree.GetEntries(), nStart + nPerJob) + last_mu_event = min(muon_tree.GetEntries(), first_mu_event + n_events) - logging.info(f"Creating output file: muonDis_{args.nJob}.root") + logging.info("Creating output file: muonDis.root") - outputFile = r.TFile.Open(f"muonDis_{args.nJob}.root", "recreate") + outputFile = r.TFile.Open("muonDis.root", "recreate") output_tree = r.TTree("DIS", "muon DIS") iMuon = r.TClonesArray("TVectorD") @@ -94,11 +128,15 @@ def makeMuonDIS(): mutype = {-13: "gamma/mu+", 13: "gamma/mu-"} myPythia.SetMSTU(11, 11) - logging.info(f"Processing events from {nStart} to {nEnd-1}...") + logging.info( + f"Processing muon events from {first_mu_event} to {last_mu_event-1}..." + ) nMade = 0 - for k in range(nStart, nEnd): + for k in range(first_mu_event, last_mu_event): + cross_sections = [] + muon_tree.GetEvent(k) imuondata = muon_tree.imuondata @@ -123,12 +161,13 @@ def makeMuonDIS(): logging.debug( f"\n\tMuon index:{k} \n\tPID = {pid}, weight = {w} \n\tpx = {px}, py = {py}, pz = {pz}, E = {E},\n\tx = {x}, y = {y}, z = {z}\n" ) - isProton = 1 xsec = 0 - mu = array("d", [pid, px, py, pz, E, x, y, z, w, isProton, xsec, time_muon]) - muPart = r.TVectorD(12, mu) + mu = array( + "d", [pid, px, py, pz, E, x, y, z, w, isProton, xsec, time_muon, args.nDIS] + ) + muPart = r.TVectorD(13, mu) myPythia.Initialize("FIXT", mutype[pid], "p+", p) # target = "p+" myPythia.Pylist(1) @@ -150,6 +189,7 @@ def makeMuonDIS(): for itrk in range(1, myPythia.GetN() + 1): xsec = myPythia.GetPARI(1) + muPart[10] = xsec did = myPythia.GetK(itrk, 2) dpx, dpy, dpz = rotate( @@ -169,9 +209,11 @@ def makeMuonDIS(): dPartDIS.Expand(nPart + 10) # dPartDIS.ConstructedAt(nPart).Use(part) #to be adapted later dPartDIS[nPart] = part + if itrk == 1: - with open(f"sigmadata_{args.nJob}.txt", "a") as fcross: - fcross.write(f"{xsec}\n") + logging.debug(f"cross section of event is {xsec}") + + cross_sections.append(xsec) dPartSoft.Clear() @@ -211,6 +253,8 @@ def makeMuonDIS(): output_tree.Fill() + final_xsec[k] = xsec + nMade += 1 if nMade % 10 == 0: @@ -220,8 +264,12 @@ def makeMuonDIS(): output_tree.Write() myPythia.SetMSTU(11, 6) logging.info( - f"DIS generated for muons (index {nStart} - {nEnd-1}) , output saved in muonDis_{args.nJob}.root, nDISPerMuon = {args.nDIS}" + f"DIS generated for muons (index {first_mu_event} - {last_mu_event-1}) , output saved in muonDis.root, nDISPerMuon = {args.nDIS}" ) + outputFile.Close() + muonFile.Close() + + update_file("muonDis.root", final_xsec) if __name__ == "__main__": diff --git a/shipgen/MuDISGenerator.cxx b/shipgen/MuDISGenerator.cxx index e51ab18b8d..026be6a1b5 100644 --- a/shipgen/MuDISGenerator.cxx +++ b/shipgen/MuDISGenerator.cxx @@ -245,11 +245,13 @@ Bool_t MuDISGenerator::ReadEvent(FairPrimaryGenerator* cpg) // incoming muon array('d',[pid,px,py,pz,E,x,y,z,w,t]) TVectorD* mu = dynamic_cast(iMuon->AddrAt(0)); LOG(DEBUG) << "muon DIS Generator in muon " << int(mu[0][0]); - Double_t x = mu[0][5] * 100.; // come in m -> cm - Double_t y = mu[0][6] * 100.; // come in m -> cm - Double_t z = mu[0][7] * 100.; // come in m -> cm - Double_t w = mu[0][8]; - Double_t t_muon = mu[0][11]; // in ns + Double_t x = mu[0][5] * 100.; // come in m -> cm + Double_t y = mu[0][6] * 100.; // come in m -> cm + Double_t z = mu[0][7] * 100.; // come in m -> cm + Double_t w = mu[0][8]; // weight of the original muon ( normalised to a spill) + Double_t cross_sec = mu[0][10]; // in mbarns + Double_t t_muon = mu[0][11]; // in ns + Double_t DIS_multiplicity = 1 / mu[0][12]; // 1/nDIS // calculate start/end positions along this muon, and amout of material in between @@ -261,9 +263,11 @@ Bool_t MuDISGenerator::ReadEvent(FairPrimaryGenerator* cpg) end[1] = y - (z - end[2]) * tymu; LOG(DEBUG) << "MuDIS: mu xyz position " << x << ", " << y << ", " << z; LOG(DEBUG) << "MuDIS: mu pxyz position " << mu[0][1] << ", " << mu[0][2] << ", " << mu[0][3]; - LOG(DEBUG) << "MuDIS: mu weight " << w; + LOG(DEBUG) << "MuDIS: mu weight*DISmultiplicity " << w; + LOG(DEBUG) << "MuDIS: mu DIS cross section " << cross_sec; LOG(DEBUG) << "MuDIS: start position " << start[0] << ", " << start[1] << ", " << start[2]; LOG(DEBUG) << "MuDIS: end position " << end[0] << ", " << end[1] << ", " << end[2]; + Double_t bparam; Double_t mparam[8]; bparam = MeanMaterialBudget(start, end, mparam); @@ -326,12 +330,12 @@ Bool_t MuDISGenerator::ReadEvent(FairPrimaryGenerator* cpg) -1, false, // tracking disabled mu[0][4], - t_DIS, // shift time of the incoming muon track wrt t_muon from the input file. - w); // weight associated with a spill. + t_DIS, // shift time of the incoming muon track wrt t_muon from the input file. + w * DIS_multiplicity); // muon weight associated with a spill* DISmultiplicity // outgoing DIS particles, [did,dpx,dpy,dpz,E], put density along trajectory as weight, g/cm^2 - w = mparam[0] * mparam[4]; // modify weight, by multiplying with average densiy * track length + w = mparam[0] * mparam[4]; // modify weight, by multiplying with average density * track length for (auto&& particle : *dPart) { TVectorD* Part = dynamic_cast(particle); @@ -340,8 +344,24 @@ Bool_t MuDISGenerator::ReadEvent(FairPrimaryGenerator* cpg) << (*Part)[4]; LOG(DEBUG) << "muon DIS Generator out part pos " << xmu << " " << ymu << "" << zmu; LOG(DEBUG) << "muon DIS Generator out part w " << w; - cpg->AddTrack( - int((*Part)[0]), (*Part)[1], (*Part)[2], (*Part)[3], xmu, ymu, zmu, 0, true, (*Part)[4], t_DIS, w); + + if (int(mu[0][0]) == int((*Part)[0])) { + cpg->AddTrack(int((*Part)[0]), + (*Part)[1], + (*Part)[2], + (*Part)[3], + xmu, + ymu, + zmu, + 0, + true, + (*Part)[4], + t_DIS, + cross_sec); + } else { + cpg->AddTrack( + int((*Part)[0]), (*Part)[1], (*Part)[2], (*Part)[3], xmu, ymu, zmu, 0, true, (*Part)[4], t_DIS, w); + } } // Soft interaction tracks