Generate, digitize and reconstruct two signal events with decays

Bd->J/Psi(mumu)Kshort

This tutorial shows you how to simulate, digitize and reconstruct signal events so that they can be analyzed in DaVinci. The versions of the code quoted are available in December 2011 and replace the version numbers with more recent ones if required. This is setup to generate MC11a events which are the latest version at the time of writing but will probably be depreciated by the time you read this.

Using Gauss to simulate events

The slides that explain what the program Gauss does are here.

Setup a version of Gauss from the default installation at your site (or CERN)

Log in to lxplus or if it is installed locally any computer with access to the LHCb software.

cd SomeDirectory
SetupProject Gauss v41r0
cp $GAUSSOPTS/Gauss-Job.py ./
cp $GAUSSOPTS/Gauss-MC11a.py ./
emacs -nw Gauss-Job.py

The SetupProject sets the environment variables to run Gauss with the specific version. That line must always be executed each session before Gauss is run, although it is only needed once per terminal. SomeDirectory is a directory where the options, outputs and log files for the jobs will be stored. You should probably create this using mkdir if you do not have a suitable place already. The emacs line can be replaced with an editor of your choice.

Edit Gauss-Job.py to simulate the required events

Change the 5 in line

nEvts = 5
so that Gauss generates two events.

Change the name of the output file from the default (Gauss writes it for you) to something that will help you remember the content; for example

idFile = 'BsJPsiKs-2evt'
Remeber to uncomment the lines that set the histogram file name and OutputStream file name.

Save the file and quit emacs

Find out the event type for Bd->J/Psi(mumu)Kshort

Go into the EvtGen web page linked from that of DecFiles and find out the event type for the decay you want to simulate. You will want to generate events with the decay products in the acceptance so pick one with DecProdCut in the name (the table is at the end). For Bd->J/Psi(mumu)Kshort this means you will pick 11144103. Clicking on the number of the decay will take you to the official file in python detailing how the decay will be generated.

Run Gauss and look at the monitoring output

To run Gauss you will use the command gaudirun.py and specify that you want

  • a standard Gauss jobs with the MC11a geometry
  • the event type you have chosen
  • what is specific to your job (i.e. number of events, name of output file,...)
You will do so providing the appropriate arguments to the command, as in this case
gaudirun.py Gauss-Job.py Gauss-MC11a.py $DECFILESROOT/options/11144103.py | tee BsJPsiKs-2evt_Gauss.log 

Wait while Gauss configures itself, Pythia, EvtGen and GEANT4, then generates two events.

You should see a very long file produced, the bit where it tells you about the B mesons generated is here:

Generation.SignalRepe...   INFO *************   Signal counters   ****************
Number of events for generator level cut, before : 3, after : 0
Efficiency of the generator level cut : 0 +/- 0
Number of z-inverted events : 2

Number of events for generator anti-particle level cut, before : 3, after : 0
Efficiency of the generator anti-particle level cut : 0 +/- 0

Number of signal B0 in sample : 0 [fraction : 0 +/- 0]
Number of signal B~0 in sample : 2 [fraction : 1 +/- 0]

Number of accepted B0 : 2 [fraction : 0.66667 +/- 0.27217]
Number of accepted B+ : 0 [fraction : 0 +/- 0]
Number of accepted Bs0 : 1 [fraction : 0.33333 +/- 0.27217]
Number of accepted b-Baryon : 0 [fraction : 0 +/- 0]
Number of accepted Bc+ : 0 [fraction : 0 +/- 0]
Number of accepted anti-B0 : 0 [fraction : 0 +/- 0]
Number of accepted B- : 1 [fraction : 1 +/- 0]
Number of accepted anti-Bs0 : 0 [fraction : 0 +/- 0]
Number of accepted anti-b-Baryon : 0 [fraction : 0 +/- 0]
Number of accepted Bc- : 0 [fraction : 0 +/- 0]
Number of accepted (bb) : 0

Number of accepted D0 : 1 [fraction : 0.5 +/- 0.35355]
Number of accepted D+ : 0 [fraction : 0 +/- 0]
Number of accepted Ds+ : 1 [fraction : 0.5 +/- 0.35355]
Number of accepted c-Baryon : 0 [fraction : 0 +/- 0]
Number of accepted anti-D0 : 0 [fraction : nan +/- nan]
Number of accepted D- : 0 [fraction : nan +/- nan]
Number of accepted Ds- : 0 [fraction : nan +/- nan]
Number of accepted anti-c-Baryon : 0 [fraction : nan +/- nan]
Number of accepted (cc) : 0

Number of accepted B : 0 [fraction : 0 +/- 0]
Number of accepted B* : 2 [fraction : 1 +/- 0]
Number of accepted B** : 0 [fraction : 0 +/- 0]

Number of accepted D : 0 [fraction : 0 +/- 0]
Number of accepted D* : 2 [fraction : 1 +/- 0]
Number of accepted D** : 0 [fraction : 0 +/- 0]
showing the number of "signal" particle types created in each event. Note this is not the end of the log file and may have scrolled off the top of the screen. The | tee bit of above command means there is a copy of the output file in BsJPsiKs-2evt.log.

Look at the histogram file produced called something like BsJPsiKs-2evt-histos.root and check there are hits in the VELO, RICH etc.

Digitize the two decays made in Gauss

Setup access to a version of Boole

SetupProject Boole v23r1
cp $BOOLEOPTS/MC11-Files.py ./Boole-MC11a-Files.py

Edit Boole-MC11a-Files.py to setup Boole to digitize the events generated earlier

Change the input data file to the BsJPsiKs-2evt.sim file you generated earlier in Gauss. Also you must always make the DDDBtag and CondDBtag values match as these control the detector geometry and alignment. This must be the same in Gauss and Boole, the values used in Gauss can be found from the file Gauss-MC11a.py which you copied earlier.

So make the following changes to Boole-MC11a-Files.py

#-- Event input
LHCbApp().DDDBtag   = "MC11-20111102"
LHCbApp().CondDBtag = "sim-20111111-vc-md100"

datasetName = 'BsJPsiKs-2evt'
EventSelector().Input = ["DATAFILE='PFN:" + datasetName + ".sim'"]

Run Boole and check the output

Run Boole with gaudirun.py Boole-MC11a-Files.py | tee BsJPsiKs-2evt_Boole.log

Look at the histograms produced by Boole in BsJPsiKs-2evt-histos.root, check there are entries in the histograms.

Use Brunel to reconstruct the events digitized with Boole

Use the same instructions as for Gauss and Boole to setup Brunel version v41r2.

The Brunel options required to run the job are the values of DDDBtag and CondDBtag, the input file name and type and the required output from Brunel.

We can write a file that sets all of this:

from Gaudi.Configuration import *
from Configurables import Brunel, LHCbApp

#-- Event input
LHCbApp().DDDBtag   = "MC11-20111102"
LHCbApp().CondDBtag = "sim-20111111-vc-md100"

datasetName = 'BsJPsiKs-2evt'

EventSelector().Input = ["DATAFILE='PFN:" + datasetName + ".digi?svcClass=default' TYP='POOL_ROOTTREE' OPT='READ'"]

Brunel().DatasetName = datasetName # sets output and histogram file names
Brunel().DataType = '2011'  # sets the 2011 configuration of Brunel
Brunel().InputType = "DIGI" # input has the format digi
Brunel().WithMC    = True   # use the MC truth information in the digi file

Save the above in Brunel-MC11a-Files.py to set up Brunel.

Then run Brunel with gaudirun.py Brunel-MC11a-Files.py | tee BsJPsiKs-2evt_Brunel.log

Again check the output log file and the histograms produced.

-- DavidHutchcroft - 02-Dec-2011

Edit | Attach | Watch | Print version | History: r19 | r16 < r15 < r14 < r13 | Backlinks | Raw View | Raw edit | More topic actions...
Topic revision: r14 - 2011-12-02 - DavidHutchcroft
 
    • Cern Search Icon Cern Search
    • TWiki Search Icon TWiki Search
    • Google Search Icon Google Search

    LHCb All webs login

This site is powered by the TWiki collaboration platform Powered by PerlCopyright & 2008-2019 by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
Ideas, requests, problems regarding TWiki? Send feedback