diff --git a/MC/bin/o2dpg_sim_workflow.py b/MC/bin/o2dpg_sim_workflow.py index 09e9ecd0b..4ee64f500 100755 --- a/MC/bin/o2dpg_sim_workflow.py +++ b/MC/bin/o2dpg_sim_workflow.py @@ -1,5 +1,23 @@ #!/usr/bin/env python3 +# TODO: +# - enforce pregencollcontext (DONE) +# - move pregencollcontext outside of timeframeloop (DONE) +# - remove old way of generating collision context (DONE) +# - make sure vertexing is picked up from the collision context (DONE) +# - background events read from global coll context (DONE) +# - change name of signal sim from sgn_tf to just sgn (DONE) +# - QED handling? (DONE) +# - make several examples needing this: +# -----> (a) ordinary sgn only sim for pp **and** PbPb +# -----> (b) embedding simulation for pp **and** PbPb +# - check than non-trivial embeddings work (ratio 1:3 or the like) (DONE) +# - verify that the right vertex is applied and the right number of events extracted with non-trivial embedding patterns +# - fix getNCollisions in digitization context to be automatic (is N hadronic collisiosn --> neglect QED) +# - fix orbit early functionality with new global pregencollcontext +# - offer option to just have a global transport before any digitization +# - check that the anchoring script works + # # A script producing a consistent MC->RECO->AOD workflow # It aims to handle the different MC possible configurations @@ -55,7 +73,7 @@ parser.add_argument('--orbitsPerTF', type=int, help="Timeframe size in number of LHC orbits", default=128) parser.add_argument('--anchor-config',help="JSON file to contextualise workflow with external configs (config values etc.) for instance comping from data reco workflows.", default='') parser.add_argument('--dump-config',help="Dump JSON file with all settings used in workflow", default='user_config.json') -parser.add_argument('-ns',help='number of signal events / timeframe', default=20) +parser.add_argument('-ns',type=int,help='number of signal events / timeframe', default=20) parser.add_argument('-gen',help='generator: pythia8, extgen', default='') parser.add_argument('-proc',help='process type: inel, dirgamma, jets, ccbar, ...', default='none') parser.add_argument('-trigger',help='event selection: particle, external', default='') @@ -88,13 +106,13 @@ parser.add_argument('-colBkg',help='embedding background collision system', default='PbPb') parser.add_argument('-e',help='simengine', default='TGeant4', choices=['TGeant4', 'TGeant3', 'TFluka']) -parser.add_argument('-tf',help='number of timeframes', default=2) +parser.add_argument('-tf',type=int,help='number of timeframes', default=2) parser.add_argument('--production-offset',help='Offset determining bunch-crossing ' + ' range within a (GRID) production. This number sets first orbit to ' + 'Offset x Number of TimeFrames x OrbitsPerTimeframe (up for further sophistication)', default=0) parser.add_argument('-j', '--n-workers', dest='n_workers', help='number of workers (if applicable)', default=8, type=int) parser.add_argument('--force-n-workers', dest='force_n_workers', action='store_true', help='by default, number of workers is re-computed ' - 'for given interaction rate if --pregenCollContext is set; ' + 'for given interaction rate; ' 'pass this to avoid that') parser.add_argument('-mod',help='Active modules (deprecated)', default='--skipModules ZDC') parser.add_argument('--with-ZDC', action='store_true', help='Enable ZDC in workflow') @@ -111,7 +129,7 @@ # power features (for playing) --> does not appear in help message # help='Treat smaller sensors in a single digitization') -parser.add_argument('--pregenCollContext', action='store_true', help=argparse.SUPPRESS) # the mode where we pregenerate the collision context for each timeframe (experimental) +parser.add_argument('--pregenCollContext', action='store_true', help=argparse.SUPPRESS) # Now the default, giving this option or not makes not difference. We keep it for backward compatibility parser.add_argument('--no-combine-smaller-digi', action='store_true', help=argparse.SUPPRESS) parser.add_argument('--no-combine-dpl-devices', action='store_true', help=argparse.SUPPRESS) parser.add_argument('--no-mc-labels', action='store_true', default=False, help=argparse.SUPPRESS) @@ -356,6 +374,79 @@ def extractVertexArgs(configKeyValuesStr, finalDiamondDict): print ("Using initialisation seed: ", RNDSEED) SIMSEED = random.randint(1, 900000000 - NTIMEFRAMES - 1) # PYTHIA maximum seed is 900M for some reason +# ---- initialize global (physics variables) for signal parts ---- +ECMS=float(args.eCM) +EBEAMA=float(args.eA) +EBEAMB=float(args.eB) +NSIGEVENTS=args.ns +GENERATOR=args.gen +if GENERATOR =='': + print('o2dpg_sim_workflow: Error! generator name not provided') + exit(1) + +INIFILE='' +if args.ini!= '': + INIFILE=' --configFile ' + args.ini +PROCESS=args.proc +TRIGGER='' +if args.trigger != '': + TRIGGER=' -t ' + args.trigger + +## Pt Hat productions +WEIGHTPOW=float(args.weightPow) +PTHATMIN=float(args.ptHatMin) +PTHATMAX=float(args.ptHatMax) + +# translate here collision type to PDG +COLTYPE=args.col + +doembedding=True if args.embedding=='True' or args.embedding==True else False + +if COLTYPE == 'pp': + PDGA=2212 # proton + PDGB=2212 # proton + +if COLTYPE == 'PbPb': + PDGA=1000822080 # Pb + PDGB=1000822080 # Pb + if ECMS < 0: # assign 5.02 TeV to Pb-Pb + print('o2dpg_sim_workflow: Set CM Energy to PbPb case 5.02 TeV') + ECMS=5020.0 + +if COLTYPE == 'pPb': + PDGA=2212 # proton + PDGB=1000822080 # Pb + +if COLTYPE == 'Pbp': + PDGA=1000822080 # Pb + PDGB=2212 # proton + +# If not set previously, set beam energy B equal to A +if EBEAMB < 0 and ECMS < 0: + EBEAMB=EBEAMA + print('o2dpg_sim_workflow: Set beam energy same in A and B beams') + if COLTYPE=="pPb" or COLTYPE=="Pbp": + print('o2dpg_sim_workflow: Careful! both beam energies are the same') + +if ECMS > 0: + if COLTYPE=="pPb" or COLTYPE=="Pbp": + print('o2dpg_sim_workflow: Careful! ECM set for pPb/Pbp collisions!') + +if ECMS < 0 and EBEAMA < 0 and EBEAMB < 0: + print('o2dpg_sim_workflow: Error! CM or Beam Energy not set!!!') + exit(1) + +# Determine interaction rate +INTRATE=int(args.interactionRate) +if INTRATE <= 0: + print('o2dpg_sim_workflow: Error! Interaction rate not >0 !!!') + exit(1) +BCPATTERN=args.bcPatternFile + +# ----- global background specific stuff ------- +COLTYPEBKG=args.colBkg +havePbPb = (COLTYPE == 'PbPb' or (doembedding and COLTYPEBKG == "PbPb")) + workflow={} workflow['stages'] = [] @@ -379,7 +470,6 @@ def getDPL_global_options(bigshm=False, ccdbbackend=True): else: return common -doembedding=True if args.embedding=='True' or args.embedding==True else False usebkgcache=args.use_bkg_from!=None includeFullQC=args.include_qc=='True' or args.include_qc==True includeLocalQC=args.include_local_qc=='True' or args.include_local_qc==True @@ -402,6 +492,42 @@ def getDPL_global_options(bigshm=False, ccdbbackend=True): workflow['stages'].append(GRP_TASK) +includeQED = (COLTYPE == 'PbPb' or (doembedding and COLTYPEBKG == "PbPb")) or (args.with_qed == True) +signalprefix='sgn' + +# preproduce the collision context / timeframe structure for all timeframes at once +precollneeds=[GRP_TASK['name']] +NEventsQED=10000 # max number of QED events to simulate per timeframe +PbPbXSec=8. # expected PbPb cross section +QEDXSecExpected=35237.5 # expected magnitude of QED cross section +PreCollContextTask=createTask(name='precollcontext', needs=precollneeds, cpu='1') + +# adapt timeframeID + orbits + seed + qed +# apply max-collisision offset +# apply vertexing +interactionspecification = signalprefix + ',' + str(INTRATE) + ',' + str(1000000) + ':' + str(1000000) +if doembedding: + interactionspecification = 'bkg,' + str(INTRATE) + ',' + str(NTIMEFRAMES*args.ns) + ':' + str(args.nb) + ' ' + signalprefix + ',' + args.embeddPattern + +PreCollContextTask['cmd']='${O2_ROOT}/bin/o2-steer-colcontexttool -i ' + interactionspecification \ + + ' --show-context ' \ + + ' --timeframeID ' + str(int(args.production_offset)*NTIMEFRAMES) \ + + ' --orbitsPerTF ' + str(orbitsPerTF) \ + + ' --orbits ' + str(NTIMEFRAMES * (orbitsPerTF + math.ceil(args.orbits_early))) \ + + ' --seed ' + str(RNDSEED) \ + + ' --noEmptyTF --first-orbit ' + str(args.first_orbit - args.orbits_early) \ + + ' --extract-per-timeframe tf:sgn' \ + + ' --with-vertices kCCDB' \ + + ' --maxCollsPerTF ' + str(args.ns) + +PreCollContextTask['cmd'] += ' --bcPatternFile ccdb' # <--- the object should have been set in (local) CCDB +if includeQED: + qedrate = INTRATE * QEDXSecExpected / PbPbXSec # hadronic interaction rate * cross_section_ratio + qedspec = 'qed' + ',' + str(qedrate) + ',10000000:' + str(NEventsQED) + PreCollContextTask['cmd'] += ' --QEDinteraction ' + qedspec +workflow['stages'].append(PreCollContextTask) + + if doembedding: if not usebkgcache: # ---- do background transport task ------- @@ -412,7 +538,6 @@ def getDPL_global_options(bigshm=False, ccdbbackend=True): exit(1) PROCESSBKG=args.procBkg - COLTYPEBKG=args.colBkg ECMSBKG=float(args.eCM) EBEAMABKG=float(args.eA) EBEAMBBKG=float(args.eB) @@ -482,12 +607,14 @@ def getDPL_global_options(bigshm=False, ccdbbackend=True): # determine final configKey values for background transport CONFKEYBKG = constructConfigKeyArg(create_geant_config(args, args.confKeyBkg)) - BKGtask=createTask(name='bkgsim', lab=["GEANT"], needs=[BKG_CONFIG_task['name'], GRP_TASK['name']], cpu=NWORKERS ) + bkgsimneeds = [BKG_CONFIG_task['name'], GRP_TASK['name'], PreCollContextTask['name']] + BKGtask=createTask(name='bkgsim', lab=["GEANT"], needs=bkgsimneeds, cpu=NWORKERS) BKGtask['cmd']='${O2_ROOT}/bin/o2-sim -e ' + SIMENGINE + ' -j ' + str(NWORKERS) + ' -n ' + str(NBKGEVENTS) \ + ' -g ' + str(GENBKG) + ' ' + str(MODULES) + ' -o bkg ' + str(INIBKG) \ + ' --field ccdb ' + str(CONFKEYBKG) \ + ('',' --timestamp ' + str(args.timestamp))[args.timestamp!=-1] + ' --run ' + str(args.run) \ - + ' --vertexMode kCCDB' + + ' --vertexMode kCCDB' \ + + ' --fromCollContext collisioncontext.root:bkg' if not isActive('all'): BKGtask['cmd'] += ' --readoutDetectors ' + " ".join(activeDetectors) @@ -568,125 +695,34 @@ def getDPL_global_options(bigshm=False, ccdbbackend=True): for tf in range(1, NTIMEFRAMES + 1): TFSEED = SIMSEED + tf print("Timeframe " + str(tf) + " seed: ", TFSEED) - timeframeworkdir='tf'+str(tf) - # ---- transport task ------- - # function encapsulating the signal sim part - # first argument is timeframe id - ECMS=float(args.eCM) - EBEAMA=float(args.eA) - EBEAMB=float(args.eB) - NSIGEVENTS=args.ns - GENERATOR=args.gen - if GENERATOR =='': - print('o2dpg_sim_workflow: Error! generator name not provided') - exit(1) - - INIFILE='' - if args.ini!= '': - INIFILE=' --configFile ' + args.ini - PROCESS=args.proc - TRIGGER='' - if args.trigger != '': - TRIGGER=' -t ' + args.trigger - - ## Pt Hat productions - WEIGHTPOW=float(args.weightPow) - PTHATMIN=float(args.ptHatMin) - PTHATMAX=float(args.ptHatMax) - - # translate here collision type to PDG - COLTYPE=args.col - havePbPb = (COLTYPE == 'PbPb' or (doembedding and COLTYPEBKG == "PbPb")) - - if COLTYPE == 'pp': - PDGA=2212 # proton - PDGB=2212 # proton - - if COLTYPE == 'PbPb': - PDGA=1000822080 # Pb - PDGB=1000822080 # Pb - if ECMS < 0: # assign 5.02 TeV to Pb-Pb - print('o2dpg_sim_workflow: Set CM Energy to PbPb case 5.02 TeV') - ECMS=5020.0 - - if COLTYPE == 'pPb': - PDGA=2212 # proton - PDGB=1000822080 # Pb - - if COLTYPE == 'Pbp': - PDGA=1000822080 # Pb - PDGB=2212 # proton - - # If not set previously, set beam energy B equal to A - if EBEAMB < 0 and ECMS < 0: - EBEAMB=EBEAMA - print('o2dpg_sim_workflow: Set beam energy same in A and B beams') - if COLTYPE=="pPb" or COLTYPE=="Pbp": - print('o2dpg_sim_workflow: Careful! both beam energies are the same') - - if ECMS > 0: - if COLTYPE=="pPb" or COLTYPE=="Pbp": - print('o2dpg_sim_workflow: Careful! ECM set for pPb/Pbp collisions!') - - if ECMS < 0 and EBEAMA < 0 and EBEAMB < 0: - print('o2dpg_sim_workflow: Error! CM or Beam Energy not set!!!') - exit(1) - - # Determine interaction rate - signalprefix='sgn_' + str(tf) - INTRATE=int(args.interactionRate) - if INTRATE <= 0: - print('o2dpg_sim_workflow: Error! Interaction rate not >0 !!!') - exit(1) - BCPATTERN=args.bcPatternFile - includeQED = (COLTYPE == 'PbPb' or (doembedding and COLTYPEBKG == "PbPb")) or (args.with_qed == True) - - # preproduce the collision context - precollneeds=[GRP_TASK['name']] - NEventsQED=10000 # max number of QED events to simulate per timeframe - PbPbXSec=8. # expected PbPb cross section - QEDXSecExpected=35237.5 # expected magnitude of QED cross section - PreCollContextTask=createTask(name='precollcontext_' + str(tf), needs=precollneeds, tf=tf, cwd=timeframeworkdir, cpu='1') - PreCollContextTask['cmd']='${O2_ROOT}/bin/o2-steer-colcontexttool -i ' + signalprefix + ',' + str(INTRATE) + ',' + str(args.ns) + ':' + str(args.ns) \ - + ' --show-context ' + ' --timeframeID ' + str(tf-1 + int(args.production_offset)*NTIMEFRAMES) \ - + ' --orbitsPerTF ' + str(orbitsPerTF) + ' --orbits ' + str(orbitsPerTF + math.ceil(args.orbits_early)) \ - + ' --seed ' + str(TFSEED) + ' --noEmptyTF --first-orbit ' + str(args.first_orbit - args.orbits_early) - PreCollContextTask['cmd'] += ' --bcPatternFile ccdb' # <--- the object should have been set in (local) CCDB - if includeQED: - qedrate = INTRATE * QEDXSecExpected / PbPbXSec # hadronic interaction rate * cross_section_ratio - qedspec = 'qed_' + str(tf) + ',' + str(qedrate) + ',10000000:' + str(NEventsQED) - PreCollContextTask['cmd'] += ' --QEDinteraction ' + qedspec - workflow['stages'].append(PreCollContextTask) - + # ---- transport task ------- # produce QED background for PbPb collissions QEDdigiargs = "" if includeQED: NEventsQED=10000 # 35K for a full timeframe? - qedneeds=[GRP_TASK['name']] - if args.pregenCollContext == True: - qedneeds.append(PreCollContextTask['name']) + qedneeds=[GRP_TASK['name'], PreCollContextTask['name']] QED_task=createTask(name='qedsim_'+str(tf), needs=qedneeds, tf=tf, cwd=timeframeworkdir, cpu='1') ######################################################################################################## # # ATTENTION: CHANGING THE PARAMETERS/CUTS HERE MIGHT INVALIDATE THE QED INTERACTION RATES USED ELSEWHERE # ######################################################################################################## - QED_task['cmd'] = 'o2-sim -e TGeant3 --field ccdb -j ' + str('1') + ' -o qed_' + str(tf) \ + QED_task['cmd'] = 'o2-sim -e TGeant3 --field ccdb -j ' + str('1') + ' -o qed' \ + ' -n ' + str(NEventsQED) + ' -m PIPE ITS MFT FT0 FV0 FDD ' \ + ('', ' --timestamp ' + str(args.timestamp))[args.timestamp!=-1] + ' --run ' + str(args.run) \ + ' --seed ' + str(TFSEED) \ - + ' -g extgen --configKeyValues \"GeneratorExternal.fileName=$O2_ROOT/share/Generators/external/QEDLoader.C;QEDGenParam.yMin=-7;QEDGenParam.yMax=7;QEDGenParam.ptMin=0.001;QEDGenParam.ptMax=1.;Diamond.width[2]=6.\"' # + (' ',' --fromCollContext collisioncontext.root')[args.pregenCollContext] + + ' -g extgen --configKeyValues \"GeneratorExternal.fileName=$O2_ROOT/share/Generators/external/QEDLoader.C;QEDGenParam.yMin=-7;QEDGenParam.yMax=7;QEDGenParam.ptMin=0.001;QEDGenParam.ptMax=1.;Diamond.width[2]=6.\"' # + ' --fromCollContext collisioncontext.root' QED_task['cmd'] += '; RC=$?; QEDXSecCheck=`grep xSectionQED qedgenparam.ini | sed \'s/xSectionQED=//\'`' QED_task['cmd'] += '; echo "CheckXSection ' + str(QEDXSecExpected) + ' = $QEDXSecCheck"; [[ ${RC} == 0 ]]' # TODO: propagate the Xsecion ratio dynamically - QEDdigiargs=' --simPrefixQED qed_' + str(tf) + ' --qed-x-section-ratio ' + str(QEDXSecExpected/PbPbXSec) + QEDdigiargs=' --simPrefixQED qed' + ' --qed-x-section-ratio ' + str(QEDXSecExpected/PbPbXSec) workflow['stages'].append(QED_task) # recompute the number of workers to increase CPU efficiency - NWORKERS_TF = compute_n_workers(INTRATE, COLTYPE) if (args.pregenCollContext and not args.force_n_workers) else NWORKERS + NWORKERS_TF = compute_n_workers(INTRATE, COLTYPE) if (not args.force_n_workers) else NWORKERS # produce the signal configuration SGN_CONFIG_task=createTask(name='gensgnconf_'+str(tf), tf=tf, cwd=timeframeworkdir) @@ -736,12 +772,10 @@ def getDPL_global_options(bigshm=False, ccdbbackend=True): # transport signals # ----------------- signalneeds=[ SGN_CONFIG_task['name'], GRP_TASK['name'] ] - if (args.pregenCollContext == True): - signalneeds.append(PreCollContextTask['name']) + signalneeds.append(PreCollContextTask['name']) # add embedIntoFile only if embeddPattern does contain a '@' embeddinto= "--embedIntoFile ../bkg_MCHeader.root" if (doembedding & ("@" in args.embeddPattern)) else "" - #embeddinto= "--embedIntoFile ../bkg_MCHeader.root" if doembedding else "" if doembedding: if not usebkgcache: signalneeds = signalneeds + [ BKGtask['name'] ] @@ -779,10 +813,8 @@ def getDPL_global_options(bigshm=False, ccdbbackend=True): + ' --run ' + str(args.run) + ' ' + str(CONFKEY) + str(TRIGGER) \ + ' -g ' + str(GENERATOR) + ' ' + str(INIFILE) + ' -o genevents ' + embeddinto \ + ('', ' --timestamp ' + str(args.timestamp))[args.timestamp!=-1] \ - + ' --seed ' + str(TFSEED) + ' -n ' + str(NSIGEVENTS) - - if args.pregenCollContext == True: - SGNGENtask['cmd'] += ' --fromCollContext collisioncontext.root:' + signalprefix + + ' --seed ' + str(TFSEED) + ' -n ' + str(NSIGEVENTS) \ + + ' --fromCollContext collisioncontext.root:' + signalprefix if GENERATOR=="hepmc": SGNGENtask['cmd'] += "; RC=$?; ${O2DPG_ROOT}/UTILS/UpdateHepMCEventSkip.sh ../HepMCEventSkip.json " + str(tf) + '; [[ ${RC} == 0 ]]' if sep_event_mode == True: @@ -806,8 +838,8 @@ def getDPL_global_options(bigshm=False, ccdbbackend=True): SGNtask['cmd'] = sgncmdbase + ' -g ' + str(GENERATOR) + ' ' + str(TRIGGER) + ' --vertexMode kCCDB ' if not isActive('all'): SGNtask['cmd'] += ' --readoutDetectors ' + " ".join(activeDetectors) - if args.pregenCollContext == True: - SGNtask['cmd'] += ' --fromCollContext collisioncontext.root' + + SGNtask['cmd'] += ' --fromCollContext collisioncontext.root' workflow['stages'].append(SGNtask) # some tasks further below still want geometry + grp in fixed names, so we provide it here @@ -923,23 +955,18 @@ def putConfigValuesNew(listOfMainKeys=[], localCF = {}): PASSNAME='${ALIEN_JDL_LPMANCHORPASSNAME:-unanchored}' # This task creates the basic setup for all digitizers! all digitization configKeyValues need to be given here + # The purpose of this short task is to generate the digi INI file which all other tasks may use contextneeds = [LinkGRPFileTask['name'], SGNtask['name']] if includeQED: contextneeds += [QED_task['name']] ContextTask = createTask(name='digicontext_'+str(tf), needs=contextneeds, tf=tf, cwd=timeframeworkdir, lab=["DIGI"], cpu='1') # this is just to have the digitizer ini file - ContextTask['cmd'] = '${O2_ROOT}/bin/o2-sim-digitizer-workflow --only-context --interactionRate ' + str(INTRATE) \ - + ' ' + getDPL_global_options(ccdbbackend=False) + ' -n ' + str(args.ns) + simsoption \ - + ' --seed ' + str(TFSEED) \ - + ' ' + putConfigValuesNew({"DigiParams.maxOrbitsToDigitize" : str(orbitsPerTF)},{"DigiParams.passName" : str(PASSNAME)}) + ('',' --incontext ' + CONTEXTFILE)[args.pregenCollContext] + QEDdigiargs + ContextTask['cmd'] = '${O2_ROOT}/bin/o2-sim-digitizer-workflow --only-context --interactionRate ' + str(INTRATE) \ + + ' ' + getDPL_global_options(ccdbbackend=False) + ' -n ' + str(args.ns) + simsoption \ + + ' --seed ' + str(TFSEED) \ + + ' ' + putConfigValuesNew({"DigiParams.maxOrbitsToDigitize" : str(orbitsPerTF)},{"DigiParams.passName" : str(PASSNAME)}) \ + + ' --incontext ' + CONTEXTFILE + QEDdigiargs ContextTask['cmd'] += ' --bcPatternFile ccdb' - - # in case of embedding we engineer the context directly and allow the user to provide an embedding pattern - # The :r flag means to shuffle the background events randomly - if doembedding: - ContextTask['cmd'] += ';ln -nfs ../bkg_Kine.root .;${O2_ROOT}/bin/o2-steer-colcontexttool -i bkg,' + str(INTRATE) + ',' + str(args.ns) + ':' + str(args.nb) + ' ' + signalprefix + ',' + args.embeddPattern + ' --show-context ' + ' --timeframeID ' + str(tf-1 + int(args.production_offset)*NTIMEFRAMES) + ' --orbitsPerTF ' + str(orbitsPerTF) + ' --use-existing-kine' - ContextTask['cmd'] += ' --bcPatternFile ccdb --seed ' + str(TFSEED) + ' --orbits ' + str(orbitsPerTF) + ' --noEmptyTF --first-orbit ' + str(args.first_orbit) - workflow['stages'].append(ContextTask) # ===| TPC digi part |===