Generating events with MCatNLO 3.1

Prepare the package

• Download Package31.tar.gz and pdf file data.tar.gz from some where (your know?)

cd $HOME
cd myAnalysisCode
mkdir gensim; cd gensim
mkdir bin; cd bin; tar zxvf Package31.tar.gz; cd ..
mkdir pdf; cd pdf; tar zxvf data.tar.gz

Compilation of the package

• For LHC, overwrite MCatNLO.inputs with this file

cd $HOME/myAnalysisCode/gensim/mcnlo/bin
rm Linux -rf
. MCatNLO.inputs     # it will define needed env and some bash programs
compile NLO              # compile the package.  you can "compile MC"  too, 

Integration

• Download this file or take ttbarNLOinput in directory "Linux" and modify as you need, this file will overwrite the parameters defined in MCatNLO.inputs

cd Linux
dothelinks 
./ttbarNLO_EXE_THISLIB < Integration.Inputfile

• The integral will take a time depending on how precise you want it. it is needed to be done only once, the next time you can set the last line of Integration.Inputfile to be "0 0" to avoid integral once again.
• 3 files are the results of the integration, save them in a specific directory

mkdir $HOME/myAnalysisCode/gensim/mcnlo/cfg1
cp ttbarNLO_EXE_THISLIB  $HOME/myAnalysisCode/gensim/mcnlo/cfg1
cp ttbar_a_bs.data  ttbar_b_bs.data  ttbar.integrals  $HOME/myAnalysisCode/gensim/mcnlo/cfg1
cp $HOME/myAnalysisCode/gensim/mcnlo/pdf/cteq6m.tbl  $HOME/myAnalysisCode/gensim/mcnlo/cfg1/cteq6m
cp Integration.Inputfile  $HOME/myAnalysisCode/gensim/mcnlo/cfg1

• You can also do as following to unite compilation and integration into one step. Integration.Inputfile will be substituted by ttbarNLOinput produced according to the definition in MCatNLO.inputs

. MCatNLO.inputs 
runNLO

Event Generation

• In fact, the integration process has generated some events, too, depending on the number of events you defined.
• To generate more events:
  • Once you have copied several needed file to cfg1 (above), you can just go to that directory, reset number of events and set last line of Integration.Inputfile to be "0 0" to avoid integral process.
  • Run command ./ttbarNLO_EXE_THISLIB < Integration.Inputfile
  • it's very quick, and you get ttbar.events ttbar_b.events ttbar_a.events as the generating results.

top quark decay and showering

• with herwig http://hepwww.rl.ac.uk/theory/seymour/herwig/hw65_manual.htm

method1

• In any work directory, compose athena_setup.sh and source athena_setup.sh
  • Make DC3.005200.T1_McAtNlo_Jimmy.py according to the original one copied from athena standard.
  • get_files PDGTABLE.MeV
  • compose file inparmMcAtNlo.dat
  • copy ttbar.events produced in last step here and replace the full line CTEQ   56 .... with HWLHAPDF  10000
  • run athena DC3.005200.T1_McAtNlo_Jimmy.py

method2

• Still you need to have the XXX.events and XXX.dat files firstly same as above
• Look at DC3.005200.T1_McAtNlo_Jimmy.py , check the value for evgenConfig.inputfilebase = 'mcatnlo31.005200.ttbar'
• change the XXX.events and XXX.dat file to be consistent with mcatnlo31.005200.ttbar.events and mcatnlo31.005200.ttbar.dat, more over, don't forget to change the .events name in mcatnlo31.005200.ttbar.dat
• tar zcvf mcatnlo31.005200.ttbar.tar.gz mcatnlo31.005200.ttbar.dat mcatnlo31.005200.ttbar.events
• csc_evgen_trf.py  runNumber=5200 firstEvent=1 maxEvents=5000 randomSeed=234324 jobConfig./DC3.005200.T1_McAtNlo_Jimmy.py outputEvgenFile=event.pool.root histogramFile=NONE ntupleFile=Ntuple.root inputGeneratorFile=./mcatnlo31.005200.ttbar.tar.gz=
  • have to set run number to be 5200, XXX.events should have more events contained.
  • event number have to be 5000 or more
  • official DC3.005200.T1_McAtNlo_Jimmy.py will be searched and used, if there is no ./ before it. This file in mehtod 1 has been modified.

Atlfast process

• compose file atlfast.py
• if you want to change the input showered file name get_files Atlfast_ReadPOOLFile.py and change in it
• If you want to change the output AOD file name get_files Atlfast_MakeAOD.py and change in it
• athena atlfast.py

• or you can use csc_atlfast_trf.py evgen.pool.root AOD.pool.root NONE 5000 0 to finish this step
• or you can use csc_genAtlfastTwoStep_trf.py to finish the step and last step together.

Testing with 100 events

• Do integral with low precision (1,1), and produced 100 events : ttbar.events
• Showering these 100 events, to get file DC3_005200_McAtNlo_Jimmy_02.pool.root, the log file is mnnlo_herwig_showering.log
  • 100 events, but in the showering, event 4 29 32 61 66 are killed, finally 95 events are showered (0-94) ( There is one subtlety. when HERWIG reconstructs the initial state parton shower from the chosen (x′,Q′2) pair, it might reject it due to mismatches between the Q2 ordering in the evolution equations of the employed parton distribution functions and the angular ordering in HERWIG’s evolution. This is taken care of by rejecting the event.
  • It still try showering events (eventnumber 95) but because no events left, it produce a bad event and cause the problem below.
  • so now the real problem is in inparmMcAtNlo.dat, we set Nevent in this file to tell how many events will be processed, but it doesn't work. you can put any number, but athena process all events. We can set theApp.EvtMax in job option, but we don't know how many events left each time.
• Use atlfast to produce the AOD.pool.root, and use AOD analysis code to readout the truth decay chain truthinfo_my100events.txt
  • problem: the program crash in eventnumber 95, only eventnumber=0 to 94 works well

There is one subtlety. when HERWIG reconstructs the initial state parton shower from the chosen (x′,Q′2) pair, it might reject it due to mismatches between the Q2 ordering in the evolution equations of the employed parton distribution functions and the angular ordering in HERWIG’s evolution. This is taken care of by rejecting the event.

Generating events with MCatNLO 3.31

Prepare the package 3.31

• Download Package33.tar.gz and pdf file data.tar.gz from http://www.hep.phy.cam.ac.uk/theory/webber/MCatNLO/
• unpack v3.31 patch to the same directory where v3.3 untarred.

cd $HOME/myAnalysisCode
mkdir -p gensim/mcnlo3.31; cd gensim/mcnlo3.31
mkdir bin; cd bin; tar zxvf Package33.tar.gz; tar Addendum331.tar.gz; cd ..
mkdir pdf; cd pdf; tar zxvf data.tar.gz

Compilation integration and production

• Open file MCatNLO.inputs and comment the command runMCatNLO. *the file will only define funcitons. it also defined many variables but many are overwritten by the parameter file later.

• cd $HOME/myAnalysisCode/gensim/mcnlo3.31/bin
• rm Linux -rf
• . MCatNLO.inputs
• compile NLO # compile the package ==> compile log file
• cd Linux
• ./ttbarNLO_EXE_THISLIB < czhu.Integration.Inputfile # ( do or not do integration according to the last line of Integration.Inputfile) and produce events

• a few events can be produced.

More Event Generation with 3.31

• To generate more events, Run command ./ttbarNLO_EXE_THISLIB < Integration.Inputfile again, with more event numbers set
  • Problem:
    • in subroutine getspincorr(jproc0), an almost dead loop emerged, which make <100 events production comsume 24 hours.
    • line 4983 in file mcatnlo_qqmain.f (v3.31) make the loop again and again.
  • Solution:
    • The reason is to set top width to 0. Now set width as 0.000000001 to bypass this prolbem.

Saving configuration for specific generating

• • copy ttbarNLO_EXE_THISLIB  ttbar_a_bs.data  ttbar_b_bs.data  ttbar.integrals   needed pdf file cteq6m.tbl  cteq6m Integration.Inputfile to a standalone directory
  • Once you have copied the needed file to a standalone directory, you can just go to that directory, reset number of events and set last line of Integration.Inputfile to be "0 0" to avoid integral process.

events file to evgen file

• Use Atlas Release 14.0.10 and set the env.
  • Look at CSC.005200.T1_McAtNlo_Jimmy.py , check the value for evgenConfig.inputfilebase = 'mcatnlo31.005200.ttbar'
• Rename XXX.events to be mcatnlo31.005200.ttbar.events, open the file and change cteq  56 to LHAPDF    10000 (4 spaces between LHAPDF and 10000) or HWLHAPDF  10000 (2 spaces between HWLHAPDF and 10000)
• Use this file as mcatnlo31.005200.ttbar.dat to be consistent with athena 14.0.10
• tar zcvf mcatnlo31.005200.ttbar.tar.gz mcatnlo31.005200.ttbar.dat mcatnlo31.005200.ttbar.events
• csc_evgen_trf.py  5200 1 5000 234324 CSC.005200.T1_McAtNlo_Jimmy.py evgen.pool.root NONE NONE ./mcatnlo31.005200.ttbar.tar.gz
  • running log
  • have to set run number to be 5200, and maxevent >5000.
  • event in mcatnlo31.005200.ttbar.events have to > maxevent
  • official CSC.005200.T1_McAtNlo_Jimmy.py will be searched and used, if there is no ./ before it.

evgen passes ATLFAST simulation

• csc_atlfast_trf.py inputEvgenFile='evgen.pool.root' outputAODFile='aod.pool.root' ntupleFile=NONE maxEvents=5000 skipEvents=0
  • error message
  • change to 14.1.0.2 on lxplus, Still problem: error message
    • use csc_evgen08_trf.py + CSC.005200.T1_McAtNlo_Jimmy.py, events can't pass energy conservation check "WARNING ENERGY BALANCE FAILED : E-difference = 4e+06 MeV" --> error message
    • use csc_evgen08_trf.py + MC8.105200.T1_McAtNlo_Jimmy.py in 14.1.0.3, same error as above

Use athena 14.1.0.2 on lxplus.cern.ch succeed with the 2 steps above

• Set up 14.1.0.2 environment: source $HOME/setup_athena/cmthome/setup.sh -tag=14.1.0.2,AtlasProduction,releases,runtime,slc4,32(?gcc34?)

• Run evgen: csc_evgen_trf.py runNumber=5200 firstEvent=0 maxEvents=5000 randomSeed=12345 jobConfig=my.evgen.py outputEvgenFile=EVNT.root histogramFile=NONE ntupleFile=NONE inputGeneratorFile=mcatnlo331.005200.ttbar.tar.gz
  • my.evgen.py = CSC.005200.T1_McAtNlo_Jimmy.py except change input file to mcatnlo331.005200.ttbar
  • Old TTbarLeptonFilter rejected leptons which came from W decayed by MC@NLO.
    • New TTbarWToLeptonFilter in GeneratorFilters-00-00-97 fixes the problem.

• Run atlfast: csc_atlfast_trf.py inputEvgenFile=EVNT.root outputAODFile=AOD._00001.root ntupleFile=NTUP._00001.root maxEvents=5000 skipEvents=0

• used .dat and .events files are:

Show Hide

'mcatnlo31.005200.ttbar.events'        ! event file
 500                        ! number of events
 1                        ! 0->Herwig PDFs, 1 otherwise
'P'  'P'               ! hadron types
 7000.00 7000.00               ! beam momenta
 -1706                          ! -1705/1706=bb/tt
 'LHAPDF'                      ! PDF group (1)
 10000                       ! PDF id number (1).
 'LHAPDF'                      ! PDF group (2).
 10000                       ! PDF id number (2).
 -1                     ! Lambda_5, < 0 for default
 175 80.4                 ! M_t, M_W
 0.32 0.32 0.5 1.5 5 0.75 ! quark and gluon masses


 0.1400E+05 0.1000E+01 0.1000E+01 0.1000E+01 0.1000E+01 --> CM energy, muR/mu0[NLO], muF/mu0[NLO], muR/mu0[MC], muF/mu0[MC]
   1706 --> 1705/1706=b/t
 0.1750E+03 0.1000E-08 --> M_top, Gamma_top
 0.8040E+02 0.2124E+01 --> M_W, Gamma_W
  4  5 --> IL1, IL2 (0..7)
 0.3200E+00 0.3200E+00 0.5000E+00 0.1500E+01 0.5000E+01 0.7500E+00 --> quark and gluon masses
 P    P    --> colliding particles
 LHAPDF    10000
 0.2260E+00 MS --> Lambda_5, scheme
 P,M --> Format of v3.1 and higher
     8000
 407  4  9
  21  21  21 -13  14   5   1  -2  -5
 0.40205383E+17
 0.00000000E+00 0.00000000E+00 0.92184023E+02 0.75000000E+00 0.00000000E+00 0.00000000E+00 -.11693000E+04 0.75000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 -.88104968E+00 0.14612050E+02 -.23479250E+02 0.10565836E+00 0.82095989E+02 0.30725785E+02 -.82662758E+03 0.00000000E+00 0.13810531E+02 -.61199529E+02 -.18900338E+03 0.50000000E+01 -.32023458E+02 0.77235047E+02 0.26255570E+01 0.32000000E+00 -.40444789E+02 -.11728112E+02 0.11955713E+02 0.32000000E+00 -.22557223E+02 -.49645242E+02 -.52587072E+02 0.50000000E+01

• Problem with produced AOD file. AOD._00001.root is produced with only EventInfo_p2_McEventInfo and an unwanted AOD.pool.root is produced with full infomation but can't pass my AOD processing code (error message)
  • first understanding of this problem:
    • Use valgrind to monitor the above processing shows that program Jump to the invalid address stated on the next line at 0x4E455645: ???  by 0x2F4005E8: SG::DataBucket<DataVector<IParticle, DataVector_detail::NoBase> >::~DataBucket(). DataBucket is a tool which takes ownership of a data object to be put in the DataStore and allows to access the transient object address as a void*, this means that when the program try to point to the data in memory, it failed, which may be caused by the uncorrect structure of the data.
  • possible sollution: find other way to produce the correct data.

Generating Events With AcerMC

Install the package

• Go to http://borut.home.cern.ch/borut/, download AcerMC 3.4
• bunzip2 AcerMC-3.4.tar.bz2
• Go to $HOME/myAnalysisCode/gensim/
  • tar xvf AcerMC-3.4.tar
• got directory AcerMC-3.4
  • make demo_py
  • make run_py
  • directory prod will be changed, all files concerning running is in it.

Generating events

• just go to $HOME/myAnalysisCode/gensim/AcerMC-3.4/prod
  • modify run.card and acermc.card according to your need, especially set ACER 2 in run.card, which will ouput AcerMC_pXXX_rYYYYYY.events file, where XXX is physics process ID and YYYYYY is random number.
  • modify prefix=./lhapdf in lhapdf-config accoding to where lhapdf is
  • use b shell and put the "prod" directory into PATH then run ../demo_py.e
  • *.out, AcerMC_pXXX_rYYYYYY.events and AcerMC_pXXX_rYYYYYY.inparm are the outputs

Simulate with Atlfast and making AOD

method 1

• Reference: RunningAtlfast
• Use this file DC3_generation_skeleton.py. you can change some if you understand what it is
• copy AcerMC_pXXX_rYYYYYY.inparm here, rename or make a symbol link to inparmAcerMC.dat
• copy AcerMC_pXXX_rYYYYYY.events here
• Run athena DC3_generation_skeleton.py
• AOD file is generated
• Attention: Setting theApp.EvtMax in the last line of MAIN job option files is important!! because, this may have been set many times in the included job option files, if you set it before these job option files, your setting will be overwritten.

method 2

• Producing evgen.pool.root ( *.events file should at least containes >5000 events, otherwise csc_evgen_trf.py will complaint you input so few events to trouble it )

      mv AcerMC_pXXX_rYYYYYY.inparm AcerMC_pXXX_rYYYYYY.dat
      tar zcvf AcerMC_pXXX_rYYYYYY.tar.gz  AcerMC_pXXX_rYYYYYY.dat AcerMC_pXXX_rYYYYYY.events
      csc_evgen_trf.py  5205 1 5000 234234 ./DC3.005205.AcerMCttbar.py evgen.pool.root 'NONE' 'NONE'  AcerMC_pXXX_rYYYYYY.tar.gz

• Each option's mean, please reference here: https://twiki.cern.ch/twiki/bin/view/Main/JobTr
• Producing AOD file from evgen (this produced AOD is not checked)
  • csc_atlfast_trf.py evgen.pool.root AOD.pool.root NONE 5000 0
• Merge the 2 steps together with csc_genAtlfastTwoStep_trf.py use -h to see help.

Hardronization interface in athena Acermc_i (6.5.0 and later)

Script for automatic production

• auto_acermc.sh

Check the data

truth information of 1 event: one particle is from other 22 particles in the decay chain.

• • ReadTruthParticles   INFO |  |  |  |-child: 0 PDGID: 92 nChildren: 61 status: 2 bc: 267 bc Parents: 22 [ 110 111 112 113 114 115 116 117 118 119 121 122 123 124 125 126 127 107 108 109 120 128 ] pxpypze: 13141.4 : 300.892 : -1.90694e+06 : 2.01103e+06
  • in 500 ttbar AcerMC events, totally 4600 particles (4595 * pdgid92 + 5 *pdgid91) are from more than 2 mother particles. they are from 3, 4, 5 ... 60 particles as following plot. This process is string fragmentation model in Pythia, 92 represents a fragmentation process.:
    

The truth infomation problem

• The energy and momtum of the AcerMC generated hard-process-particles (see "show" 1st column) are different with the corresponding particles from the AOD truth container (2nd column). the 3rd column is the ratio of the first 2 numbers. see:
  • The truth information in different stage of simulation: [[https://twiki.cern.ch/twiki/bin/view/Atlas/McEventCollection][McEventCollection]

Show Hide

 231.21   315.19   1.36
317.51   235.12   0.74
293.03   194.29   0.66
426.34   268.81   0.63
201.7   2118.93   10.51
1088.76   1839.81   1.69
1708.88   1127.18   0.66
1013.15   521.03   0.51
1586.21   3152.51   1.99
220.89   247.75   1.12
485.75   628.46   1.29
379.03   244.48   0.65
571.13   303.29   0.53
1125.99   1591.58   1.41
276.42   297.59   1.08
304.05   284.05   0.93
463.59   326.88   0.71
648.42   803.64   1.24
912.91   517.05   0.57
284.14   204.99   0.72
699.8   966.45   1.38
239.26   217.17   0.91
295.87   309.83   1.05
441.41   319.51   0.72
803.38   316.97   0.39
465.33   476.1   1.02
1401.41   302.36   0.22
601.62   563.16   0.94
759.4   329.32   0.43
456.85   817.9   1.79
247.89   262.54   1.06
437.71   313.47   0.72
525.18   1322.24   2.52
272.23   208.83   0.77
1064.21   427.07   0.4
468.71   603.17   1.29
832.05   1751.4   2.1
186.08   194.78   1.05
557.47   1310.54   2.35
343.63   312.37   0.91
2130.59   772.31   0.36
257.01   412.28   1.6
382.24   262.1   0.69
529.56   455.09   0.86
206.99   350.4   1.69
725.51   599.87   0.83
225.34   202.34   0.9
221.66   235.8   1.06
778.06   968.75   1.25
279.24   395.6   1.42

• Reason: AcerMC produced hard process is input to Pthia.
  • Pythia does ISR, which re-organize the energy and momentum of the hard-process-particles. If ISR is closed, the energy of the hard-process-particles will be similar with the input ones. but the difference in momentum is still larger.
  • The multiple scattering: Pythia does multiple scattering. If it is closed, the momentum of the hard-process-particles get similar with the input ones
  • There is still little difference (~1GeV), don't know the reason.

incorrect spin correlation parameter A in truth

• the spin correlation factor truth A computed from AOD truth information are not 0.42(SM expect), but 0.49.
• Test Summanry:
  • Using AcerMC3.1 process 11(90%) + 12(10%), with or without ISR and multiple scattering, got the same results, they are not the reason
  • Using AcerMC3.4 and process 23 which mixes the process 11 and 12, correct value is retrieved as below plot.
  • 
  • Look at the pythia output file (work after AcerMC3.4), the particles after radiation and multiple scattering ... are still much close to the hard process, except px and py are much different. AcerMC 3.1 has a bug on this issue?
• Is it the reason (different gg and qq production ratio)
  • gg->tt 100% : A=0.605033, 158508 events,
  • gg->tt 90% and qq->tt 10%: A=0.491676 similar number of events
  • gg->tt 86% and qq->tt 14% ?? in AcerMC3.4, A=0.42

A PDGID list

• PDG.py.txt

ZhuPrivateDataValidation

Simulation with FastCaloSim

• ZhuFastSim2

• log4.txt: log4.txt