Simulation using SIMC
---Working on simc version taken from Github---
---This is Barak Schmookler's activity Log---
Sep 28, 2015:
1) Created 'simc_gfortran' fork on GitHub
2) Cloned to GMp workdisk
Sep 29, 2015:
1) Changed Makefile for 64bit library compatibility
2) Created "save_sep15" directory in "infiles" directory. Moved initial input files here. Began to create GMp input files in the "infiles" directory. Specifically I created a "point hydrogen" inclusive elastic input file and a 15cm long hydrogen inclusive elastic input file.
3) Copied correct forward and reconstruction maps (COSY) to hrsr and hrsl directories. Also copied .fox file used to generate said matrices and .fox file that can be used to plot rays (within COSY).
Sep 30, 2015:
1) Added input file option (col_flag) which allows one to choose an open, collimator, or sieve configuration. Previously this had to be done in the spectrometer code itself, and then recompiled. I currently only implemented the option for the HRSs (not for Hall C spectrometers).
Nov 4, 2015:
1) Added line to "enerloss_new" subroutine to print out the mode ('typeflag' in subroutine) when Eloss > Ekin. Also made to print values of Eloss and Ekin.
Nov 10, 2015:
1) Adjusted "run_simc" script to dump output to 'outfiles' directory
Nov 21, 2015:
1) Implemented GMp target geometry. The GMp target will be used if "target%can" is set to 3.
1) Created version of SIMC ("simc_sos") where the normal Q1 on both hrs spectrometers is replaced by the SOS quad.
2) Created single-arm versions called "phase_space_hrs" and "phase_space_sos". This version has an option to dump all events into the ntuple. Multiple scattering and energy loss is not implemented.
3) Two external sieves have been implemented. One is called "GMp Sieve", which is the one that was used on the Right HRS during Spring 2016.
March 30, 2016:
1) Implemented number generation using a random state file. This is done so that random seed is different for different runs. In the standard SIMC, this option already existed. In the single arm version, I had to write some code.
2) Set extra_hydrogen input files to use flat raster pattern. In single-arm versions, only flat raster option is available. I've also dumped the raster values into the ntuple for all versions (already present for y direction.)
March 31, 2016:
1) Implemented material for multiple scattering in single-arm versions. I used the RHRS Spring 16 geometry. I also set a "beer can" target as the only option.
April 11, 2016:
1) For the scattered electron, made it so the forward-going ionization energy loss is protected with the "using_Eloss" flag. See lines 1611-1620 in simc.f.
2) Corrected line 1631 in simc.f to use the angle of the electron arm. Previously, it incorrectly used the angle of the proton arm.
3) In the single-arm version, corrected a formatting error when writing event loss information to the .hist file in the outfiles subdirectory.
4) In the single-arm version, made a rough calculation of reconstructed z-tar. The vertex resolution is now calculated (and printed) correctly.
5) In the "extreme_trip_thru_target" subroutine in target.f, the diameter of the "beer can" target was incorrect. See lines 387-403 and lines 539-553.
April 12, 2016:
1) Implemented "dump all events" option and event-loss flag/location for all versions. If one chooses to put all generated events into the Ntuple, a branch named "success" ("ok_spec") tells whether the event reached the hut and was reconstructed. The cross-sections/weights are calculated for every generated event. The information in the .hist file is not affected, however. The variables "x_stop", "y_stop", and "stopwhen" provide info on the events that are not successful. Reconstructed quantities will only be updated on a success; the event-loss variables ("x_stop","y_stop","stopwhen") will only be updated on a failure; so it's important to cut on the "success" branch. If the events are run through both spectrometers, the event is first run through the hadron spectrometer. So, only hadron arm event-loss information will be written if the event fails there. For GMp, we only run the events through the electron arm, but someone may want to change this at some point. (N.B. All of this is only implemented for HRS at the point).
April 25, 2016:
1) Dumped real and reconstructed ionization energy loss information for the incoming electron and scattered particles.
2) Corrected sign error in transformation from hall coordinate system to spectrometer coordinate system. See line 1421 for hadron and line 1646 for electron.
3) Dumped generated and reconstructed z-vertex into ntuple.
4) Made reconstructed energy loss calculations z-vertex dependent. When considering the incoming electron energy, I use the average reconstructed vertex of the electron and hadron. If the events are only run through one arm, I obviously only use that arm's reconstructed vertex.
May 4, 2016:
1) In single-arm version added option to read in generated events from a text file. Hat tip to Erez Cohen on this one.
May 25, 2016:
1) Put correct geometry for Left HRS and Right HRS for Spring 2016 run into a directory called "simc_spring16". Also, I included the correct material before spectrometer.
2) Created a single-arm version for Spring 2016 RHRS, called "phase_space_spring16".
June 23, 2016:
1) Put in vdc offsets from Analyzer database for Left HRS and Right HRS into mc_hrsr(l).f; as well as into detector hut subroutine.
2) In "phase_space_spring16", put correct material before spectrometer (assuming RHRS). I assume the target is shaped as a "beer can".
August 30, 2016:
1) When using benhar spectral function, previously the limits on missing momentum were not being set. So, a small change was made in dbase.f (see lines 509-510). Probably a better fix than this could be developed.
September 2, 2016:
1) Updated the Right HRS detector stack with the corrected material and detector positions. I noticed that the previous version of SIMC had the focal plane in a slightly incorrect position. So when I put the correct position in SIMC, I had to also use this correct position in COSY.
October 4, 2016:
1) For both HRS spectrometers, rotated the vacuum exit windows so that they lie in the horrizontal plane. This aperture does not seem to define the acceptance at all. So it's not a very important change.
October 23, 2016:
1) Updated the Left HRS detector stack with the corrected material and detector positions. Same comment on focal plane as for Right HRS.
2) Made the raster correction account for the average beam position.
October 28, 2016:
1) Created a simc_fall16" version, with updated hrs, detector hut, target info.
November 6, 2016:
1) Created a single-arm version for Spring 2016 RHRS, called "phase_space_fall16".
January 26, 2017:
1) Made option to run events through Left HRS in all 4 phase-space versions.
February 10, 2017:
1) Added ionization energy loss option for outgoing particle in phase-space-hrs version. Also added some additional variables to output file. Will do the same for other 3 versions.
To do (may be extended if necessary...):
1) Merge all (or most) of the different versions (probably won't happen for awhile).
2) Make ionization energy loss correction consistent with what is done in Analyzer.
3) Add some information about the changes made to the .hist file in the outfiles directory. Can also do this for the single-arm version.
1) When calculating "geometrical" jacobian, comments are wrong but calculation is right. (i.e. cos(dtheta) should be inverted and dtheta is not theta - theta0.)
2) When a spectrometer is not used, the y_tar variable is set to zero. The other 3 spectrometer variables (delta, x', y') are set to approximately the generated values. Only these 3 variables need to be within the limits to add to the written-out resolution and accumulated weight (as well as to calculate physics variables like p_miss), so this is acceptable. But I guess at some point someone may want to do something similar for the y_tar variable to aid check the reconstruction of the reaction vertex...and maybe for other reasons.
3) In the single arm version, the comments say that +xbeam is right (looking downstream), but it is actually left.
4) I was initially a little confused about the following... For extended multiple scatters (in the detector hut), the "project" subroutine is called as well as the "musc_ext" routine. This is because the "musc_ext" routine only accounts for the effect of multiple scattering; it does not do the drift itself.
Some documents concerning above changes:
1) Schematics of HRS spectrometer:
2) Target Material, etc. for GMp experiment:
3) HRS sieve designs:
4) Table of locations (relative the the VDC) of detectors in HRS detector stack:
5) Information on HRS detector materials, thickness, etc.:
6) Calculation of reaction z position:
7) Information on Ionization energy loss corrections: