Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Timing Correction for MuonDIS generation #510

Open
wants to merge 23 commits into
base: master
Choose a base branch
from

Conversation

anupama-reghunath
Copy link
Contributor

@anupama-reghunath anupama-reghunath commented Sep 12, 2024

Attached are the SBT response plots produced for the previous MuonDIS(involving backwards travelling muons), and with the suggested corrections with the new MuonDIS. SBThits are color coded wrt to the time of DIS. New method preserves the incoming muon's time as well as its veto response, and ensures that the event is temporally consistent.

Before:
SBT_event_oldmuDIS
After:
SBT_event_newmuDIS

Copy link
Contributor

@olantwin olantwin left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Very nice improvements. Generally looking good, but some nitpicks.

Once all other issues are resolved, it would be good to run clang-format and ruff over the touched code.

I would suggest fixing issues with --fixup commits, and then at the end we can do a rebase and automatically squash the fixes.
You might want to check out git-absorb to make that process a bit simpler.

macro/run_simScript.py Outdated Show resolved Hide resolved
macro/run_simScript.py Outdated Show resolved Hide resolved
macro/run_simScript.py Outdated Show resolved Hide resolved
macro/run_simScript.py Outdated Show resolved Hide resolved
muonDIS/add_muonresponse.py Outdated Show resolved Hide resolved
shipgen/MuDISGenerator.cxx Outdated Show resolved Hide resolved
shipgen/MuDISGenerator.cxx Outdated Show resolved Hide resolved
shipgen/MuDISGenerator.cxx Outdated Show resolved Hide resolved
shipgen/MuDISGenerator.cxx Outdated Show resolved Hide resolved
shipgen/MuDISGenerator.h Outdated Show resolved Hide resolved
@olantwin
Copy link
Contributor

olantwin commented Oct 3, 2024

_ No description provided. _

Could you add some explanation to your first post, maybe with the before/after plots you sent me?

@olantwin
Copy link
Contributor

olantwin commented Oct 3, 2024

This change would also benefit from a description in the CHANGELOG (probably an entire paragraph!)

@olantwin
Copy link
Contributor

@maksymovchynnikov Could you please have a look and ideally test these changes? It would be good to merge these ASAP.

@olantwin olantwin added this to the later milestone Oct 30, 2024
@olantwin
Copy link
Contributor

@maksymovchynnikov Your monthly reminder that your feedback would be very welcome on this pull request.

muonDIS/make_nTuple_Tr.py Outdated Show resolved Hide resolved
@olantwin
Copy link
Contributor

@anupama-reghunath I also noticed that there are still a lot of open comments.

@olantwin olantwin linked an issue Nov 12, 2024 that may be closed by this pull request
@olantwin olantwin removed the request for review from maksymovchynnikov November 13, 2024 14:43
macro/run_simScript.py Outdated Show resolved Hide resolved
@olantwin
Copy link
Contributor

olantwin commented Dec 3, 2024

For any new code, we should not use TVectorD anymore, I'll make an exception for this PR, as the muon DIS was already using it.

We should consider refactoring to use STL containers.

# 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:
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Reasoning for having an external textfile to save the DIS cross section from Pythia6:

  1. As is, the DIS cross section otherwise is not saved anywhere
  2. Consider generating N >>1 DIS events for the single muon. Since Pythia generates the cross-section value events-by-event, it converges from some meaningless value for the 1st event to the proper value for the Nth event. In the proper implementation, all the events must share the same correct cross-section. As a temporary fix,we achieve this by reading the Nth DIS cross section entry for a given muon to calculate the weight.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

My plan of action to avoid this external file is:

1.Directly write the cross section to the muonDis file ( which is passed as the input for muonDIS production)
2. Set weight of the scattered muon to the cross section in the simulation (atm both the incoming and scattered have the same weight from the spill).

Perhaps also set a minimum value for the number of DIS per muon to allow for convergence?

Please let me know what you think! @olantwin

Copy link
Contributor Author

@anupama-reghunath anupama-reghunath Dec 6, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Heres a sample from 200 DIS events per muon from a test sample. There is a reset at the halfpoint mark when the fixed target changes from proton to neutron.
Cross_sections

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, this is now much clearer, thanks!

I agree, that this information should be directly in the sample. Is there a weight field in the event header? Otherwise, using the weight of the scattered muon seems like an OK workaround (maybe we can preserve the weight of the incoming muon from the target sim).

We do not use the cross section anywhere else, right?

Could you prepare that plot also relative to the final cross section? That allow us to estimate a good default number of events to generate.

Copy link
Contributor Author

@anupama-reghunath anupama-reghunath Dec 11, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree, that this information should be directly in the sample. Is there a weight field in the event header?

I couldnt find one in ShipEventHeader.

Otherwise, using the weight of the scattered muon seems like an OK workaround (maybe we can preserve the weight of the incoming muon from the target sim).

At the moment I save the cross section onto the scattered muon and also refactor the incoming muon's weight to include the DIS multiplicity such that we maintain the event weight's normalisation to a spill. Previously this refactoring was explicitly done in the analysis, but the meta data (nDIS per muon) needed to be explicitly plugged in to calculate weight for any sim file.

We do not use the cross section anywhere else, right?

Only when calculating the weight of the event. ( Done during analysis)

Could you prepare that plot also relative to the final cross section? That allow us to estimate a good default number of events to generate.

For 2000 DIS events, the convergence seems to be inversely correlated with the energy of the muon.
image

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've removed the usage of an external file, and now save the DIS cross section to MCTrack[1] instead of scattered muon(simply out of convenience for easy calculation of weights).

@olantwin
Copy link
Contributor

@anupama-reghunath How shall we proceed with this PR?

@anupama-reghunath
Copy link
Contributor Author

Summary of the latest commits:

  • An overlooked bug with saving the incoming muon's vetopoint information. The issue arises if we have multiple muons in the events and then the vetoPoint information muon hitting the SBT/Tr may not be correctly chosen. This has now been rectified.
  • Adding Carbon onto the PDG Database in run_simScript.py (From Coulomb scattering of LiSc (Linear Alkyl Benzene))
  • Lowering the MCTrack Energy Cut to 10MeV in run_simScript.py , this helps with studying veto efficiencies, but also avoids losing Tracks with vital info ( MCTrack[1].GetWeight() contains the DIS cross section of the incoming muon).

Comment on lines 378 to 380
# add Carbon to PDG
pdg = ROOT.TDatabasePDG.Instance()
pdg.AddParticle("C", "Carbon_Nucleus", 12.0, True, 0, 6.0, "nucleus", 1000060120)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe we would want to add this where we add the other particles to the database:

def configure(darkphoton=None):

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

wasnt sure where to add this, maybe we should then also add line 445: pdg.AddParticle('W','Ion', 1.71350e+02, True, 0., 74, 'XXX', 1000741840)

Comment on lines +522 to +526
if options.mudis:
EnergyCut=10.*u.MeV
else:
EnergyCut=100.*u.MeV

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe better as a ternary statement to avoid possible issues of variable scope.

@olantwin
Copy link
Contributor

Are you adding any new TVectorD's compared to what was used before?

@anupama-reghunath
Copy link
Contributor Author

anupama-reghunath commented Jan 30, 2025

Are you adding any new TVectorD's compared to what was used before?

ah yes, but this is metadata -> ( number of muons within the original MuonBack event that pass the threshold), didn't know how else to save this other than in muonDis.root , incase of recalculation of weights.

@olantwin
Copy link
Contributor

Are you adding any new TVectorD's compared to what was used before?

ah yes, but this is metadata -> ( number of muons within the original MuonBack event that pass the threshold), didn't know how else to save this other than in muonDis.root , incase of recalculation of weights.

It's fine to store it in the file, but please use a std::vector instead, see #562.

@anupama-reghunath
Copy link
Contributor Author

anupama-reghunath commented Jan 30, 2025

Are you adding any new TVectorD's compared to what was used before?

ah yes, but this is metadata -> ( number of muons within the original MuonBack event that pass the threshold), didn't know how else to save this other than in muonDis.root , incase of recalculation of weights.

It's fine to store it in the file, but please use a std::vector instead, see #562.

To be specific, there is no new TVectorD but just an additional element in the existing TVectorD. I did so to avoid making new changes to the rest of the code implementation.

@olantwin
Copy link
Contributor

Are you adding any new TVectorD's compared to what was used before?

ah yes, but this is metadata -> ( number of muons within the original MuonBack event that pass the threshold), didn't know how else to save this other than in muonDis.root , incase of recalculation of weights.

It's fine to store it in the file, but please use a std::vector instead, see #562.

To be specific, there is no new TVectorD but just an additional element in the existing TVectorD. I did so to avoid making new changes to the rest of the code implementation.

Ok, that's fine then. If we're just extending an existing object I can live with it.

@olantwin olantwin modified the milestones: Later, 25.02 Jan 30, 2025
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

muon DIS and inactivateMuonProcesses
2 participants