diff --git a/MC/bin/o2dpg_sim_workflow.py b/MC/bin/o2dpg_sim_workflow.py index 4b330fd10..1bbd30fe6 100755 --- a/MC/bin/o2dpg_sim_workflow.py +++ b/MC/bin/o2dpg_sim_workflow.py @@ -126,6 +126,8 @@ # help='Treat smaller sensors in a single digitization') 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('--data-anchoring', type=str, default='', help="Take collision contexts (from data) stored in this path") +parser.add_argument('--aod-output-folder', type=str, default='', help="Force this AOD folder in the AOD producer") +parser.add_argument('--aod-parent-file', type=str, default='', help="Link this as parent file in the AOD") 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) @@ -1796,6 +1798,10 @@ def getDigiTaskName(det): if created_by_option != '': created_by_option += ' ' + aod_creator + aod_timeframe_id = f"${{ALIEN_PROC_ID}}{aod_df_id}" if not args.run_anchored else "" + if len(args.aod_output_folder) > 0: + aod_timeframe_id = args.aod_output_folder + AODtask = createTask(name='aod_'+str(tf), needs=aodneeds, tf=tf, cwd=timeframeworkdir, lab=["AOD"], mem='4000', cpu='1') AODtask['cmd'] = ('','ln -nfs ../bkg_Kine.root . ;')[doembedding] AODtask['cmd'] += '[ -f AO2D.root ] && rm AO2D.root; ' @@ -1812,12 +1818,13 @@ def getDigiTaskName(det): "--anchor-pass ${ALIEN_JDL_LPMANCHORPASSNAME:-unknown}", "--anchor-prod ${ALIEN_JDL_LPMANCHORPRODUCTION:-unknown}", "--reco-pass ${ALIEN_JDL_LPMPASSNAME:-unknown}", + f"--aod-parent {args.aod_parent_file}", created_by_option, "--combine-source-devices" if not args.no_combine_dpl_devices else "", "--disable-mc" if args.no_mc_labels else "", "--enable-truncation 0" if environ.get("O2DPG_AOD_NOTRUNCATE") or environ.get("ALIEN_JDL_O2DPG_AOD_NOTRUNCATE") else "", "--disable-strangeness-tracker" if args.no_strangeness_tracking else "", - f"--aod-timeframe-id ${{ALIEN_PROC_ID}}{aod_df_id}" if not args.run_anchored else "", + f"--aod-timeframe-id {aod_timeframe_id}" if len(aod_timeframe_id) > 0 else "" ]) # Consider in future: AODtask['disable_alternative_reco_software'] = True # do not apply reco software here (we prefer latest aod converter) workflow['stages'].append(AODtask) diff --git a/MC/config/PWGEM/external/generator/Generator_pythia8_GapTriggered_DYLL.C b/MC/config/PWGEM/external/generator/Generator_pythia8_GapTriggered_DYll.C similarity index 100% rename from MC/config/PWGEM/external/generator/Generator_pythia8_GapTriggered_DYLL.C rename to MC/config/PWGEM/external/generator/Generator_pythia8_GapTriggered_DYll.C diff --git a/MC/config/PWGEM/ini/GeneratorDYee_GapTriggered_Gap5_PbPb5360GeV.ini b/MC/config/PWGEM/ini/GeneratorDYee_GapTriggered_Gap2_OO5360GeV.ini similarity index 72% rename from MC/config/PWGEM/ini/GeneratorDYee_GapTriggered_Gap5_PbPb5360GeV.ini rename to MC/config/PWGEM/ini/GeneratorDYee_GapTriggered_Gap2_OO5360GeV.ini index 8cc7c0ff6..a83045099 100644 --- a/MC/config/PWGEM/ini/GeneratorDYee_GapTriggered_Gap5_PbPb5360GeV.ini +++ b/MC/config/PWGEM/ini/GeneratorDYee_GapTriggered_Gap2_OO5360GeV.ini @@ -1,6 +1,6 @@ [GeneratorExternal] fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/external/generator/Generator_pythia8_GapTriggered_DYll.C -funcName = GeneratorPythia8GapTriggeredDYll(5, 1, 11, -1.5, +1.5, 1000822080, 1000822080, 5360.0) +funcName = GeneratorPythia8GapTriggeredDYll(2, 1, 11, -1.5, +1.5, 1000080160, 1000080160, 5360.0) [GeneratorPythia8] config = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/pythia8/generator/configPythiaEmpty.cfg diff --git a/MC/config/PWGEM/ini/GeneratorDYmumu_GapTriggered_Gap5_PbPb5360GeV.ini b/MC/config/PWGEM/ini/GeneratorDYmumu_GapTriggered_Gap5_PbPb5360GeV.ini deleted file mode 100644 index de8cefbe8..000000000 --- a/MC/config/PWGEM/ini/GeneratorDYmumu_GapTriggered_Gap5_PbPb5360GeV.ini +++ /dev/null @@ -1,7 +0,0 @@ -[GeneratorExternal] -fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/external/generator/Generator_pythia8_GapTriggered_DYll.C -funcName = GeneratorPythia8GapTriggeredDYll(5, 1, 13, -6, -1, 1000822080, 1000822080, 5360.0) - -[GeneratorPythia8] -config = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/pythia8/generator/configPythiaEmpty.cfg -includePartonEvent = true diff --git a/MC/config/PWGEM/ini/GeneratorDYmumu_GapTriggered_Gap5_pp13600GeV.ini b/MC/config/PWGEM/ini/GeneratorDYmumu_GapTriggered_Gap5_pp13600GeV.ini deleted file mode 100644 index e9a7c06da..000000000 --- a/MC/config/PWGEM/ini/GeneratorDYmumu_GapTriggered_Gap5_pp13600GeV.ini +++ /dev/null @@ -1,7 +0,0 @@ -[GeneratorExternal] -fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/external/generator/Generator_pythia8_GapTriggered_DYll.C -funcName = GeneratorPythia8GapTriggeredDYll(5, 1, 13, -6, -1, 2212, 2212, 13600.0) - -[GeneratorPythia8] -config = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGEM/pythia8/generator/configPythiaEmpty.cfg -includePartonEvent = true diff --git a/MC/config/PWGEM/ini/tests/GeneratorDYee_GapTriggered_Gap2_OO5360GeV.C b/MC/config/PWGEM/ini/tests/GeneratorDYee_GapTriggered_Gap2_OO5360GeV.C new file mode 100644 index 000000000..aca7bb427 --- /dev/null +++ b/MC/config/PWGEM/ini/tests/GeneratorDYee_GapTriggered_Gap2_OO5360GeV.C @@ -0,0 +1,51 @@ +int External() { + std::string path{"o2sim_Kine.root"}; + // Check that file exists, can be opened and has the correct tree + TFile file(path.c_str(), "READ"); + if (file.IsZombie()) + { + std::cerr << "Cannot open ROOT file " << path << "\n"; + return 1; + } + auto tree = (TTree *)file.Get("o2sim"); + if (!tree) + { + std::cerr << "Cannot find tree o2sim in file " << path << "\n"; + return 1; + } + std::vector *tracks{}; + tree->SetBranchAddress("MCTrack", &tracks); + + // Check if all events are filled + auto nEvents = tree->GetEntries(); + for (Long64_t i = 0; i < nEvents; ++i) + { + tree->GetEntry(i); + if (tracks->empty()) + { + std::cerr << "Empty entry found at event " << i << "\n"; + return 1; + } + } + // check if each event has at least two oxygen ions + for (int i = 0; i < nEvents; i++) + { + auto check = tree->GetEntry(i); + int count = 0; + for (int idxMCTrack = 0; idxMCTrack < tracks->size(); ++idxMCTrack) + { + auto track = tracks->at(idxMCTrack); + if (track.GetPdgCode() == 1000080160) + { + count++; + } + } + if (count < 2) + { + std::cerr << "Event " << i << " has less than 2 oxygen ions\n"; + return 1; + } + } + + return 0; +} diff --git a/MC/config/PWGEM/ini/tests/GeneratorDYee_GapTriggered_Gap5_OO5360GeV.C b/MC/config/PWGEM/ini/tests/GeneratorDYee_GapTriggered_Gap5_OO5360GeV.C new file mode 100644 index 000000000..aca7bb427 --- /dev/null +++ b/MC/config/PWGEM/ini/tests/GeneratorDYee_GapTriggered_Gap5_OO5360GeV.C @@ -0,0 +1,51 @@ +int External() { + std::string path{"o2sim_Kine.root"}; + // Check that file exists, can be opened and has the correct tree + TFile file(path.c_str(), "READ"); + if (file.IsZombie()) + { + std::cerr << "Cannot open ROOT file " << path << "\n"; + return 1; + } + auto tree = (TTree *)file.Get("o2sim"); + if (!tree) + { + std::cerr << "Cannot find tree o2sim in file " << path << "\n"; + return 1; + } + std::vector *tracks{}; + tree->SetBranchAddress("MCTrack", &tracks); + + // Check if all events are filled + auto nEvents = tree->GetEntries(); + for (Long64_t i = 0; i < nEvents; ++i) + { + tree->GetEntry(i); + if (tracks->empty()) + { + std::cerr << "Empty entry found at event " << i << "\n"; + return 1; + } + } + // check if each event has at least two oxygen ions + for (int i = 0; i < nEvents; i++) + { + auto check = tree->GetEntry(i); + int count = 0; + for (int idxMCTrack = 0; idxMCTrack < tracks->size(); ++idxMCTrack) + { + auto track = tracks->at(idxMCTrack); + if (track.GetPdgCode() == 1000080160) + { + count++; + } + } + if (count < 2) + { + std::cerr << "Event " << i << " has less than 2 oxygen ions\n"; + return 1; + } + } + + return 0; +} diff --git a/MC/config/PWGEM/ini/tests/GeneratorDYee_GapTriggered_Gap5_pp13600GeV.C b/MC/config/PWGEM/ini/tests/GeneratorDYee_GapTriggered_Gap5_pp13600GeV.C new file mode 100644 index 000000000..0332e54fc --- /dev/null +++ b/MC/config/PWGEM/ini/tests/GeneratorDYee_GapTriggered_Gap5_pp13600GeV.C @@ -0,0 +1,51 @@ +int External() { + std::string path{"o2sim_Kine.root"}; + // Check that file exists, can be opened and has the correct tree + TFile file(path.c_str(), "READ"); + if (file.IsZombie()) + { + std::cerr << "Cannot open ROOT file " << path << "\n"; + return 1; + } + auto tree = (TTree *)file.Get("o2sim"); + if (!tree) + { + std::cerr << "Cannot find tree o2sim in file " << path << "\n"; + return 1; + } + std::vector *tracks{}; + tree->SetBranchAddress("MCTrack", &tracks); + + // Check if all events are filled + auto nEvents = tree->GetEntries(); + for (Long64_t i = 0; i < nEvents; ++i) + { + tree->GetEntry(i); + if (tracks->empty()) + { + std::cerr << "Empty entry found at event " << i << "\n"; + return 1; + } + } + // check if each event has at least two oxygen ions + for (int i = 0; i < nEvents; i++) + { + auto check = tree->GetEntry(i); + int count = 0; + for (int idxMCTrack = 0; idxMCTrack < tracks->size(); ++idxMCTrack) + { + auto track = tracks->at(idxMCTrack); + if (track.GetPdgCode() == 2212) + { + count++; + } + } + if (count < 2) + { + std::cerr << "Event " << i << " has less than 2 oxygen ions\n"; + return 1; + } + } + + return 0; +} diff --git a/MC/config/PWGEM/ini/tests/GeneratorDYee_GapTriggered_Gap5_pp5360GeV.C b/MC/config/PWGEM/ini/tests/GeneratorDYee_GapTriggered_Gap5_pp5360GeV.C new file mode 100644 index 000000000..0332e54fc --- /dev/null +++ b/MC/config/PWGEM/ini/tests/GeneratorDYee_GapTriggered_Gap5_pp5360GeV.C @@ -0,0 +1,51 @@ +int External() { + std::string path{"o2sim_Kine.root"}; + // Check that file exists, can be opened and has the correct tree + TFile file(path.c_str(), "READ"); + if (file.IsZombie()) + { + std::cerr << "Cannot open ROOT file " << path << "\n"; + return 1; + } + auto tree = (TTree *)file.Get("o2sim"); + if (!tree) + { + std::cerr << "Cannot find tree o2sim in file " << path << "\n"; + return 1; + } + std::vector *tracks{}; + tree->SetBranchAddress("MCTrack", &tracks); + + // Check if all events are filled + auto nEvents = tree->GetEntries(); + for (Long64_t i = 0; i < nEvents; ++i) + { + tree->GetEntry(i); + if (tracks->empty()) + { + std::cerr << "Empty entry found at event " << i << "\n"; + return 1; + } + } + // check if each event has at least two oxygen ions + for (int i = 0; i < nEvents; i++) + { + auto check = tree->GetEntry(i); + int count = 0; + for (int idxMCTrack = 0; idxMCTrack < tracks->size(); ++idxMCTrack) + { + auto track = tracks->at(idxMCTrack); + if (track.GetPdgCode() == 2212) + { + count++; + } + } + if (count < 2) + { + std::cerr << "Event " << i << " has less than 2 oxygen ions\n"; + return 1; + } + } + + return 0; +} diff --git a/MC/config/PWGEM/ini/tests/GeneratorDYmumu_GapTriggered_Gap5_OO5360GeV.C b/MC/config/PWGEM/ini/tests/GeneratorDYmumu_GapTriggered_Gap5_OO5360GeV.C new file mode 100644 index 000000000..aca7bb427 --- /dev/null +++ b/MC/config/PWGEM/ini/tests/GeneratorDYmumu_GapTriggered_Gap5_OO5360GeV.C @@ -0,0 +1,51 @@ +int External() { + std::string path{"o2sim_Kine.root"}; + // Check that file exists, can be opened and has the correct tree + TFile file(path.c_str(), "READ"); + if (file.IsZombie()) + { + std::cerr << "Cannot open ROOT file " << path << "\n"; + return 1; + } + auto tree = (TTree *)file.Get("o2sim"); + if (!tree) + { + std::cerr << "Cannot find tree o2sim in file " << path << "\n"; + return 1; + } + std::vector *tracks{}; + tree->SetBranchAddress("MCTrack", &tracks); + + // Check if all events are filled + auto nEvents = tree->GetEntries(); + for (Long64_t i = 0; i < nEvents; ++i) + { + tree->GetEntry(i); + if (tracks->empty()) + { + std::cerr << "Empty entry found at event " << i << "\n"; + return 1; + } + } + // check if each event has at least two oxygen ions + for (int i = 0; i < nEvents; i++) + { + auto check = tree->GetEntry(i); + int count = 0; + for (int idxMCTrack = 0; idxMCTrack < tracks->size(); ++idxMCTrack) + { + auto track = tracks->at(idxMCTrack); + if (track.GetPdgCode() == 1000080160) + { + count++; + } + } + if (count < 2) + { + std::cerr << "Event " << i << " has less than 2 oxygen ions\n"; + return 1; + } + } + + return 0; +} diff --git a/MC/config/PWGEM/ini/tests/GeneratorDYmumu_GapTriggered_Gap5_pp5360GeV.C b/MC/config/PWGEM/ini/tests/GeneratorDYmumu_GapTriggered_Gap5_pp5360GeV.C new file mode 100644 index 000000000..0332e54fc --- /dev/null +++ b/MC/config/PWGEM/ini/tests/GeneratorDYmumu_GapTriggered_Gap5_pp5360GeV.C @@ -0,0 +1,51 @@ +int External() { + std::string path{"o2sim_Kine.root"}; + // Check that file exists, can be opened and has the correct tree + TFile file(path.c_str(), "READ"); + if (file.IsZombie()) + { + std::cerr << "Cannot open ROOT file " << path << "\n"; + return 1; + } + auto tree = (TTree *)file.Get("o2sim"); + if (!tree) + { + std::cerr << "Cannot find tree o2sim in file " << path << "\n"; + return 1; + } + std::vector *tracks{}; + tree->SetBranchAddress("MCTrack", &tracks); + + // Check if all events are filled + auto nEvents = tree->GetEntries(); + for (Long64_t i = 0; i < nEvents; ++i) + { + tree->GetEntry(i); + if (tracks->empty()) + { + std::cerr << "Empty entry found at event " << i << "\n"; + return 1; + } + } + // check if each event has at least two oxygen ions + for (int i = 0; i < nEvents; i++) + { + auto check = tree->GetEntry(i); + int count = 0; + for (int idxMCTrack = 0; idxMCTrack < tracks->size(); ++idxMCTrack) + { + auto track = tracks->at(idxMCTrack); + if (track.GetPdgCode() == 2212) + { + count++; + } + } + if (count < 2) + { + std::cerr << "Event " << i << " has less than 2 oxygen ions\n"; + return 1; + } + } + + return 0; +} diff --git a/MC/config/PWGEM/pythia8/generator/pythia8_pp_5360_bbbar_forceddecayscharmbeauty.cfg b/MC/config/PWGEM/pythia8/generator/pythia8_pp_5360_bbbar_forceddecayscharmbeauty.cfg index 375c72768..2cc3486fe 100644 --- a/MC/config/PWGEM/pythia8/generator/pythia8_pp_5360_bbbar_forceddecayscharmbeauty.cfg +++ b/MC/config/PWGEM/pythia8/generator/pythia8_pp_5360_bbbar_forceddecayscharmbeauty.cfg @@ -30,97 +30,73 @@ MultiPartonInteractions:pT0Ref = 2.15 BeamRemnants:remnantMode = 1 BeamRemnants:saturation =5 -### only semileptonic decays from https://pdg.lbl.gov/2025/ +### only semileptonic decays for charm ### D+ -411:oneChannel = 1 0.08810 0 -311 -11 12 -#411:addChannel = 1 0.04020 0 -321 211 -11 12 # may be double counting -411:addChannel = 1 0.05400 0 -313 -11 12 -411:addChannel = 1 0.00372 0 111 -11 12 -411:addChannel = 1 0.00111 0 221 -11 12 -411:addChannel = 1 0.00187 0 113 -11 12 -411:addChannel = 1 0.00169 0 223 -11 12 -411:addChannel = 1 0.00020 0 331 -11 12 +411:oneChannel = 1 0.087 0 -311 -11 12 +411:addChannel = 1 0.040 0 -321 211 -11 12 +411:addChannel = 1 0.037 0 -313 -11 12 ### D0 -421:oneChannel = 1 0.0354 0 -321 -11 12 -421:addChannel = 1 0.0216 0 -323 -11 12 -#421:addChannel = 1 0.0160 0 -321 111 -11 12 # may be double counting -#421:addChannel = 1 0.0144 0 -311 211 -11 12 # may be double counting -421:addChannel = 1 0.00291 0 -211 -11 12 -421:addChannel = 1 0.00145 0 -211 111 -11 12 -421:addChannel = 1 0.00146 0 -213 -11 12 +421:oneChannel = 1 0.035 0 -321 -11 12 +421:addChannel = 1 0.022 0 -323 -11 12 +421:addChannel = 1 0.016 0 -321 111 -11 12 +421:addChannel = 1 0.014 0 -311 -211 -11 12 ### Ds -431:oneChannel = 1 0.0234 0 333 -11 12 -431:addChannel = 1 0.0227 0 221 -11 12 -431:addChannel = 1 0.0081 0 331 -11 12 -431:addChannel = 1 0.0029 0 311 -11 12 -431:addChannel = 1 0.0021 0 313 -11 12 +431:oneChannel = 1 0.025 0 333 -11 12 +431:addChannel = 1 0.027 0 221 -11 12 ### Lambdac -4122:oneChannel = 1 0.0356 0 3122 -11 12 -4122:oneChannel = 1 0.0009 0 2212 -321 -11 12 -4122:oneChannel = 1 0.0010 0 3124 -11 12 -4122:oneChannel = 1 0.0004 0 3126 -11 12 -### Xi_{c}^{+} +4122:oneChannel = 1 0.036 0 3122 -11 12 +### chi_{c}^{+} 4232:oneChannel = 1 0.07 0 3322 -11 12 -### Xi_{c}^{0} -4132:oneChannel = 1 0.0105 0 3312 -11 12 +### chi_{c}^{0} +4132:oneChannel = 1 0.014 0 3312 -11 12 ### Omega_{c} 4332:oneChannel = 1 0.01224 0 3334 -11 12 -### only semileptonic decays for beauty from https://pdg.lbl.gov/2025/ +### only semileptonic decays for beauty ### B0 -511:oneChannel = 1 0.021000 0 12 -11 -411 -511:addChannel = 1 0.048700 0 12 -11 -413 -511:addChannel = 1 0.002300 0 12 -11 -415 -511:addChannel = 1 0.000150 0 12 -11 -211 -511:addChannel = 1 0.000294 0 12 -11 -213 -#511:addChannel = 1 0.004500 0 12 -11 -10411 -#511:addChannel = 1 0.005200 0 12 -11 -10413 -#511:addChannel = 1 0.008300 0 12 -11 -20413 -511:addChannel = 1 0.003640 0 12 -11 -421 -211 -511:addChannel = 1 0.005440 0 12 -11 -423 -211 -511:addChannel = 1 0.001450 0 12 -11 -411 -211 211 -511:addChannel = 1 0.000510 0 12 -11 -413 -211 211 +511:oneChannel = 1 0.0207000 0 12 -11 -411 +511:addChannel = 1 0.0570000 0 12 -11 -413 +511:addChannel = 1 0.0023000 0 12 -11 -415 +511:addChannel = 1 0.0001330 0 12 -11 -211 +511:addChannel = 1 0.0002690 0 12 -11 -213 +511:addChannel = 1 0.0045000 0 12 -11 -10411 +511:addChannel = 1 0.0052000 0 12 -11 -10413 +511:addChannel = 1 0.0083000 0 12 -11 -20413 ### B+ -521:oneChannel = 1 0.000078 0 12 -11 111 -521:addChannel = 1 0.000158 0 12 -11 113 -521:addChannel = 1 0.000035 0 12 -11 221 -521:addChannel = 1 0.000119 0 12 -11 223 -521:addChannel = 1 0.000024 0 12 -11 331 -521:addChannel = 1 0.022600 0 12 -11 -421 -521:addChannel = 1 0.052600 0 12 -11 -423 -521:addChannel = 1 0.001590 0 12 -11 -425 -#521:addChannel = 1 0.000900 0 12 -11 -10421 -#521:addChannel = 1 0.002840 0 12 -11 -10423 -#521:addChannel = 1 0.001700 0 12 -11 -20423 -521:addChannel = 1 0.003820 0 12 -11 -411 211 -521:addChannel = 1 0.005420 0 12 -11 -413 211 -521:addChannel = 1 0.001730 0 12 -11 -421 -211 211 -521:addChannel = 1 0.000700 0 12 -11 -423 -211 211 -521:addChannel = 1 0.000300 0 12 -11 -431 321 -521:addChannel = 1 0.000290 0 12 -11 -433 321 +521:oneChannel = 1 0.0000720 0 12 -11 111 +521:addChannel = 1 0.0001450 0 12 -11 113 +521:addChannel = 1 0.0000840 0 12 -11 221 +521:addChannel = 1 0.0001450 0 12 -11 223 +521:addChannel = 1 0.0000840 0 12 -11 331 +521:addChannel = 1 0.0224000 0 12 -11 -421 +521:addChannel = 1 0.0617000 0 12 -11 -423 +521:addChannel = 1 0.0030000 0 12 -11 -425 +521:addChannel = 1 0.0049000 0 12 -11 -10421 +521:addChannel = 1 0.0056000 0 12 -11 -10423 +521:addChannel = 1 0.0090000 0 12 -11 -20423 ### Bs -531:oneChannel = 1 0.000106 0 12 -11 -321 -531:addChannel = 1 0.000300 0 12 -11 -323 -531:addChannel = 1 0.022900 0 12 -11 -431 -531:addChannel = 1 0.052000 0 12 -11 -433 -531:addChannel = 1 0.007000 0 12 -11 -435 -531:addChannel = 1 0.000300 0 12 -11 -10323 -531:addChannel = 1 0.004000 0 12 -11 -10431 -531:addChannel = 1 0.007000 0 12 -11 -10433 -531:addChannel = 1 0.000200 0 12 -11 -20323 -531:addChannel = 1 0.004000 0 12 -11 -20433 +531:oneChannel = 1 0.0002000 0 12 -11 -321 +531:addChannel = 1 0.0003000 0 12 -11 -323 +531:addChannel = 1 0.0210000 0 12 -11 -431 +531:addChannel = 1 0.0490000 0 12 -11 -433 +531:addChannel = 1 0.0070000 0 12 -11 -435 +531:addChannel = 1 0.0003000 0 12 -11 -10323 +531:addChannel = 1 0.0040000 0 12 -11 -10431 +531:addChannel = 1 0.0070000 0 12 -11 -10433 +531:addChannel = 1 0.0002000 0 12 -11 -20323 +531:addChannel = 1 0.0040000 0 12 -11 -20433 ### Lambdab 5122:oneChannel = 1 0.0546000 0 -12 11 4122 5122:addChannel = 1 0.0096000 0 -12 11 4124 5122:addChannel = 1 0.0128000 0 -12 11 14122 -### Xi_{b}^{-} +### Chi_{b}^{-} 5132:oneChannel = 1 0.1080010 0 -12 11 4 3101 5132:addChannel = 1 0.0020000 0 -12 11 2 3101 -### Xi_{b}^{0} +### Chi_{b}^{0} 5232:oneChannel = 1 0.1080010 0 -12 11 4 3201 5232:addChannel = 1 0.0020000 0 -12 11 2 3201 ### Omega_{b}^{-} @@ -132,4 +108,3 @@ BeamRemnants:saturation =5 4332:tau0 = 0.08000000000 # Correct Lb decay length (wrong in PYTHIA8 decay table) 5122:tau0 = 4.41000e-01 - diff --git a/MC/config/PWGEM/pythia8/generator/pythia8_pp_5360_cr2_forceddecayscharm.cfg b/MC/config/PWGEM/pythia8/generator/pythia8_pp_5360_cr2_forceddecayscharm.cfg index 2ee27cad6..7f5a2bcab 100644 --- a/MC/config/PWGEM/pythia8/generator/pythia8_pp_5360_cr2_forceddecayscharm.cfg +++ b/MC/config/PWGEM/pythia8/generator/pythia8_pp_5360_cr2_forceddecayscharm.cfg @@ -36,39 +36,24 @@ MultiPartonInteractions:pT0Ref = 2.15 BeamRemnants:remnantMode = 1 BeamRemnants:saturation =5 -### only semileptonic decays from https://pdg.lbl.gov/2025/ + ### only semileptonic decays ### D+ -411:oneChannel = 1 0.08810 0 -311 -11 12 -#411:addChannel = 1 0.04020 0 -321 211 -11 12 # may be double counting -411:addChannel = 1 0.05400 0 -313 -11 12 -411:addChannel = 1 0.00372 0 111 -11 12 -411:addChannel = 1 0.00111 0 221 -11 12 -411:addChannel = 1 0.00187 0 113 -11 12 -411:addChannel = 1 0.00169 0 223 -11 12 -411:addChannel = 1 0.00020 0 331 -11 12 +411:oneChannel = 1 0.087 0 -311 -11 12 +411:addChannel = 1 0.040 0 -321 211 -11 12 +411:addChannel = 1 0.037 0 -313 -11 12 ### D0 -421:oneChannel = 1 0.0354 0 -321 -11 12 -421:addChannel = 1 0.0216 0 -323 -11 12 -#421:addChannel = 1 0.0160 0 -321 111 -11 12 # may be double counting -#421:addChannel = 1 0.0144 0 -311 211 -11 12 # may be double counting -421:addChannel = 1 0.00291 0 -211 -11 12 -421:addChannel = 1 0.00145 0 -211 111 -11 12 -421:addChannel = 1 0.00146 0 -213 -11 12 +421:oneChannel = 1 0.035 0 -321 -11 12 +421:addChannel = 1 0.022 0 -323 -11 12 +421:addChannel = 1 0.016 0 -321 111 -11 12 ### Ds -431:oneChannel = 1 0.0234 0 333 -11 12 -431:addChannel = 1 0.0227 0 221 -11 12 -431:addChannel = 1 0.0081 0 331 -11 12 -431:addChannel = 1 0.0029 0 311 -11 12 -431:addChannel = 1 0.0021 0 313 -11 12 +431:oneChannel = 1 0.025 0 333 -11 12 +431:addChannel = 1 0.027 0 221 -11 12 ### Lambdac -4122:oneChannel = 1 0.0356 0 3122 -11 12 -4122:oneChannel = 1 0.0009 0 2212 -321 -11 12 -4122:oneChannel = 1 0.0010 0 3124 -11 12 -4122:oneChannel = 1 0.0004 0 3126 -11 12 -### Xi_{c}^{+} +4122:oneChannel = 1 0.036 0 3122 -11 12 +### chi_{c}^{+} 4232:oneChannel = 1 0.07 0 3322 -11 12 -### Xi_{c}^{0} -4132:oneChannel = 1 0.0105 0 3312 -11 12 +### chi_{c}^{0} +4132:oneChannel = 1 0.014 0 3312 -11 12 ### Omega_{c} 4332:oneChannel = 1 0.01224 0 3334 -11 12 diff --git a/MC/config/PWGHF/external/generator/generator_pythia8_embed_charmnuclei.C b/MC/config/PWGHF/external/generator/generator_pythia8_embed_charmnuclei.C new file mode 100644 index 000000000..9e2667e4f --- /dev/null +++ b/MC/config/PWGHF/external/generator/generator_pythia8_embed_charmnuclei.C @@ -0,0 +1,335 @@ +#include "FairGenerator.h" +#include "Generators/GeneratorPythia8.h" +#include "Pythia8/Pythia.h" +#include "TRandom.h" +#include "TF1.h" +#include "TMath.h" +#include +#include +#include +#include + +#include "Math/Vector3D.h" +#include "Math/Vector4D.h" + +R__ADD_INCLUDE_PATH($O2DPG_MC_CONFIG_ROOT) +#include "MC/config/common/external/generator/CoalescencePythia8.h" + +using namespace Pythia8; + +class GeneratorPythia8HFEmbedCharmNuclei : public o2::eventgen::GeneratorPythia8 +{ + public: + + /// constructor + GeneratorPythia8HFEmbedCharmNuclei(int pdgCode = 2010010020, float lifetime = 1.f, int nCharmNucleiPerEvent = 10, float yMin = -1.f, float yMax = 1.f, float ptMax = 25.f, bool trivialCoal = false, float coalMomentum = 0.2, float fracFromB = 0.f) + { + nNumberOfCharmNucleiPerEvent = nCharmNucleiPerEvent; + mRapidityMinCharmNuclei = yMin; + mRapidityMaxCharmNuclei = yMax; + mPtMaxCharmNuclei = ptMax; + mTrivialCoal = trivialCoal; + mCoalMomentum = coalMomentum; + mFractionFromBeauty = fracFromB; + mPdgCharmNucleus = pdgCode; + mSign = 1; + if (std::abs(mPdgCharmNucleus) == 2010010020) { + mMassCharmNucleus = 3.226f; + } else { + LOG(fatal) << "********** [GeneratorPythia8HFEmbedCharmNuclei] Only c-deuteron (pdg=2010010020) currently supported! Exit **********"; + } + mLifetimeCharmNucleus = lifetime; + mDecayDistr = new TF1("mDecayDistr", "TMath::Exp(-x * 1./[0])", 0., mLifetimeCharmNucleus * 100); + mDecayDistr->SetNpx(10000); + mDecayDistr->SetParameter(0, mLifetimeCharmNucleus); + mDecayDistrLb = new TF1("mDecayDistrLb", "TMath::Exp(-x * 1./[0])", 0., 44.f); + mDecayDistrLb->SetParameter(0, 0.440f); // lifetime of Lambda_b in mm/c + mDecayDistrLb->SetNpx(10000); + mPtDistrLb = new TF1("mPtDistrLb","[0]*x/TMath::Power((1+TMath::Power(x/[1],[3])),[2])",0.,100.); + mPtDistrLb->SetParameters(1000., 6.97355, 3.20721, 1.71678); + mPtDistrLb->SetNpx(10000); + + Print(); + + auto& param = o2::eventgen::GeneratorPythia8Param::Instance(); + LOG(info) << "Init \'GeneratorPythia8HFEmbedCharmNuclei\' with following parameters"; + LOG(info) << param; + if (param.config.empty()) { + LOG(fatal) << "Failed to init \'GeneratorPythia8\': problems with configuration file "; + } + std::string cfg = gSystem->ExpandPathName(param.config.c_str()); + LOG(info) << "GeneratorPythia8HFEmbedCharmNuclei Reading configuration from file: " << cfg; + if (!mPythiaGun.readFile(cfg, true)) { + LOG(fatal) << "Failed to init \'GeneratorPythia8\': problems with configuration file " << cfg; + } + + if (!mPythiaGun.init()) { + LOG(fatal) << "Failed to init \'GeneratorPythia8\': init returned with error"; + } + } + + /// Destructor + ~GeneratorPythia8HFEmbedCharmNuclei() = default; + + /// Print the input + void Print() + { + LOG(info) << "********** GeneratorPythia8HFEmbedCharmNuclei configuration dump **********"; + LOG(info) << Form("* PDG code of charm nuclei to be injected: %d", mPdgCharmNucleus); + LOG(info) << Form("* Mass of charm nuclei to be injected (GeV/c2): %f", mMassCharmNucleus); + LOG(info) << Form("* Lifetime of charm nuclei to be injected (mm): %f", mLifetimeCharmNucleus); + LOG(info) << Form("* Number of charm nuclei injected per event: %d", nNumberOfCharmNucleiPerEvent); + LOG(info) << Form("* Charmed nucleus rapidity: %f - %f", mRapidityMinCharmNuclei, mRapidityMaxCharmNuclei); + LOG(info) << Form("* Charmed nucleus pT max (prompt): %f", mPtMaxCharmNuclei); + LOG(info) << Form("* Trivial coalescence: %d", mTrivialCoal); + LOG(info) << Form("* Coalescence momentum: %f", mCoalMomentum); + LOG(info) << Form("* Fraction from beauty: %f", mFractionFromBeauty); + LOG(info) << "***********************************************************************"; + } + + void setHadronRapidity(float yMin, float yMax) + { + mRapidityMinCharmNuclei = yMin; + mRapidityMaxCharmNuclei = yMax; + }; + + void setUsedSeed(unsigned int seed) + { + mUsedSeed = seed; + }; + + unsigned int getUsedSeed() const + { + return mUsedSeed; + }; + + //__________________________________________________________________ + bool generateEvent() override + { + // we start from an empty event + mPythia.event.reset(); + + // we simulate c-deuteron decays + for (int iCharmNuclei{0}; iCharmNuclei 0) ? mSign = -1 : mSign = 1; + if (nNumberOfCharmNucleiPerEvent % 2 != 0 && iCharmNuclei == nNumberOfCharmNucleiPerEvent - 1) { + if (gRandom->Rndm() < 0.5) { + mSign = 1; + } + } + + int pdgToGen = mPdgCharmNucleus; + float massToGen = mMassCharmNucleus; + float lifetimeToGen = 0.f; + float minRapToGen = mRapidityMinCharmNuclei; + float maxRapToGen = mRapidityMaxCharmNuclei; + bool isFromB = gRandom->Rndm() < mFractionFromBeauty; + // we determine if it's prompt or non-prompt + if (isFromB) { + pdgToGen = 5122; // we generate a Lambda_b and we let it decay into the charmed nucleus, no other beauty hadrons are considered + massToGen = 5.61940f; // mass of Lambda_b (GeV/c2) + lifetimeToGen = mDecayDistrLb->GetRandom(); + minRapToGen *= 2; + maxRapToGen *= 2; + } else { + lifetimeToGen = mDecayDistr->GetRandom(); + } + + auto pt = (!isFromB) ? gRandom->Uniform(0., mPtMaxCharmNuclei) : mPtDistrLb->GetRandom(); + auto y = gRandom->Uniform(minRapToGen, maxRapToGen); + auto phi = gRandom->Uniform(0, TMath::TwoPi()); + auto px = pt * TMath::Cos(phi); + auto py = pt * TMath::Sin(phi); + auto mt = TMath::Sqrt(massToGen * massToGen + pt * pt); + auto pz = mt * TMath::SinH(y); + auto p = TMath::Sqrt(pt * pt + pz * pz); + auto e = TMath::Sqrt(massToGen * massToGen + p * p); + + Particle particle; + particle.id(mSign * pdgToGen); + particle.status(83); + particle.m(massToGen); + particle.px(px); + particle.py(py); + particle.pz(pz); + particle.e(e); + particle.xProd(0.f); + particle.yProd(0.f); + particle.zProd(0.f); + particle.tau(lifetimeToGen); + mPythiaGun.particleData.mayDecay(5122, true); // force decay + mPythiaGun.particleData.mayDecay(mPdgCharmNucleus, true); // force decay + + bool isCoalSuccess{false}; + int nTrials{0}; + while(!isCoalSuccess || nTrials > 1e4) { + mPythiaGun.event.reset(); + mPythiaGun.event.append(particle); + mPythiaGun.moreDecays(); + std::array dausToCoal = {-1, -1}; + std::vector pdgShortLivedResos = {313, 2224, 102134}; + bool isResoFound{false}; + int idxCharmNucleus{-1}; + for (int iPart{0}; iPart= idxCharmNucleus) { + // we need to change the indices of the daughter particles to point to the charmed nucleus + auto dauList = part.daughterList(); + for (auto const& dau : dauList) { + mPythiaGun.event[dau].mother1(idxCharmNucleus); + } + mPythiaGun.event.remove(iPart, iPart, true); + isResoFound=true; + } + } + if (isResoFound) { // we have to reset all the particles as daughters of the charm nucleus + std::vector idxDausCharmNucleus{}; + for (int iPart{0}; iPart newPartList{}; + std::vector idxToRemove{}; + for (int iPart{idxDausCharmNucleus[0]}; iPart updatedMothers{}; + for (int iPart{0}; iPart{1000010020}, mTrivialCoal, mCoalMomentum, dausToCoal[0], dausToCoal[1], 10.); + if (isCoalSuccess) { + restoreEnergyConservation(mPythiaGun.event, idxCharmNucleus); + int offset = mPythia.event.size(); // we need to rescale the indices of mothers and daughters, accounting for the particles that are already appended to the event + for (int iPart{0}; iPart 0) { + part.mother1(mother1 + offset - 1); + } + if (mother2 > 0) { + part.mother2(mother2 + offset - 1); + } + if (daughter1 > 0) { + part.daughter1(daughter1 + offset - 1); + } + if (daughter2 > 0) { + part.daughter2(daughter2 + offset - 1); + } + mPythia.event.append(part); + } + } + nTrials++; + } + } + + return true; + } + + +private: + + void restoreEnergyConservation(Pythia8::Event& event, int idxCharmNucleus, float targetMassTolerance = 1e-5) { + /// In the coalescence, the energy is not conserved, we rescale the momentum of the charmed nuclei daughters to restore it to avoid changes in the invariant mass of the charmed nucleus + + float scale = 1.f; + float invMass{0.f}; + while (abs(invMass - mMassCharmNucleus) > targetMassTolerance) { + ROOT::Math::PxPyPzMVector fourVecCharmNucleus; + for (int iDau{event[idxCharmNucleus].daughter1()}; iDau<=event[idxCharmNucleus].daughter2(); ++iDau) { + auto dau = event[iDau]; + fourVecCharmNucleus += ROOT::Math::PxPyPzMVector(dau.px() * scale, dau.py() * scale, dau.pz() * scale, dau.m()); + } + invMass = fourVecCharmNucleus.M(); + scale *= mMassCharmNucleus / invMass; + } + + for (int iDau{event[idxCharmNucleus].daughter1()}; iDau<=event[idxCharmNucleus].daughter2(); ++iDau) { + event[iDau].px(event[iDau].px() * scale); + event[iDau].py(event[iDau].py() * scale); + event[iDau].pz(event[iDau].pz() * scale); + } + event[idxCharmNucleus].px(event[idxCharmNucleus].px() * scale); + event[idxCharmNucleus].py(event[idxCharmNucleus].py() * scale); + event[idxCharmNucleus].pz(event[idxCharmNucleus].pz() * scale); + } + + // Properties of selection + float mMassCharmNucleus; /// mass of the charmed nucleus + int mPdgCharmNucleus; /// pdg code of the charmed nucleus + float mLifetimeCharmNucleus; /// lifetime of the charmed nucleus + int nNumberOfCharmNucleiPerEvent; /// number of charmed nuclei injected per event + float mRapidityMinCharmNuclei; /// rapidity min of the generated charmed nuclei + float mRapidityMaxCharmNuclei; /// rapidity max of the generated charmed nuclei + float mPtMaxCharmNuclei; /// pT max of the generated charmed nuclei + unsigned int mUsedSeed; /// seed + bool mTrivialCoal; /// if true, the coalescence is done without checking the distance in the phase space of the nucleons + float mCoalMomentum; /// coalescence momentum + Pythia8::Pythia mPythiaGun; /// Gun generator with decay support + TF1* mDecayDistr; /// Lifetime distribution + TF1* mDecayDistrLb; /// Lifetime distribution for Lb + TF1* mPtDistrLb; /// pt distribution for Lb (power-law fit to FONLL) + float mFractionFromBeauty; /// fraction of charmed nuclei coming from beauty hadrons + int mSign; /// sign of the charmed nuclei to be generated, if 0 they are generated with 50% of probability as particle or antiparticle +}; + + +///___________________________________________________________ +FairGenerator *GenerateHFEmbedCDeuteron(float lifetime = 1.f, int nCharmNucleiPerEvent = 10, float yMin = -1.f, float yMax = 1.f, float ptMax = 25.f, bool trivialCoal = false, float coalMomentum = 0.2f, float fracFromB = 0.25f) +{ + auto myGen = new GeneratorPythia8HFEmbedCharmNuclei(2010010020, lifetime, nCharmNucleiPerEvent, yMin, yMax, ptMax, trivialCoal, coalMomentum, fracFromB); + auto seed = (gRandom->TRandom::GetSeed() % 900000000); + myGen->readString("Random:setSeed on"); + myGen->readString("Random:seed " + std::to_string(seed)); + return myGen; +} diff --git a/MC/config/PWGHF/ini/GeneratorHFTrigger_BToDeuteron.ini b/MC/config/PWGHF/ini/GeneratorHFTrigger_BToDeuteron.ini new file mode 100644 index 000000000..95bf68b19 --- /dev/null +++ b/MC/config/PWGHF/ini/GeneratorHFTrigger_BToDeuteron.ini @@ -0,0 +1,8 @@ +#NEV_TEST> 10 +#O2DPG_TEST> ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/ini/tests/GeneratorHFTrigger_BToDeuteron.C +### The external generator derives from GeneratorPythia8. +[GeneratorExternal] +fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/external/generator/generator_pythia8_hfhadron_to_nuclei.C +funcName=generateHFHadToNuclei(3, {521}, {1000010020}, true, 0.5) +[GeneratorPythia8] +config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/pythia8/generator/pythia8_bminus_CRmode2.cfg diff --git a/MC/config/PWGHF/ini/GeneratorHF_CDeuteron_Injected.ini b/MC/config/PWGHF/ini/GeneratorHF_CDeuteron_Injected.ini new file mode 100644 index 000000000..2a574f340 --- /dev/null +++ b/MC/config/PWGHF/ini/GeneratorHF_CDeuteron_Injected.ini @@ -0,0 +1,10 @@ +#NEV_TEST> 100 +### The external generator derives from GeneratorPythia8. +[GeneratorExternal] +fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/external/generator/generator_pythia8_embed_charmnuclei.C +funcName=GenerateHFEmbedCDeuteron(0.5, 10, -1., 1., 25., 0, 0.4, 0.25) + +### includePartonEvent is needed to keep the c-deuteron in the event record, even if there are no partons in the event +[GeneratorPythia8] +config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/pythia8/generator/pythia8_cdeuteron_to_dkpi.cfg +includePartonEvent=true diff --git a/MC/config/PWGHF/ini/GeneratorHF_bbbar_fwd_gap4_natural_inel_Mode2_muTrig.ini b/MC/config/PWGHF/ini/GeneratorHF_bbbar_fwd_gap4_natural_inel_Mode2_muTrig.ini new file mode 100644 index 000000000..b126fceed --- /dev/null +++ b/MC/config/PWGHF/ini/GeneratorHF_bbbar_fwd_gap4_natural_inel_Mode2_muTrig.ini @@ -0,0 +1,8 @@ +#NEV_TEST> 20 +### The external generator derives from GeneratorPythia8. +[GeneratorExternal] +fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/external/generator/generator_pythia8_gaptriggered_hf.C +funcName=GeneratorPythia8GapTriggeredBeauty(4, -5, -1.5, -4.3, -2.2, {13}) +[GeneratorPythia8] +config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/pythia8/generator/pythia8_inel_Mode2.cfg +includePartonEvent=true \ No newline at end of file diff --git a/MC/config/PWGHF/ini/tests/GeneratorHFTrigger_BToDeuteron.C b/MC/config/PWGHF/ini/tests/GeneratorHFTrigger_BToDeuteron.C new file mode 100644 index 000000000..bd56f98e6 --- /dev/null +++ b/MC/config/PWGHF/ini/tests/GeneratorHFTrigger_BToDeuteron.C @@ -0,0 +1,48 @@ +int External() +{ + std::string path{"o2sim_Kine.root"}; + std::vector checkPdgHadron{521}; + std::vector nucleiDauPdg{1000010020}; + + TFile file(path.c_str(), "READ"); + if (file.IsZombie()) + { + std::cerr << "Cannot open ROOT file " << path << "\n"; + return 1; + } + + auto tree = (TTree *)file.Get("o2sim"); + std::vector *tracks{}; + tree->SetBranchAddress("MCTrack", &tracks); + + int nSignals{}, nSignalGoodDecay{}; + auto nEvents = tree->GetEntries(); + + for (int i = 0; i < nEvents; i++) + { + tree->GetEntry(i); + for (auto &track : *tracks) + { + auto pdg = track.GetPdgCode(); + if (std::find(checkPdgHadron.begin(), checkPdgHadron.end(), std::abs(pdg)) != checkPdgHadron.end()) + { + if(std::abs(track.GetRapidity()) > 1.5) continue; + nSignals++; + for (int j{track.getFirstDaughterTrackId()}; j <= track.getLastDaughterTrackId(); ++j) + { + auto pdgDau = tracks->at(j).GetPdgCode(); + if (std::find(nucleiDauPdg.begin(), nucleiDauPdg.end(), std::abs(pdgDau)) != nucleiDauPdg.end()) + { + nSignalGoodDecay++; + } + } + } + } + } + std::cout << "--------------------------------\n"; + std::cout << "# Events: " << nEvents << "\n"; + std::cout << "# signal hadrons: " << nSignals << "\n"; + std::cout << "# signal hadrons decaying into nuclei: " << nSignalGoodDecay << "\n"; + + return 0; +} \ No newline at end of file diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_CDeuteron_Injected.C b/MC/config/PWGHF/ini/tests/GeneratorHF_CDeuteron_Injected.C new file mode 100644 index 000000000..56d22c72a --- /dev/null +++ b/MC/config/PWGHF/ini/tests/GeneratorHF_CDeuteron_Injected.C @@ -0,0 +1,127 @@ +int External() { + std::string path{"o2sim_Kine.root"}; + + int pdgCDeuteron{2010010020}; + int checkNumberOfCDeuteronPerEvent{10}; + float checkLifetimeCDeuteron{0.05f}; + std::map>> checkDecays{ + {2010010020, {{-321, 211, 1000010020}}} // c-deuteron -> K- + pi+ + deuteron + }; + float checkFracCDeuteronFromBeauty{0.25}; + + TFile file(path.c_str(), "READ"); + if (file.IsZombie()) { + std::cerr << "Cannot open ROOT file " << path << "\n"; + return 1; + } + + auto tree = (TTree *)file.Get("o2sim"); + std::vector *tracks{}; + tree->SetBranchAddress("MCTrack", &tracks); + + int nCDeuteron{}, nCDeuteronGoodDecay{}, nCDeuteronFromBeauty{}; + std::array averageLifetimeCDeuteron{0.f, 0.f}; // prompt and non-prompt + float massCDeuteron = 3.226f; + auto nEvents = tree->GetEntries(); + + for (int i = 0; i < nEvents; i++) { + tree->GetEntry(i); + + for (auto &track : *tracks) { + auto pdg = track.GetPdgCode(); + auto absPdg = std::abs(pdg); + // std::cout << "Event " << i << ": found particle with PDG " << pdg << std::endl; + + if (absPdg == pdgCDeuteron) { // found signal + nCDeuteron++; // count signal PDG + + // std::cout << "Event " << i << ": found c-deuteron with PDG " << pdg << std::endl; + + std::vector pdgsDecay{}; + std::vector pdgsDecayAntiPart{}; + if (track.getFirstDaughterTrackId() >= 0 && track.getLastDaughterTrackId() >= 0) { + for (int j{track.getFirstDaughterTrackId()}; j <= track.getLastDaughterTrackId(); ++j) { + // std::cout << "Fetching daughter track with ID " << j << std::endl; + auto pdgDau = tracks->at(j).GetPdgCode(); + // std::cout << "PDG of daughter track " << j << ": " << pdgDau << std::endl; + pdgsDecay.push_back(pdgDau); + if (pdgDau != 333 && pdgDau != 111) { // phi and pi0 are antiparticles of themselves + pdgsDecayAntiPart.push_back(-pdgDau); + } else { + pdgsDecayAntiPart.push_back(pdgDau); + } + } + } + // std::cout << "Daughters fetched" << std::endl; + + auto mother = track.getMotherTrackId(); + bool isFromBeauty{false}; + if (mother >= 0 && std::abs(tracks->at(mother).GetPdgCode()) == 5122) { // check if c-deuteron comes from Lb + nCDeuteronFromBeauty++; + isFromBeauty = true; + } + + auto dauTrack = tracks->at(track.getFirstDaughterTrackId()); + float decayLength = std::sqrt((track.GetStartVertexCoordinatesX() - dauTrack.GetStartVertexCoordinatesX()) * (track.GetStartVertexCoordinatesX() - dauTrack.GetStartVertexCoordinatesX()) + (track.GetStartVertexCoordinatesY() - dauTrack.GetStartVertexCoordinatesY()) * (track.GetStartVertexCoordinatesY() - dauTrack.GetStartVertexCoordinatesY()) + (track.GetStartVertexCoordinatesZ() - dauTrack.GetStartVertexCoordinatesZ()) * (track.GetStartVertexCoordinatesZ() - dauTrack.GetStartVertexCoordinatesZ())); + if (!isFromBeauty) { + averageLifetimeCDeuteron[0] += decayLength * massCDeuteron / track.GetP(); + } else { + averageLifetimeCDeuteron[1] += decayLength * massCDeuteron / track.GetP(); + } + + std::sort(pdgsDecay.begin(), pdgsDecay.end()); + std::sort(pdgsDecayAntiPart.begin(), pdgsDecayAntiPart.end()); + + for (auto &decay : checkDecays[std::abs(pdg)]) { + if (pdgsDecay == decay || pdgsDecayAntiPart == decay) { + nCDeuteronGoodDecay++; + break; + } + } + // std::cout << "Daughters checked " << std::endl; + } + } + } + + averageLifetimeCDeuteron[0] /= nCDeuteron - nCDeuteronFromBeauty; + averageLifetimeCDeuteron[1] /= nCDeuteronFromBeauty; + + std::cout << "--------------------------------\n"; + std::cout << "# Events: " << nEvents << "\n"; + std::cout <<"# signal c-deuteron: " << nCDeuteron << "\n"; + std::cout <<"# signal c-deuteron decaying in the correct channel: " << nCDeuteronGoodDecay << "\n"; + std::cout <<"# signal c-deuteron from beauty: " << nCDeuteronFromBeauty << "\n"; + std::cout <<"Average lifetime of c-deuteron (prompt): " << averageLifetimeCDeuteron[0] << " (cm) \n"; + std::cout <<"Average lifetime of c-deuteron (non-prompt): " << averageLifetimeCDeuteron[1] << " (cm) \n"; + + float numberOfCDeuteronPerEvent = float(nCDeuteron) / nEvents; + float fracCDeuteronGoodDecay = float(nCDeuteronGoodDecay) / nCDeuteron; + float fracCDeuteronFromBeauty = float(nCDeuteronFromBeauty) / nCDeuteron; + + if (std::abs(numberOfCDeuteronPerEvent - checkNumberOfCDeuteronPerEvent) / numberOfCDeuteronPerEvent > 0.05) { // we put some tolerance since the number of generated events is small + std::cerr << "Number of C-deuterons per event " << numberOfCDeuteronPerEvent << " different than expected " << checkNumberOfCDeuteronPerEvent << "\n"; + return 1; + } + + if (fracCDeuteronGoodDecay < 0.95) { // we put some tolerance since the number of generated events is small + std::cerr << "Fraction of signals decaying into the correct channel " << fracCDeuteronGoodDecay << " lower than expected\n"; + return 1; + } + + if (std::abs(fracCDeuteronFromBeauty - checkFracCDeuteronFromBeauty) / checkFracCDeuteronFromBeauty > 0.10) { // we put some tolerance since the number of generated events is small + std::cerr << "Fraction of signals from beauty " << fracCDeuteronFromBeauty << " different than expected " << checkFracCDeuteronFromBeauty << "\n"; + return 1; + } + + if (std::abs(averageLifetimeCDeuteron[0] - checkLifetimeCDeuteron) / checkLifetimeCDeuteron > 0.10) { // we put some tolerance since the number of generated events is small + std::cerr << "Lifetime for prompt c-deuteron " << averageLifetimeCDeuteron[0] << " different than expected " << checkLifetimeCDeuteron << "\n"; + return 1; + } + + if (std::abs(averageLifetimeCDeuteron[1] - checkLifetimeCDeuteron) / checkLifetimeCDeuteron > 0.10) { // we put some tolerance since the number of generated events is small + std::cerr << "Lifetime for non-prompt c-deuteron " << averageLifetimeCDeuteron[1] << " different than expected " << checkLifetimeCDeuteron << "\n"; + return 1; + } + + return 0; +} diff --git a/MC/config/PWGHF/ini/tests/GeneratorHF_bbbar_fwd_gap4_natural_inel_Mode2_muTrig.C b/MC/config/PWGHF/ini/tests/GeneratorHF_bbbar_fwd_gap4_natural_inel_Mode2_muTrig.C new file mode 100644 index 000000000..34072820e --- /dev/null +++ b/MC/config/PWGHF/ini/tests/GeneratorHF_bbbar_fwd_gap4_natural_inel_Mode2_muTrig.C @@ -0,0 +1,103 @@ +int External() { + + int checkPdgDecayMuon = 13; + int checkPdgQuark = 5; + + float ratioTrigger = 1. / 4; // one event triggered out of 4 + + std::string path{"o2sim_Kine.root"}; + + TFile file(path.c_str(), "READ"); + if (file.IsZombie()) { + std::cerr << "Cannot open ROOT file" << path << "\n"; + return 1; + } + + auto tree = (TTree *)file.Get("o2sim"); + if (!tree) { + std::cerr << "Cannot find tree o2sim in file" << path << "\n"; + return 1; + } + + std::vector *tracks{}; + tree->SetBranchAddress("MCTrack", &tracks); + + o2::dataformats::MCEventHeader *eventHeader = nullptr; + tree->SetBranchAddress("MCEventHeader.", &eventHeader); + + int nEventsMB{}; + int nEventsInj{}; + int nQuarks{}; + int nMuons{}; + + int nMuonsInAcceptance{}; + + auto nEvents = tree->GetEntries(); + + for (int i = 0; i < nEvents; i++) { + tree->GetEntry(i); + // check subgenerator information + if (eventHeader->hasInfo(o2::mcgenid::GeneratorProperty::SUBGENERATORID)) { + bool isValid = false; + int subGeneratorId = eventHeader->getInfo( + o2::mcgenid::GeneratorProperty::SUBGENERATORID, isValid); + if (subGeneratorId == 0) { + nEventsMB++; + } else if (subGeneratorId == checkPdgQuark) { + nEventsInj++; + } + } // if event header + + int nmuonsev = 0; + int nmuonsevinacc = 0; + + for (auto &track : *tracks) { + auto pdg = track.GetPdgCode(); + if (std::abs(pdg) == checkPdgQuark) { + nQuarks++; + continue; + } // pdgquark + auto y = track.GetRapidity(); + if (std::abs(pdg) == checkPdgDecayMuon) { + int igmother = track.getMotherTrackId(); + auto gmTrack = (*tracks)[igmother]; + int gmpdg = gmTrack.GetPdgCode(); + if (int(std::abs(gmpdg) / 100.) == 5 || + int(std::abs(gmpdg) / 1000.) == 5) { + nMuons++; + nmuonsev++; + if (-4.3 < y && y < -2.2) { + nMuonsInAcceptance++; + nmuonsevinacc++; + } + } // gmpdg + + } // pdgdecay + + } // loop track + // std::cout << "#muons per event: " << nmuonsev << "\n"; + // std::cout << "#muons in acceptance per event: " << nmuonsev << "\n"; + } // events + + std::cout << "#events: " << nEvents << "\n"; + std::cout << "# MB events: " << nEventsMB << "\n"; + std::cout << Form("# events injected with %d quark pair: ", checkPdgQuark) + << nEventsInj << "\n"; + if (nEventsMB < nEvents * (1 - ratioTrigger) * 0.95 || + nEventsMB > nEvents * (1 - ratioTrigger) * + 1.05) { // we put some tolerance since the number of + // generated events is small + std::cerr << "Number of generated MB events different than expected\n"; + return 1; + } + if (nEventsInj < nEvents * ratioTrigger * 0.95 || + nEventsInj > nEvents * ratioTrigger * 1.05) { + std::cerr << "Number of generated events injected with " << checkPdgQuark + << " different than expected\n"; + return 1; + } + std::cout << "#muons: " << nMuons << "\n"; + std::cout << "#muons in acceptance: " << nMuonsInAcceptance << "\n"; + + return 0; +} // external diff --git a/MC/config/PWGHF/pythia8/generator/pythia8_bminus_CRmode2.cfg b/MC/config/PWGHF/pythia8/generator/pythia8_bminus_CRmode2.cfg new file mode 100644 index 000000000..285f42a79 --- /dev/null +++ b/MC/config/PWGHF/pythia8/generator/pythia8_bminus_CRmode2.cfg @@ -0,0 +1,81 @@ +### beams +Beams:idA 2212 # proton +Beams:idB 2212 # proton +Beams:eCM 13600. # GeV + +### processes +SoftQCD:inelastic on # all inelastic processes + +### decays +ParticleDecays:limitTau0 on +ParticleDecays:tau0Max 10. + +### switching on Pythia Mode2 +ColourReconnection:mode 1 +ColourReconnection:allowDoubleJunRem off +ColourReconnection:m0 0.3 +ColourReconnection:allowJunctions on +ColourReconnection:junctionCorrection 1.20 +ColourReconnection:timeDilationMode 2 +ColourReconnection:timeDilationPar 0.18 +StringPT:sigma 0.335 +StringZ:aLund 0.36 +StringZ:bLund 0.56 +StringFlav:probQQtoQ 0.078 +StringFlav:ProbStoUD 0.2 +StringFlav:probQQ1toQQ0join 0.0275,0.0275,0.0275,0.0275 +MultiPartonInteractions:pT0Ref 2.15 +BeamRemnants:remnantMode 1 +BeamRemnants:saturation 5 +BeamRemnants:allowBeamJunction = off +BeamRemnants:beamJunction = off +ColourReconnection:allowDiquarkJunctionCR = off + +# Correct decay lengths +# B+ +521:tau0 = 0.4914 # PDG: 1.638 ps -> c*tau = 0.4914 mm + +## B+ (521) decay modes – at least 2 nucleons in the final states +521:onMode = off +521:oneChannel = 1 3e-07 0 -2112 -2112 111 2112 2212 +521:addChannel = 1 3e-07 0 -3122 -2212 213 2112 2212 +521:addChannel = 1 3e-07 0 -2212 -2112 111 323 2112 2212 +521:addChannel = 1 3e-07 0 -2212 -2112 -211 211 321 2112 2212 +521:addChannel = 1 2e-07 0 -3122 -2212 211 2112 2212 +521:addChannel = 1 2e-07 0 -3122 -2212 -211 111 111 211 211 221 2112 2212 +521:addChannel = 1 2e-07 0 -3122 -2112 111 2112 2212 +521:addChannel = 1 2e-07 0 -3122 -2112 -211 111 211 2112 2212 +521:addChannel = 1 2e-07 0 -2214 -2212 111 211 321 2112 2212 +521:addChannel = 1 2e-07 0 -2212 -2212 211 211 311 2112 2212 +521:addChannel = 1 2e-07 0 -2212 -2212 113 211 321 2112 2212 +521:addChannel = 1 2e-07 0 -2212 -2212 111 111 211 211 2112 2212 +521:addChannel = 1 2e-07 0 -2212 -2112 211 311 2112 2212 +521:addChannel = 1 2e-07 0 -2212 -2112 113 321 2112 2212 +521:addChannel = 1 2e-07 0 -2212 -2112 111 221 321 2112 2212 +521:addChannel = 1 2e-07 0 -2212 -2112 111 113 211 2112 2212 +521:addChannel = 1 2e-07 0 -2112 -2112 -211 211 311 2112 2212 +521:addChannel = 1 2e-07 0 -2212 -2112 -211 113 211 211 2112 2212 +521:addChannel = 1 2e-07 0 -2212 -2112 -211 111 211 211 2112 2212 +521:addChannel = 1 2e-07 0 -2112 -2112 221 223 2112 2212 +521:addChannel = 1 2e-07 0 -2112 -2112 111 221 2112 2212 +521:onIfMatch = -2112 -2112 111 2112 2212 +521:onIfMatch = -3122 -2212 213 2112 2212 +521:onIfMatch = -2212 -2112 111 323 2112 2212 +521:onIfMatch = -2212 -2112 -211 211 321 2112 2212 +521:onIfMatch = -3122 -2212 211 2112 2212 +521:onIfMatch = -3122 -2212 -211 111 111 211 211 221 2112 2212 +521:onIfMatch = -3122 -2112 111 2112 2212 +521:onIfMatch = -3122 -2112 -211 111 211 2112 2212 +521:onIfMatch = -2214 -2212 111 211 321 2112 2212 +521:onIfMatch = -2212 -2212 211 211 311 2112 2212 +521:onIfMatch = -2212 -2212 113 211 321 2112 2212 +521:onIfMatch = -2212 -2212 111 111 211 211 2112 2212 +521:onIfMatch = -2212 -2112 211 311 2112 2212 +521:onIfMatch = -2212 -2112 113 321 2112 2212 +521:onIfMatch = -2212 -2112 111 221 321 2112 2212 +521:onIfMatch = -2212 -2112 111 113 211 2112 2212 +521:onIfMatch = -2112 -2112 -211 211 311 2112 2212 +521:onIfMatch = -2212 -2112 -211 113 211 211 2112 2212 +521:onIfMatch = -2212 -2112 -211 111 211 211 2112 2212 +521:onIfMatch = -2112 -2112 221 223 2112 2212 +521:onIfMatch = -2112 -2112 111 221 2112 2212 diff --git a/MC/config/PWGHF/pythia8/generator/pythia8_cdeuteron_to_dkpi.cfg b/MC/config/PWGHF/pythia8/generator/pythia8_cdeuteron_to_dkpi.cfg new file mode 100644 index 000000000..788e56fe3 --- /dev/null +++ b/MC/config/PWGHF/pythia8/generator/pythia8_cdeuteron_to_dkpi.cfg @@ -0,0 +1,102 @@ +### author: Fabrizio Grosa (fabrizio.grosa@cern.ch) +### since: May 2026 + +### beams (not relevant) +Beams:idA 2212 # proton +Beams:idB 2212 # proton +Beams:eCM 13600. # GeV + +### processes +SoftQCD:inelastic on # all inelastic processes + +### decays +ParticleDecays:limitTau0 on +ParticleDecays:tau0Max 10. + +### add the c-deuteron +### id:all = name antiName spinType chargeType colType m0 mWidth mMin mMax tau0 +2010010020:all = CDeuteron AntiCDeuteron 1 3 0 3.226 0 3.226 3.226 0.5 +2010010020:mayDecay = on + +# same as Λc+ -> p K- π+ including resonances + a neutron +2010010020:oneChannel = 1 0.14400 0 2112 2212 -321 211 ### cd+ -> n p K- π+ (non-resonant) 3.5% +2010010020:addChannel = 1 0.08100 100 2112 2212 -313 ### cd+ -> n p K*0(892) 1.96% +2010010020:addChannel = 1 0.04500 100 2112 2224 -321 ### cd+ -> n Delta++ K- 1.08% +2010010020:addChannel = 1 0.09000 100 2112 102134 211 ### cd+ -> n Lambda(1520) K- 2.20e-3 + +### K*0(892) -> K- π+ +313:onMode = off +313:onIfAll = 321 211 +### for Λc -> Delta++ K- +2224:onMode = off +2224:onIfAll = 2212 211 +### for Λc -> Lambda(1520) K- +102134:onMode = off +102134:onIfAll = 2212 321 + +# Λb0 -> cd+ decays taken from PYTHIA +# NB: cd+ must always be the last daughter, not to screw up replacements in the generator +# 15% 2 body +5122:oneChannel = 1 0.15 0 2010010020 -2212 ### Λb0 -> cd+ anti-p +# 27% 3 body +5122:addChannel = 1 0.054 100 -2212 111 2010010020 ### Λb0 -> cd+ anti-p pi0 +5122:addChannel = 1 0.027 100 -2212 113 2010010020 ### Λb0 -> cd+ anti-p rho0 +5122:addChannel = 1 0.027 100 -2212 221 2010010020 ### Λb0 -> cd+ anti-p eta +5122:addChannel = 1 0.027 100 -2212 223 2010010020 ### Λb0 -> cd+ anti-p omega +5122:addChannel = 1 0.060 100 -2112 -211 2010010020 ### Λb0 -> cd+ anti-n0 pi- +5122:addChannel = 1 0.030 100 -2112 -213 2010010020 ### Λb0 -> cd+ anti-n0 rho- +5122:addChannel = 1 0.010 100 -2214 111 2010010020 ### Λb0 -> cd+ DeltaBar- pi0 +5122:addChannel = 1 0.005 100 -2214 113 2010010020 ### Λb0 -> cd+ DeltaBar- rho0 +5122:addChannel = 1 0.005 100 -2214 221 2010010020 ### Λb0 -> cd+ DeltaBar- eta +5122:addChannel = 1 0.005 100 -2214 223 2010010020 ### Λb0 -> cd+ DeltaBar- omega +5122:addChannel = 1 0.007 100 -2114 -211 2010010020 ### Λb0 -> cd+ DeltaBar0 pi- +5122:addChannel = 1 0.003 100 -2114 -213 2010010020 ### Λb0 -> cd+ DeltaBar0 rho- +5122:addChannel = 1 0.007 100 -2224 211 2010010020 ### Λb0 -> cd+ DeltaBar-- pi+ +5122:addChannel = 1 0.003 100 -2224 213 2010010020 ### Λb0 -> cd+ DeltaBar-- rho+ +5122:addChannel = 1 0.010 0 -3122 -321 2010010020 ### Λb0 -> cd+ Lambdabar0 K- +# 31% 4 body +5122:addChannel = 1 0.054 100 -2212 111 111 2010010020 ### Λb0 -> cd+ anti-p pi0 pi0 +5122:addChannel = 1 0.027 100 -2212 211 -211 2010010020 ### Λb0 -> cd+ anti-p pi+ pi- +5122:addChannel = 1 0.013 100 -2212 213 -211 2010010020 ### Λb0 -> cd+ anti-p rho+ pi- +5122:addChannel = 1 0.018 100 -2212 113 111 2010010020 ### Λb0 -> cd+ anti-p rho0 pi0 +5122:addChannel = 1 0.009 100 -2212 113 113 2010010020 ### Λb0 -> cd+ anti-p rho0 rho0 +5122:addChannel = 1 0.018 100 -2212 221 111 2010010020 ### Λb0 -> cd+ anti-p eta pi0 +5122:addChannel = 1 0.009 100 -2212 221 113 2010010020 ### Λb0 -> cd+ anti-p eta rho0 +5122:addChannel = 1 0.018 100 -2212 223 111 2010010020 ### Λb0 -> cd+ anti-p omega pi0 +5122:addChannel = 1 0.009 100 -2212 223 113 2010010020 ### Λb0 -> cd+ anti-p omega rho0 +5122:addChannel = 1 0.040 100 -2112 -211 111 2010010020 ### Λb0 -> cd+ anti-n0 pi- pi0 +5122:addChannel = 1 0.020 100 -2112 -211 113 2010010020 ### Λb0 -> cd+ anti-n0 pi- rho0 +5122:addChannel = 1 0.020 100 -2112 -213 111 2010010020 ### Λb0 -> cd+ anti-n0 rho- pi0 +5122:addChannel = 1 0.010 100 -2112 -213 113 2010010020 ### Λb0 -> cd+ anti-n0 rho- rho0 +5122:addChannel = 1 0.010 100 -2214 111 111 2010010020 ### Λb0 -> cd+ DeltaBar- pi0 pi0 +5122:addChannel = 1 0.005 100 -2214 113 111 2010010020 ### Λb0 -> cd+ DeltaBar- rho0 pi0 +5122:addChannel = 1 0.005 100 -2214 221 111 2010010020 ### Λb0 -> cd+ DeltaBar- eta pi0 +5122:addChannel = 1 0.005 100 -2214 223 111 2010010020 ### Λb0 -> cd+ DeltaBar- omega pi0 +5122:addChannel = 1 0.007 100 -2114 -211 111 2010010020 ### Λb0 -> cd+ DeltaBar0 pi- pi0 +5122:addChannel = 1 0.003 100 -2114 -213 111 2010010020 ### Λb0 -> cd+ DeltaBar0 rho- pi0 +5122:addChannel = 1 0.007 100 -2224 211 111 2010010020 ### Λb0 -> cd+ DeltaBar-- pi+ pi0 +5122:addChannel = 1 0.003 100 -2224 213 111 2010010020 ### Λb0 -> cd+ DeltaBar-- rho+ pi0 +5122:addChannel = 1 0.010 0 -3122 -321 111 2010010020 ### Λb0 -> cd+ Lambdabar0 K- pi0 +# 19% 5 body +5122:addChannel = 1 0.033 100 -2212 111 111 111 2010010020 ### Λb0 -> cd+ anti-p pi0 pi0 pi0 +5122:addChannel = 1 0.016 100 -2212 211 -211 111 2010010020 ### Λb0 -> cd+ anti-p pi+ pi- pi0 +5122:addChannel = 1 0.008 100 -2212 213 -211 111 2010010020 ### Λb0 -> cd+ anti-p rho+ pi- pi0 +5122:addChannel = 1 0.012 100 -2212 113 111 111 2010010020 ### Λb0 -> cd+ anti-p rho0 pi0 pi0 +5122:addChannel = 1 0.006 100 -2212 113 113 111 2010010020 ### Λb0 -> cd+ anti-p rho0 rho0 pi0 +5122:addChannel = 1 0.012 100 -2212 221 111 111 2010010020 ### Λb0 -> cd+ anti-p eta pi0 pi0 +5122:addChannel = 1 0.006 100 -2212 221 113 111 2010010020 ### Λb0 -> cd+ anti-p eta rho0 pi0 +5122:addChannel = 1 0.012 100 -2212 223 111 111 2010010020 ### Λb0 -> cd+ anti-p omega pi0 pi0 +5122:addChannel = 1 0.006 100 -2212 223 113 111 2010010020 ### Λb0 -> cd+ anti-p omega rho0 pi0 +5122:addChannel = 1 0.030 100 -2112 -211 111 111 2010010020 ### Λb0 -> cd+ anti-n0 pi- pi0 pi0 +5122:addChannel = 1 0.011 100 -2112 -211 113 111 2010010020 ### Λb0 -> cd+ anti-n0 pi- rho0 pi0 +5122:addChannel = 1 0.011 100 -2112 -213 111 111 2010010020 ### Λb0 -> cd+ anti-n0 rho- pi0 pi0 +5122:addChannel = 1 0.006 100 -2112 -213 113 111 2010010020 ### Λb0 -> cd+ anti-n0 rho- rho0 pi0 +5122:addChannel = 1 0.006 100 -2214 111 111 111 2010010020 ### Λb0 -> cd+ DeltaBar- pi0 pi0 pi0 +5122:addChannel = 1 0.002 100 -2214 113 111 111 2010010020 ### Λb0 -> cd+ DeltaBar- rho0 pi0 pi0 +5122:addChannel = 1 0.002 100 -2214 221 111 111 2010010020 ### Λb0 -> cd+ DeltaBar- eta pi0 pi0 +5122:addChannel = 1 0.002 100 -2214 223 111 111 2010010020 ### Λb0 -> cd+ DeltaBar- omega pi0 pi0 +5122:addChannel = 1 0.005 100 -2114 -211 111 111 2010010020 ### Λb0 -> cd+ DeltaBar0 pi- pi0 pi0 +5122:addChannel = 1 0.001 100 -2114 -213 111 111 2010010020 ### Λb0 -> cd+ DeltaBar0 rho- pi0 pi0 +5122:addChannel = 1 0.005 100 -2224 211 111 111 2010010020 ### Λb0 -> cd+ DeltaBar-- pi+ pi0 pi0 +5122:addChannel = 1 0.001 100 -2224 213 111 111 2010010020 ### Λb0 -> cd+ DeltaBar-- rho+ pi0 pi0 +5122:addChannel = 1 0.006 0 -3122 -321 111 111 2010010020 ### Λb0 -> cd+ Lambdabar0 K- pi0 pi0 diff --git a/MC/config/PWGUD/external/generator/GeneratorStarlight.C b/MC/config/PWGUD/external/generator/GeneratorStarlight.C index b2069c16b..be7e4d367 100644 --- a/MC/config/PWGUD/external/generator/GeneratorStarlight.C +++ b/MC/config/PWGUD/external/generator/GeneratorStarlight.C @@ -337,6 +337,16 @@ class GeneratorStarlight_class : public Generator return true; } + //Store couple of event informations to the MC event header + void updateHeader(o2::dataformats::MCEventHeader* eventHeader)override + { + using Key = o2::dataformats::MCInfoKeys; + + eventHeader->putInfo(Key::generator, "STARlight"); + eventHeader->putInfo(Key::processName, mSelectedConfiguration); + eventHeader->putInfo("photonEnergy", getPhotonEnergy()); + } + protected: float eCM = 5020; //CMS energy int projA=208; //Beam diff --git a/MC/config/PWGUD/external/generator/Generator_nOOn.C b/MC/config/PWGUD/external/generator/Generator_nOOn.C index dc19f9504..798c0d6b2 100644 --- a/MC/config/PWGUD/external/generator/Generator_nOOn.C +++ b/MC/config/PWGUD/external/generator/Generator_nOOn.C @@ -1,3 +1,4 @@ +R__ADD_INCLUDE_PATH($nOOn_ROOT/include) R__LOAD_LIBRARY(NeutronGenerator_cxx.so) #include "GeneratorStarlight.C" #include "NeutronGenerator.h" diff --git a/MC/config/common/external/generator/CoalescencePythia8.h b/MC/config/common/external/generator/CoalescencePythia8.h index b1f392ba4..7e4abc53e 100644 --- a/MC/config/common/external/generator/CoalescencePythia8.h +++ b/MC/config/common/external/generator/CoalescencePythia8.h @@ -104,7 +104,7 @@ bool doCoal(Pythia8::Event& event, int charge, int pdgCode, float mass, bool tri return true; } -bool CoalescencePythia8(Pythia8::Event& event, std::vector inputPdgList = {}, bool trivialCoal = false, double coalMomentum = 0.4, int firstDauID = -1, int lastDauId = -1) +bool CoalescencePythia8(Pythia8::Event& event, std::vector inputPdgList = {}, bool trivialCoal = false, double coalMomentum = 0.4, int firstDauID = -1, int lastDauId = -1, float maxRapidity = 1.) { const double coalescenceRadius{0.5 * 1.122462 * coalMomentum}; // if coalescence from a heavy hadron, loop only between firstDauID and lastDauID @@ -131,7 +131,7 @@ bool CoalescencePythia8(Pythia8::Event& event, std::vector inputPd // fill nucleon pools std::vector protons[2], neutrons[2], lambdas[2]; for (auto iPart{loopStart}; iPart <= loopEnd; ++iPart) { - if (std::abs(event[iPart].y()) > 1.) // skip particles with y > 1 + if (std::abs(event[iPart].y()) > maxRapidity) // skip particles with y > ymax { continue; } diff --git a/MC/config/common/ini/pythia8_PbPb_rescattering_536.ini b/MC/config/common/ini/pythia8_PbPb_rescattering_536.ini new file mode 100644 index 000000000..1ed1278ef --- /dev/null +++ b/MC/config/common/ini/pythia8_PbPb_rescattering_536.ini @@ -0,0 +1,10 @@ +#NEV_TEST> 1 +[Diamond] +width[2]=6.0 + +[GeneratorExternal] +fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/ALICE3/pythia8/generator_pythia8_ALICE3.C +funcName=generator_pythia8_ALICE3() + +[GeneratorPythia8] +config=${O2DPG_MC_CONFIG_ROOT}/MC/config/common/pythia8/generator/pythia8_PbPb_rescattering_536.cfg diff --git a/MC/config/common/ini/tests/pythia8_PbPb_rescattering_536.C b/MC/config/common/ini/tests/pythia8_PbPb_rescattering_536.C new file mode 100644 index 000000000..1c28040e2 --- /dev/null +++ b/MC/config/common/ini/tests/pythia8_PbPb_rescattering_536.C @@ -0,0 +1,24 @@ +int External() { + std::string path{"o2sim_Kine.root"}; + + TFile file(path.c_str(), "READ"); + if (file.IsZombie()) { + std::cerr << "Cannot open ROOT file " << path << "\n"; + return 1; + } + + auto tree = (TTree *)file.Get("o2sim"); + if (!tree) { + std::cerr << "Cannot find tree o2sim in file " << path << "\n"; + return 1; + } + std::vector *tracks{}; + tree->SetBranchAddress("MCTrack", &tracks); + + auto nEvents = tree->GetEntries(); + if (nEvents == 0) { + std::cerr << "No event of interest\n"; + return 1; + } + return 0; +} diff --git a/MC/config/common/pythia8/generator/pythia8_PbPb_rescattering_536.cfg b/MC/config/common/pythia8/generator/pythia8_PbPb_rescattering_536.cfg new file mode 100644 index 000000000..d4158bdcf --- /dev/null +++ b/MC/config/common/pythia8/generator/pythia8_PbPb_rescattering_536.cfg @@ -0,0 +1,24 @@ +### OO beams +Beams:idA 1000822080 # Pb +Beams:idB 1000822080 # Pb +Beams:eCM = 5360.0 ### energy + +Beams:frameType = 1 +ParticleDecays:limitTau0 = on +ParticleDecays:tau0Max = 10. ### match alice: 1cm/c = 10.0mm/c + +### Initialize the Angantyr model to fit the total and semi-includive +### cross sections in Pythia within some tolerance. +HeavyIon:SigFitErr = 0.02,0.02,0.1,0.05,0.05,0.0,0.1,0.0 + +### These parameters are typicall suitable for sqrt(S_NN)=5TeV +HeavyIon:SigFitDefPar = 17.24,2.15,0.33,0.0,0.0,0.0,0.0,0.0 + +### enable hadronic rescattering +HadronLevel:Rescatter = on # default = off +Fragmentation:setVertices = on # default = off +PartonVertex:setVertex = on # default = off +Rescattering:nearestNeighbours = off # default = on (but "require a larger retuning effort") +Rescattering:inelastic = on # default = on + +Random:setSeed = on diff --git a/MC/run/ANCHOR/anchorMC_DataEmbedding.sh b/MC/run/ANCHOR/anchorMC_DataEmbedding.sh index 6108e59eb..80cb1daae 100755 --- a/MC/run/ANCHOR/anchorMC_DataEmbedding.sh +++ b/MC/run/ANCHOR/anchorMC_DataEmbedding.sh @@ -329,7 +329,12 @@ for external_context in collission_context_*.root; do # extract timeframe from name anchoring_tf="${external_context#collission_context_}" # remove prefix 'collision_context_' anchoring_tf="${anchoring_tf%.root}" # remove suffix '.root' - echo "Treating timeframe ${anchoring_tf}" + # now we have a string with DF_FOLDERNAME:TIMEFRAMEID + + df_folder="${anchoring_tf%%:*}" + anchoring_tf="${anchoring_tf##*:}" + + echo "Treating timeframe ${anchoring_tf} coming from folder ${df_folder}" # we do it in a separate workspace workspace="TF_${anchoring_tf}" @@ -366,7 +371,7 @@ for external_context in collission_context_*.root; do remainingargs="${remainingargs} ${ALIEN_JDL_CCDB_CONDITION_NOT_AFTER:+--condition-not-after ${ALIEN_JDL_CCDB_CONDITION_NOT_AFTER}}" # add external collision context injection if [ "${ALIEN_JDL_MC_DATA_EMBEDDING_AO2D}" ]; then - remainingargs="${remainingargs} --data-anchoring ${PWD}/../${external_context}" + remainingargs="${remainingargs} --data-anchoring ${PWD}/../${external_context} --aod-output-folder ${df_folder##*_} --aod-parent-file ${ALIEN_JDL_MC_DATA_EMBEDDING_AO2D}" fi echo_info "baseargs passed to o2dpg_sim_workflow_anchored.py: ${baseargs}" diff --git a/MC/run/PWGEM/runAnchoredDYee_OO_5360_Gap2.sh b/MC/run/PWGEM/runAnchoredDYee_OO_5360_Gap2.sh new file mode 100644 index 000000000..9b36bcf51 --- /dev/null +++ b/MC/run/PWGEM/runAnchoredDYee_OO_5360_Gap2.sh @@ -0,0 +1,52 @@ +#!/bin/bash + +# +# Steering script for HF enhanced dielectron MC anchored to LHC25ae +# + +# example anchoring + +export ALIEN_JDL_LPMANCHORPASSNAME=apass2 +export ALIEN_JDL_MCANCHOR=apass2 +export ALIEN_JDL_CPULIMIT=8 +export ALIEN_JDL_LPMRUNNUMBER=564356 +export ALIEN_JDL_LPMPRODUCTIONTYPE=MC +export ALIEN_JDL_LPMINTERACTIONTYPE=pp +export ALIEN_JDL_LPMPRODUCTIONTAG=LHC26a2 +export ALIEN_JDL_LPMANCHORRUN=564356 +export ALIEN_JDL_LPMANCHORPRODUCTION=LHC25ae +export ALIEN_JDL_LPMANCHORYEAR=2025 +export ALIEN_JDL_OUTPUT=*.dat@disk=1,*.txt@disk=1,*.root@disk=2 + +export NTIMEFRAMES=1 +export NSIGEVENTS=20 +export SPLITID=100 +export PRODSPLIT=153 +export CYCLE=0 + +# on the GRID, this is set and used as seed; when set, it takes precedence over SEED +#export ALIEN_PROC_ID=2963436952 +export SEED=0 + +# for pp and 50 events per TF, we launch only 4 workers. +export NWORKERS=2 + +# define the generator via ini file +# use 30/70 sampling for different generators +# No forced beauty decays as we have observed biases + +# generate random number +RNDSIG=$(($RANDOM % 100)) + +CONFIGNAME="GeneratorDYee_GapTriggered_Gap2_OO5360GeV.ini" + +export ALIEN_JDL_ANCHOR_SIM_OPTIONS="-gen external -ini $O2DPG_ROOT/MC/config/PWGEM/ini/$CONFIGNAME" + +# run the central anchor steering script; this includes +# * derive timestamp +# * derive interaction rate +# * extract and prepare configurations (which detectors are contained in the run etc.) +# * run the simulation (and QC) +# To disable QC, uncomment the following line +#export DISABLE_QC=1 +${O2DPG_ROOT}/MC/run/ANCHOR/anchorMC.sh diff --git a/MC/run/PWGEM/runAnchoredDYee_OO_5360_Gap5.sh b/MC/run/PWGEM/runAnchoredDYee_OO_5360_Gap5.sh new file mode 100644 index 000000000..8767cae2d --- /dev/null +++ b/MC/run/PWGEM/runAnchoredDYee_OO_5360_Gap5.sh @@ -0,0 +1,52 @@ +#!/bin/bash + +# +# Steering script for HF enhanced dielectron MC anchored to LHC25ae +# + +# example anchoring + +export ALIEN_JDL_LPMANCHORPASSNAME=apass2 +export ALIEN_JDL_MCANCHOR=apass2 +export ALIEN_JDL_CPULIMIT=8 +export ALIEN_JDL_LPMRUNNUMBER=564356 +export ALIEN_JDL_LPMPRODUCTIONTYPE=MC +export ALIEN_JDL_LPMINTERACTIONTYPE=pp +export ALIEN_JDL_LPMPRODUCTIONTAG=LHC26a2 +export ALIEN_JDL_LPMANCHORRUN=564356 +export ALIEN_JDL_LPMANCHORPRODUCTION=LHC25ae +export ALIEN_JDL_LPMANCHORYEAR=2025 +export ALIEN_JDL_OUTPUT=*.dat@disk=1,*.txt@disk=1,*.root@disk=2 + +export NTIMEFRAMES=1 +export NSIGEVENTS=20 +export SPLITID=100 +export PRODSPLIT=153 +export CYCLE=0 + +# on the GRID, this is set and used as seed; when set, it takes precedence over SEED +#export ALIEN_PROC_ID=2963436952 +export SEED=0 + +# for pp and 50 events per TF, we launch only 4 workers. +export NWORKERS=2 + +# define the generator via ini file +# use 30/70 sampling for different generators +# No forced beauty decays as we have observed biases + +# generate random number +RNDSIG=$(($RANDOM % 100)) + +CONFIGNAME="GeneratorDYee_GapTriggered_Gap5_OO5360GeV.ini" + +export ALIEN_JDL_ANCHOR_SIM_OPTIONS="-gen external -ini $O2DPG_ROOT/MC/config/PWGEM/ini/$CONFIGNAME" + +# run the central anchor steering script; this includes +# * derive timestamp +# * derive interaction rate +# * extract and prepare configurations (which detectors are contained in the run etc.) +# * run the simulation (and QC) +# To disable QC, uncomment the following line +#export DISABLE_QC=1 +${O2DPG_ROOT}/MC/run/ANCHOR/anchorMC.sh diff --git a/MC/run/PWGEM/runAnchoredDYee_pp_13600_Gap5.sh b/MC/run/PWGEM/runAnchoredDYee_pp_13600_Gap5.sh new file mode 100644 index 000000000..77450c320 --- /dev/null +++ b/MC/run/PWGEM/runAnchoredDYee_pp_13600_Gap5.sh @@ -0,0 +1,50 @@ +#!/bin/bash + +# +# Steering script for HF enhanced dielectron MC anchored to LHC24ap,aq apass1 +# + +# example anchoring +# taken from https://its.cern.ch/jira/browse/O2-5670 +export ALIEN_JDL_LPMANCHORPASSNAME=apass1 +export ALIEN_JDL_MCANCHOR=apass1 +export ALIEN_JDL_CPULIMIT=8 +export ALIEN_JDL_LPMRUNNUMBER=559348 +export ALIEN_JDL_LPMPRODUCTIONTYPE=MC +export ALIEN_JDL_LPMINTERACTIONTYPE=pp +export ALIEN_JDL_LPMPRODUCTIONTAG=LHC24a6 +export ALIEN_JDL_LPMANCHORRUN=559348 +export ALIEN_JDL_LPMANCHORPRODUCTION=LHC26ac +export ALIEN_JDL_LPMANCHORYEAR=2026 +export ALIEN_JDL_OUTPUT=*.dat@disk=1,*.txt@disk=1,*.root@disk=2 + +export NTIMEFRAMES=1 +export NSIGEVENTS=20 +export SPLITID=100 +export PRODSPLIT=153 +export CYCLE=0 + +# on the GRID, this is set and used as seed; when set, it takes precedence over SEED +#export ALIEN_PROC_ID=2963436952 +export SEED=0 + +# for pp and 50 events per TF, we launch only 4 workers. +export NWORKERS=2 + +# define the generator via ini file +# use 20/40/40 sampling for different generators +# generate random number +#RNDSIG=$(($RANDOM % 100)) + +CONFIGNAME="GeneratorDYee_GapTriggered_Gap5_pp13600GeV.ini" + +export ALIEN_JDL_ANCHOR_SIM_OPTIONS="-gen external -ini $O2DPG_ROOT/MC/config/PWGEM/ini/$CONFIGNAME" + +# run the central anchor steering script; this includes +# * derive timestamp +# * derive interaction rate +# * extract and prepare configurations (which detectors are contained in the run etc.) +# * run the simulation (and QC) +# To disable QC, uncomment the following line +#export DISABLE_QC=1 +${O2DPG_ROOT}/MC/run/ANCHOR/anchorMC.sh diff --git a/MC/utils/AODBcRewriter.C b/MC/utils/AODBcRewriter.C index 710b643bd..4c0240fdb 100644 --- a/MC/utils/AODBcRewriter.C +++ b/MC/utils/AODBcRewriter.C @@ -14,8 +14,9 @@ // // This tool fixes all three problems in one pass per DF_ directory: // -// Stage 0 — Sort & deduplicate the BC table. Build BC permutation map: -// bcPerm[oldBCrow] = newBCrow. +// Stage 0 — Deduplicate the BC table in place (order-preserving; the input +// must already be globalBC-sorted). Build BC permutation map: +// bcPerm[oldBCrow] = newBCrow (monotonic non-decreasing). // // Stage 1 — Process every table that carries fIndexBCs / fIndexBC. // Remap the index via bcPerm, sort rows by the new index, and @@ -501,18 +502,31 @@ static PermMap rewriteTable(TTree *src, TDirectory *dirOut, } // ============================================================================ -// SECTION 4 — Stage 0: BC table sort + deduplication +// SECTION 4 — Stage 0: BC table deduplication (order-preserving) // ============================================================================ // -// Reads fGlobalBC from the BC tree, sorts rows, drops exact-duplicate BC -// values, and writes the compacted table. Returns bcPerm[oldRow] = newRow. +// Reads fGlobalBC from the BC tree, drops exact-duplicate BC values IN PLACE +// (preserving input row order), and writes the compacted table. Returns +// bcPerm[oldRow] = newRow. +// +// The dedup is deliberately order-preserving so that bcPerm is monotonic +// non-decreasing. This matters because every BC-indexed table (collisions, +// FT0/FV0/FDD/Zdc, ...) is sliced per BC and must stay sorted by its fIndexBCs, +// and collisions are in turn the grouping anchor for tracks (sorted by +// fIndexCollisions). A non-order-preserving BC remap would force a full reorder +// cascade BC -> collisions -> tracks to keep all those groupings valid; keeping +// bcPerm monotonic means none of those tables need to be reordered at all. +// +// This REQUIRES the input BC table to already be sorted by fGlobalBC (the +// standard AO2D invariant; also asserted on the output by validateDF check #1). +// We assert it loudly rather than silently emit a non-monotonic BC table. struct BCStage0Result { - PermMap bcPerm; // bcPerm[oldRow] = newRow in sorted/deduped BC table + PermMap bcPerm; // bcPerm[oldRow] = newRow in the deduped BC table Long64_t nUnique = 0; }; -static BCStage0Result stage0_sortBCs(TTree *treeBCs, TDirectory *dirOut) { +static BCStage0Result stage0_dedupBCs(TTree *treeBCs, TDirectory *dirOut) { BCStage0Result res; Long64_t n = treeBCs->GetEntries(); if (n == 0) return res; @@ -525,19 +539,29 @@ static BCStage0Result stage0_sortBCs(TTree *treeBCs, TDirectory *dirOut) { std::vector gbcs(n); for (Long64_t i = 0; i < n; ++i) { treeBCs->GetEntry(i); gbcs[i] = gbc; } - // Sort row indices by fGlobalBC - std::vector order(n); - std::iota(order.begin(), order.end(), 0); - std::stable_sort(order.begin(), order.end(), - [&](Long64_t a, Long64_t b){ return gbcs[a] < gbcs[b]; }); + // The BC table must already be sorted by fGlobalBC: the dedup below is + // order-preserving (merges only adjacent equal-globalBC rows), which keeps + // bcPerm monotonic and avoids a reorder cascade through collisions/tracks. + // A non-monotonic input would silently break that guarantee, so abort loudly. + for (Long64_t i = 1; i < n; ++i) { + if (gbcs[i] < gbcs[i - 1]) { + std::cerr << "FATAL: O2bc_* table is not sorted by fGlobalBC (row " << i + << " globalBC=" << gbcs[i] << " < row " << (i - 1) + << " globalBC=" << gbcs[i - 1] << ").\n" + << " AODBcRewriter requires a globalBC-sorted BC table so that\n" + << " BC deduplication is order-preserving; aborting.\n"; + std::abort(); + } + } - // Build deduplicated row list and the permutation + // Build the deduplicated row list and the (monotonic) permutation in source + // row order: adjacent rows sharing a globalBC collapse onto one output row. res.bcPerm.assign(n, -1); std::vector rowOrder; // source rows to keep, in output order - ULong64_t prev = ULong64_t(-1); + ULong64_t prev = 0; Int_t newRow = -1; - for (Long64_t srcRow : order) { - if (gbcs[srcRow] != prev) { + for (Long64_t srcRow = 0; srcRow < n; ++srcRow) { + if (newRow < 0 || gbcs[srcRow] != prev) { ++newRow; prev = gbcs[srcRow]; rowOrder.push_back(srcRow); @@ -547,7 +571,7 @@ static BCStage0Result stage0_sortBCs(TTree *treeBCs, TDirectory *dirOut) { } res.nUnique = rowOrder.size(); - std::cout << " BC stage: " << n << " rows -> " << res.nUnique << " unique\n"; + std::cout << " BC stage: " << n << " rows -> " << res.nUnique << " unique (in-place dedup)\n"; // Write the BC table (no index remapping needed for the table itself) rewriteTable(treeBCs, dirOut, rowOrder, /*indexBranch=*/"", /*parentPerm=*/{}); @@ -858,6 +882,98 @@ static const PermMap *findPermByPrefix( return nullptr; } +// ============================================================================ +// SECTION 9b — Stage 1b: Collision-grouped track tables +// ============================================================================ +// +// The primary track tables (O2track_iu, O2mfttrack_*, O2fwdtrack) are GROUPED +// by collision. O2's slicing cache (ArrowTableSlicingCache::validateOrder) +// requires every fIndexCollisions group — including the "-1" ambiguous group — +// to be a single contiguous run; otherwise it aborts with +// "Table ... index fIndexCollisions has a group with index -1 that is split". +// +// When several MC sub-timeframes are merged into one DF_ folder (data-embedding +// anchoring, which stores MC timeframes under the same DF_ as the parent data +// file), each sub-frame contributes its own [collision-grouped][-1 ambiguous] +// block. Concatenating them splits the -1 group into N runs, so the table is +// no longer sliceable. Stage 1 only reorders BC-indexed tables, and tracks are +// otherwise written in input row order, so the split survives into the output. +// +// This stage re-establishes the grouping: it reorders each collision-grouped +// track table by its remapped fIndexCollisions (stable, with -1 sinking to the +// end so the ambiguous group is one contiguous run — matching the Stage 1 +// convention) and publishes the resulting row permutation. Downstream: +// * paste-join children (O2trackextra, O2trackcov_iu, O2mctracklabel, ...) +// follow the published parent permutation; +// * every fIndexTracks* / fIndexMFTTracks / fIndexFwdTracks reference is +// remapped through it in processPasteJoinTables. +static bool isCollGroupedTrackTable(const std::string &tname) { + static const char *kPrefixes[] = {"O2track_iu", "O2track", + "O2mfttrack", "O2fwdtrack"}; + for (auto *p : kPrefixes) + if (TString(tname.c_str()).BeginsWith(p)) return true; + return false; +} + +static void stage1b_reorderTrackTables( + TDirectory *dirIn, TDirectory *dirOut, + std::unordered_map &allPerms, + std::unordered_set &written) { + + const PermMap *collPermP = findPermByPrefix(allPerms, "O2collision_"); + if (!collPermP) return; // no collisions present — nothing to regroup against + + TIter it(dirIn->GetListOfKeys()); + while (TKey *key = static_cast(it())) { + if (TString(key->GetClassName()) != "TTree") continue; + std::unique_ptr obj(key->ReadObj()); + TTree *src = dynamic_cast(obj.get()); + if (!src) continue; + + std::string tname = src->GetName(); + if (written.count(tname)) continue; // BC-indexed tracks etc. already done + if (!isCollGroupedTrackTable(tname)) continue; + if (isPasteJoinChild(tname)) continue; // children follow their parent below + if (!src->GetBranch("fIndexCollisions")) continue; + + std::cout << " Stage1b [coll-grouped]: " << tname << "\n"; + + Long64_t nSrc = src->GetEntries(); + TBranch *inIdxBr = src->GetBranch("fIndexCollisions"); + TLeaf *idxLeaf = static_cast(inIdxBr->GetListOfLeaves()->At(0)); + ScalarTag idxTag = tagOf(idxLeaf); + std::vector idxBuf(byteSize(idxTag), 0); + inIdxBr->SetAddress(idxBuf.data()); + + struct SortEntry { Long64_t newColl; Long64_t srcRow; }; + std::vector entries; + entries.reserve(nSrc); + for (Long64_t i = 0; i < nSrc; ++i) { + inIdxBr->GetEntry(i); + Long64_t oldColl = readAsInt(idxBuf.data(), idxTag); + Long64_t newColl = (oldColl >= 0 && oldColl < (Long64_t)collPermP->size()) + ? (*collPermP)[oldColl] : -1; + entries.push_back({newColl, i}); + } + // Stable-sort by remapped collision; the ambiguous group (-1) sinks to the + // end as a single contiguous run. Stable keeps the within-collision order. + std::stable_sort(entries.begin(), entries.end(), + [](const SortEntry &a, const SortEntry &b){ + if (a.newColl < 0 && b.newColl >= 0) return false; + if (a.newColl >= 0 && b.newColl < 0) return true; + return a.newColl < b.newColl; + }); + std::vector rowOrder; + rowOrder.reserve(nSrc); + for (auto &e : entries) rowOrder.push_back(e.srcRow); + + // Reorder rows and remap fIndexCollisions values through collPerm. + PermMap perm = rewriteTable(src, dirOut, rowOrder, "fIndexCollisions", *collPermP); + allPerms[tname] = std::move(perm); + written.insert(tname); + } +} + static void processPasteJoinTables( TDirectory *dirIn, TDirectory *dirOut, const std::unordered_map &allPerms, @@ -870,6 +986,12 @@ static void processPasteJoinTables( const PermMap *mcParticlePerm = findPermByPrefix(allPerms, "O2mcparticle"); const PermMap *mcCollPermP = findPermByPrefix(allPerms, "O2mccollision_"); const PermMap *collPermP = findPermByPrefix(allPerms, "O2collision_"); + // Track tables reordered in Stage 1b: every reference into them must be + // remapped through their permutation (null if the table is absent / wasn't + // reordered, in which case no remap is needed). + const PermMap *trkPerm = findPermByPrefix(allPerms, "O2track_iu"); + const PermMap *mftPerm = findPermByPrefix(allPerms, "O2mfttrack"); + const PermMap *fwdPerm = findPermByPrefix(allPerms, "O2fwdtrack"); // bcPermP is passed in from processDF (the BC table is the only stage // whose permutation isn't already published in allPerms). @@ -918,6 +1040,23 @@ static void processPasteJoinTables( extraRemaps.push_back({"fIndexBC", bcPermP}); } + // Track-pointing indices: the track tables may have been reordered in + // Stage 1b, so every reference into them must be remapped through the + // corresponding permutation. (No-op when the perm is null / absent.) + auto addTrkRemap = [&](const char *br, const PermMap *pm) { + if (pm && src->GetBranch(br)) extraRemaps.push_back({br, pm}); + }; + addTrkRemap("fIndexTracks", trkPerm); + addTrkRemap("fIndexTracks_0", trkPerm); + addTrkRemap("fIndexTracks_1", trkPerm); + addTrkRemap("fIndexTracks_2", trkPerm); + addTrkRemap("fIndexTracks_Pos", trkPerm); + addTrkRemap("fIndexTracks_Neg", trkPerm); + addTrkRemap("fIndexTracks_ITS", trkPerm); + addTrkRemap("fIndexMFTTracks", mftPerm); + addTrkRemap("fIndexFwdTracks", fwdPerm); + addTrkRemap("fIndexFwdTracks_MatchMCHTrack", fwdPerm); + // Find a paste-join parent for this table (kPasteJoins lookup). const PermMap *parentPerm = nullptr; std::string parentName; @@ -1020,10 +1159,10 @@ static void processDF(TDirectory *dirIn, TDirectory *dirOut) { return; } - // ---- Stage 0: sort & deduplicate BCs ---- + // ---- Stage 0: deduplicate BCs (order-preserving) ---- std::cout << "-- Stage 0: BCs --\n"; dirOut->cd(); - BCStage0Result s0 = stage0_sortBCs(treeBCs, dirOut); + BCStage0Result s0 = stage0_dedupBCs(treeBCs, dirOut); if (treeFlags) stage0_copyBCFlags(treeFlags, dirOut, s0.bcPerm); // Track which tree names have been written so we don't double-write @@ -1056,6 +1195,13 @@ static void processDF(TDirectory *dirIn, TDirectory *dirOut) { std::cout << " (no MCCollision table found — skipping stage 2)\n"; } + // ---- Stage 1b: regroup collision-grouped track tables ---- + // Must run after Stage 1 (needs the collision permutation) and before the + // paste-join stage (so children follow the new track order and fIndexTracks* + // references are remapped). Publishes track permutations into stage1Perms. + std::cout << "-- Stage 1b: collision-grouped track tables --\n"; + stage1b_reorderTrackTables(dirIn, dirOut, stage1Perms, written); + // ---- Paste-join tables + unrelated tables ---- std::cout << "-- Paste-join and unrelated tables --\n"; processPasteJoinTables(dirIn, dirOut, stage1Perms, written, &s0.bcPerm); diff --git a/MC/utils/o2dpg_data_embedding_utils.py b/MC/utils/o2dpg_data_embedding_utils.py index dad1769a9..b152d991e 100644 --- a/MC/utils/o2dpg_data_embedding_utils.py +++ b/MC/utils/o2dpg_data_embedding_utils.py @@ -1,6 +1,13 @@ # Set of python modules/util functions for the MC-to-DATA embedding # Mostly concerning extraction of MC collision context from existing data AO2D.root +import warnings +warnings.filterwarnings( + "ignore", + message="pandas.Int64Index is deprecated", + category=FutureWarning, +) + import ROOT import uproot import pandas as pd @@ -240,7 +247,7 @@ def fetch_bccoll_to_localFile(alien_file, local_filename): return True -def convert_to_digicontext(aod_timeframe=None, timeframeID=-1): +def convert_to_digicontext(aod_timeframe, df_folder, timeframeID=-1): """ converts AOD collision information from AO2D to collision context which can be used for MC @@ -282,13 +289,16 @@ def convert_to_digicontext(aod_timeframe=None, timeframeID=-1): # set the bunch filling ---> NEED to fetch it from CCDB # digicontext.setBunchFilling(bunchFillings[0]); + + # TODO: set the interaction rate (for TPC loopers) + # digicontext.mDigitizerInteractionRate = ... prefixes = ROOT.std.vector("std::string")(); prefixes.push_back("sgn") digicontext.setSimPrefixes(prefixes); digicontext.printCollisionSummary(); - digicontext.saveToFile(f"collission_context_{timeframeID}.root") + digicontext.saveToFile(f"collission_context_{df_folder}:{timeframeID}.root") def process_data_AO2D(file_name, run_number, upper_limit = -1): @@ -315,7 +325,8 @@ def process_data_AO2D(file_name, run_number, upper_limit = -1): break tf = row['timeframeID'] cols = row['position_vectors'] - convert_to_digicontext(cols, tf) + df = key.split('@')[0] # this is the DF folder name + convert_to_digicontext(cols, df, tf) counter = counter + 1 @@ -324,7 +335,7 @@ def main(): parser.add_argument("--run-number", type=int, help="Run number to anchor to", required=True) parser.add_argument("--aod-file", type=str, help="Data AO2D file (can be on AliEn)", required=True) - parser.add_argument("--limit", type=int, default=-1, help="Upper limit of timeframes to be extracted") + parser.add_argument("--limit", type=int, default=-1, help="Upper limit of timeframes to be extracted (-1 is no limit)") args = parser.parse_args() process_data_AO2D(args.aod_file, args.run_number, args.limit)