CbmRoot
Loading...
Searching...
No Matches
LmvmFastSim.cxx
Go to the documentation of this file.
1/* Copyright (C) 2023-2025 UGiessen, Giessen
2 SPDX-License-Identifier: GPL-3.0-only
3 Authors: Cornelius Feier-Riesen [committer] */
4
5#include "LmvmFastSim.h"
6
7#include "CbmHistManager.h"
8#include "CbmUtils.h"
9#include "LmvmDef.h"
10#include "TF1.h"
11#include "TFile.h"
12#include "TH1.h"
13#include "TH1D.h"
14#include "TH2D.h"
15#include "TLorentzVector.h"
16#include "TRandom.h"
17#include "TStopwatch.h"
18#include "TSystem.h"
19
20#include <algorithm>
21#include <iomanip>
22#include <iostream>
23#include <string>
24
25using namespace std;
26using namespace Cbm;
27
28
29void LmvmFastSim::DoFastSim(vector<string> files, const string& dataDir, const string& outFile, int nofFastSimEv,
30 const string& bgHistName)
31{
32 /*************************************************** QUICK MANUAL ****************************************************
33* HISTOGRAM NAMES *
34* The names of the FastSim histograms that are delivered from analysis to LmvmFastSim must be composed of six *
35* sections which are seperated by the underscore '_'. *
36* The 1st section must be "hfsc" or "hfsp" ("hfs" stands for 'Histogram for Fast Simulation', and "c" ("p") for *
37* charge-based (particle-based). Plutos have to be stored in "hfsc" histograms, thus charge-based in considered as *
38* default. *
39* The 2nd section classifies if the histogram is a multiplicity ("mult") or a momentum histogram ("mom"). *
40* The 3rd section classifies if it is a UrQMD ("urqmd") or a Pluto ("pluto") particle. *
41* The 4th section is the name of the particle (for charge-based: just "urqmd"; thus "urqmd" occurs twice in histname)*
42* The 5th section classifies the charge ("plus" or "minus"). Optionally, for the charged-based procedure the charge *
43* strings "plusEl" and "minusEl" can be used to indicate that the corresponding particle is a true electron *
44* (abs(pdg) == 11) (not available for particle-based procedure). *
45* The 6th section gives space for own categories, as e.g. analysis steps. *
46* *
47* Examples: *
48* > "hfsc_mult_pluto_phi_minus_elid": multiplicity histogram of an electron from Pluto-phi after electron-ID *
49* > "hfsc_mom_urqmd_urqmd_plus_gamma": momentum histogram of a positive charged UrQMD particle after gamma cut *
50* > "hfsc_mom_urqmd_urqmd_plusEl_gamma": momentum histogram of a UrQMD positron (true electron!) after gamma cut *
51* > "hfsp_mom_urqmd_etaEl_plus_elid": momentum histogram of a UrQMD eta-positron (true electron!) after elid cut*
52* *
53* *
54* NUMBER OF DIFFERENT PLUTO PARTICLES *
55* If files.size() > 1, it is assumed that only one Pluto particle is embedded per event and each file corresponds to *
56* another Pluto particle. In this case the number of Pluto particles corresponds to the number of different files, *
57* i.e. files.size(). If files.size() == 1, it is assumed that per event of each Pluto particle exactly one is *
58* embedded. In this case, the number of different Pluto particles will be extracted by the number of different names *
59* in the 4th section that are not "urqmd". For this, make sure the correctness of the names! *
60* *
61* CreateMeanHists() *
62* Mean histograms are histograms that are merged of UrQMD histograms from all files with different Pluto signals. *
63* Thus, if files.size() = 1, the mean histogram is simply the single UrQMD histogram of that file. *
64* *
65* -> SEE MANUAL FOR MORE DETAILED DESCRIPTION *
66*********************************************************************************************************************/
67
68 fOutputDir = dataDir;
69 fOutFileName = outFile;
70 fNofFastSimEv = nofFastSimEv;
71 fBgHistName = bgHistName;
72 TFile* oldFile = gFile;
73 TDirectory* oldDir = gDirectory;
74
75 // Load fH[iFile]
76 fH.resize(files.size());
77 for (size_t iF = 0; iF < fH.size(); iF++) {
78 fH[iF] = new LmvmHist();
79 TFile* file = new TFile(files[iF].c_str()); // FIXME: file is not closed; needs to remain until procedure ends
80 fH[iF]->fHM.ReadFromFile(file);
81 int nofEventsSig = (int) fH[iF]->H1("hEventNumber")->GetEntries();
82 fNofSimEv += nofEventsSig;
83 LOG(info) << "Signal:" << fHMean.fSignalNames[iF] << " nofEvents:" << nofEventsSig << endl;
84 }
85 LOG(info) << "Nof simulated Events = " << fNofSimEv;
86 LOG(info) << "Nof FastSim Events per Job = " << fNofFastSimEv / 1e6 << " Mio" << endl;
87 LOG(info) << "fMixingDepth = " << fMixingDepth;
88 if (fNofFastSimEv < fMixingDepth)
89 LOG(warn) << "LmvmFastSim::DoFastSim: 'fNofFastSimEv' is smaller than 'fMixingDepth'. No histograms for mixed "
90 "events will be filled.";
91
92 // Get Pluto Mode: one or all Plutos per event?
93 if (files.size() == 1)
94 fPlutoMode = "all"; // per event, exactly one of every Pluto particle species is implemented
95 else if (files.size() > 1)
96 fPlutoMode = "one"; // only one Pluto particle per event, i.e., seperate files per Pluto
97 else
98 LOG(fatal) << "LmvmFastSim::DoFastSim(): 'vector<string> files' has no content!";
99 LOG(info) << "Pluto mode: '" << fPlutoMode << "'";
100
101 CheckHistogramNames();
102 GetChargeCategoryVector();
103 GetIndividualCategories();
104 GetPlutoAndUrqmdNames(files);
105 CreateHistos();
106 fHMean.FillH1("hfs_nofJobs", 0.5);
107
108 // Declare vectors that will be filled with either 'plus' or 'minus' particles
109 // Charge-based
110 vector<LmvmDataFastSim> C_Plus, C_Minus; // standard with all (PLUTO + UrQMD)
111 vector<LmvmDataFastSim> C_PlusEl, C_MinusEl; // electrons only (PLUTO + UrQMD)
112 vector<LmvmDataFastSim> C_PlusU, C_MinusU; // only UrQMD
113 vector<LmvmDataFastSim> C_PlusUEl, C_MinusUEl; // only UrQMD electrons
114 vector<LmvmDataFastSim> C_PlusMix, C_MinusMix; // mixed events
115
116 // Particle-based
117 vector<LmvmDataFastSim> P_Plus, P_Minus;
118 vector<LmvmDataFastSim> P_PlusEl, P_MinusEl;
119 vector<LmvmDataFastSim> P_PlusMix, P_MinusMix;
120
121 // Generate fast simulation events
122 LOG(info) << "Generating Fast Simulation Background ...";
123 TStopwatch timer;
124 timer.Start();
125 gRandom->SetSeed(0);
126
127 for (int iCat = 0; iCat < fNofCats; iCat++) {
128 C_PlusMix.clear();
129 C_MinusMix.clear();
130 P_PlusMix.clear();
131 P_MinusMix.clear();
132 for (int iN = 1; iN <= fNofFastSimEv; iN++) {
133 if (iN % (int) 1e6 == 0) LOG(info) << "Category '" << fCatNames[iCat] << "': iN = " << iN / 1e6 << " Mio ...";
134 if (iCat == 0) fHMean.FillH1("hfs_EventNumber", 0.5); // has to be filled only for one category
135 C_Plus.clear();
136 C_Minus.clear();
137 C_PlusEl.clear();
138 C_MinusEl.clear();
139 C_PlusU.clear();
140 C_MinusU.clear();
141 C_PlusUEl.clear();
142 C_MinusUEl.clear();
143 P_Plus.clear();
144 P_Minus.clear();
145 P_PlusEl.clear();
146 P_MinusEl.clear();
147
148 // Filling LmvmDataFastSim vectors
149
150 // Plutos
151 for (int iPluto = 0; iPluto < fNofPlutos; iPluto++) {
152 if (fPlutoMode == "one" && iN % fNofPlutos != 0) break; // only one Pluto per event
153 for (const string charge : {"plus", "minus"}) {
154 string hMultNameCorr = "hfsc_multCorr_pluto_" + fPlutoNames[iPluto] + "_" + charge + "_" + fCatNames[iCat];
155 string hMomName = "hfsc_mom_pluto_" + fPlutoNames[iPluto] + "_" + charge + "_" + fCatNames[iCat];
156 string hMultNameRndm = "hfsc_multRndm_pluto_" + fPlutoNames[iPluto] + "_" + charge + "_" + fCatNames[iCat];
157 string hMomNameRndm = "hfsc_momRndm_pluto_" + fPlutoNames[iPluto] + "_" + charge + "_" + fCatNames[iCat];
158 int mult = (int) fHMean.H1(hMultNameCorr)->GetRandom();
159 fHMean.H1(hMultNameRndm)->Fill(mult);
160 if (mult < 1) continue;
161 for (int iMult = 1; iMult <= mult; iMult++) {
162 double fastsimMom[3] = {0.};
163 fH[iPluto]->fHM.HnSparse(hMomName)->GetRandom(fastsimMom);
164 fHMean.fHM.HnSparse(hMomNameRndm)->Fill(fastsimMom);
165 int iEvent = (fPlutoMode == "one") ? iN - iPluto : iN;
166 LmvmDataFastSim mom(fastsimMom[0], fastsimMom[1], fastsimMom[2], iEvent, fPlutoNames[iPluto]);
167 if (charge == "plus") {
168 C_Plus.push_back(mom);
169 C_PlusMix.push_back(mom);
170 if (fDoParticleBased) {
171 P_Plus.push_back(mom);
172 P_PlusMix.push_back(mom);
173 }
174 if (fDoTrueEl) {
175 C_PlusEl.push_back(mom);
176 if (fDoParticleBased) P_PlusEl.push_back(mom);
177 }
178 }
179 else if (charge == "minus") {
180 C_Minus.push_back(mom);
181 C_MinusMix.push_back(mom);
182 if (fDoParticleBased) {
183 P_Minus.push_back(mom);
184 P_MinusMix.push_back(mom);
185 }
186 if (fDoTrueEl) {
187 C_MinusEl.push_back(mom);
188 if (fDoParticleBased) P_MinusEl.push_back(mom);
189 }
190 }
191 }
192 }
193 }
194
195 // UrQMD
196 // Charge-based
197 for (const string& charge : fChargeCats) {
198 string hMultName = "hfsc_mult_urqmd_urqmd_" + charge + "_" + fCatNames[iCat];
199 string hMomName = "hfsc_mom_urqmd_urqmd_" + charge + "_" + fCatNames[iCat];
200 string hMultNameRndm = "hfsc_multRndm_urqmd_urqmd_" + charge + "_" + fCatNames[iCat];
201 string hMomNameRndm = "hfsc_momRndm_urqmd_urqmd_" + charge + "_" + fCatNames[iCat];
202 int mult = (int) fHMean.H1(hMultName)->GetRandom();
203 fHMean.H1(hMultNameRndm)->Fill(mult);
204 if (mult < 1) continue;
205 for (int iMult = 1; iMult <= mult; iMult++) {
206 double fastsimMom[3] = {0.};
207 fHMean.fHM.HnSparse(hMomName)->GetRandom(fastsimMom);
208 fHMean.fHM.HnSparse(hMomNameRndm)->Fill(fastsimMom);
209 LmvmDataFastSim mom(fastsimMom[0], fastsimMom[1], fastsimMom[2], iN, "");
210 if (charge == "plus") {
211 C_Plus.push_back(mom);
212 C_PlusMix.push_back(mom);
213 C_PlusU.push_back(mom);
214 }
215 else if (charge == "minus") {
216 C_Minus.push_back(mom);
217 C_MinusMix.push_back(mom);
218 C_MinusU.push_back(mom);
219 }
220 else if (charge == "plusEl") {
221 C_PlusEl.push_back(mom);
222 C_PlusUEl.push_back(mom);
223 }
224 else if (charge == "minusEl") {
225 C_MinusEl.push_back(mom);
226 C_MinusUEl.push_back(mom);
227 }
228 }
229 } // ..charge
230 CalcMinvAndFillHist(C_Plus, C_Minus, "hfsc_MinvBg_" + fCatNames[iCat]);
231 CalcMinvAndFillHist(C_PlusU, C_MinusU, "hfsc_MinvBg_urqmdAll_" + fCatNames[iCat]);
232 CalcMinvAndFillHist(C_Plus, C_Minus, "hfsc_MinvPM_sameEv_" + fCatNames[iCat]);
233 CalcMinvAndFillHist(C_Plus, C_Plus, "hfsc_MinvPP_sameEv_" + fCatNames[iCat], "likesign");
234 CalcMinvAndFillHist(C_Minus, C_Minus, "hfsc_MinvMM_sameEv_" + fCatNames[iCat], "likesign");
235 if (fDoTrueEl) {
236 CalcMinvAndFillHist(C_PlusEl, C_MinusEl, "hfsc_MinvBg_trueEl_" + fCatNames[iCat]);
237 CalcMinvAndFillHist(C_PlusUEl, C_MinusUEl, "hfsc_MinvBg_urqmdEl_" + fCatNames[iCat]);
238 }
239 if (fDoMixEvents && iN % fMixingDepth == 0) { // Fill histograms for mixed event data
240 CalcMinvAndFillHist(C_PlusMix, C_MinusMix, "hfsc_MinvPM_mixedEv_" + fCatNames[iCat], "mix");
241 CalcMinvAndFillHist(C_PlusMix, C_PlusMix, "hfsc_MinvPP_mixedEv_" + fCatNames[iCat], "mix likesign");
242 CalcMinvAndFillHist(C_MinusMix, C_MinusMix, "hfsc_MinvMM_mixedEv_" + fCatNames[iCat], "mix likesign");
243 C_PlusMix.clear();
244 C_MinusMix.clear();
245 }
246
247 // Particle-based
248 if (fDoParticleBased) {
249 for (const string& ptcl : fUrqmdNames) {
250 for (const string charge : {"plus", "minus"}) {
251 string hMultName = "hfsp_mult_urqmd_" + ptcl + "_" + charge + "_" + fCatNames[iCat];
252 string hMultNameRndm = "hfsp_multRndm_urqmd_" + ptcl + "_" + charge + "_" + fCatNames[iCat];
253 string hMomName = "hfsp_mom_urqmd_" + ptcl + "_" + charge + "_" + fCatNames[iCat];
254 string hMomNameRndm = "hfsp_momRndm_urqmd_" + ptcl + "_" + charge + "_" + fCatNames[iCat];
255 int mult =
256 fHMean.H1(hMultName)
257 ->GetRandom(); //TODO urgent: check why this gives error? Still present after updating hist naming?
258 fHMean.H1(hMultNameRndm)->Fill(mult);
259 if (mult < 1) continue;
260 for (int iMult = 1; iMult <= mult; iMult++) {
261 double fastsimMom[3] = {0.};
262 fHMean.fHM.HnSparse(hMomName)->GetRandom(fastsimMom);
263 fHMean.fHM.HnSparse(hMomNameRndm)->Fill(fastsimMom);
264 LmvmDataFastSim mom(fastsimMom[0], fastsimMom[1], fastsimMom[2], iN, "");
265 if (charge == "plus") {
266 P_Plus.push_back(mom);
267 P_PlusMix.push_back(mom);
268 }
269 else if (charge == "minus") {
270 P_Minus.push_back(mom);
271 P_MinusMix.push_back(mom);
272 }
273 if (ptcl.find("El") != string::npos) {
274 if (charge == "plus")
275 P_PlusEl.push_back(mom);
276 else if (charge == "minus")
277 P_MinusEl.push_back(mom);
278 }
279 }
280 }
281 } // ..ptcl
282 CalcMinvAndFillHist(P_Plus, P_Minus, "hfsp_MinvBg_" + fCatNames[iCat]);
283 CalcMinvAndFillHist(P_PlusEl, P_MinusEl, "hfsp_MinvBg_trueEl_" + fCatNames[iCat]);
284 CalcMinvAndFillHist(P_Plus, P_Minus, "hfsp_MinvPM_sameEv_" + fCatNames[iCat]);
285 CalcMinvAndFillHist(P_Plus, P_Plus, "hfsp_MinvPP_sameEv_" + fCatNames[iCat], "likesign");
286 CalcMinvAndFillHist(P_Minus, P_Minus, "hfsp_MinvMM_sameEv_" + fCatNames[iCat], "likesign");
287 if (fDoMixEvents && iN % fMixingDepth == 0) { // Fill histograms for mixed event data
288 CalcMinvAndFillHist(P_PlusMix, P_MinusMix, "hfsp_MinvPM_mixedEv_" + fCatNames[iCat], "mix");
289 CalcMinvAndFillHist(P_PlusMix, P_PlusMix, "hfsp_MinvPP_mixedEv_" + fCatNames[iCat], "mix likesign");
290 CalcMinvAndFillHist(P_MinusMix, P_MinusMix, "hfsp_MinvMM_mixedEv_" + fCatNames[iCat], "mix likesign");
291 P_PlusMix.clear();
292 P_MinusMix.clear();
293 }
294 } // ..fDoParticleBased
295 } // ..iN
296 } // ..iCat
297
298 double bW = fHMean.H1("hfsc_MinvBg_" + fCatNames[0])->GetBinWidth(1);
299 fHMean.fHM.ScaleByPattern("hfs.*_Minv.*", 1. / bW);
300
301 LOG(info) << "[DONE]." << endl;
302
303 timer.Stop();
304 fHMean.FillH1("hfs_time", 1e6 * timer.RealTime() / (fNofFastSimEv * fNofCats));
305 LOG(info) << "Real time " << timer.RealTime() << " s, CPU time " << timer.CpuTime() << " s";
306 LOG(info) << "Time per event (average): " << 1e6 * timer.RealTime() / (fNofFastSimEv * fNofCats) << " mus";
307
308 static ProcInfo_t info;
309 const float toMB = 1.f / 1024.f;
310 gSystem->GetProcInfo(&info);
311 fHMean.FillH1("hfs_memory_res", info.fMemResident * toMB);
312 fHMean.FillH1("hfs_memory_vir", info.fMemVirtual * toMB);
313
314 SaveHist();
315 CheckMemory("Memory after 'SaveHists'");
316
317 // delete fH[] elements
318 for (LmvmHist* fHElement : fH) {
319 delete fHElement;
320 }
321
322 CheckMemory("Memory after deleting fH[] elements");
323
324 gFile = oldFile;
325 gDirectory = oldDir;
326}
327
328void LmvmFastSim::CalcMinvAndFillHist(const vector<LmvmDataFastSim>& Cand1, const vector<LmvmDataFastSim>& Cand2,
329 const string& hName)
330{
331 CalcMinvAndFillHist(Cand1, Cand2, hName, "");
332}
333
334void LmvmFastSim::CalcMinvAndFillHist(const vector<LmvmDataFastSim>& Cand1, const vector<LmvmDataFastSim>& Cand2,
335 const string& hName, const string& opt)
336{
337 if (opt != "" && opt.find("mix") == string::npos && opt.find("likesign") == string::npos)
338 LOG(error) << "CalcMinvAndFillHist: Option " << opt << " is not available!";
339 for (size_t iCand1 = 0; iCand1 < Cand1.size(); iCand1++) {
340 for (size_t iCand2 = 0; iCand2 < Cand2.size(); iCand2++) {
341 bool isSameEvent = (Cand1[iCand1].fEventNumber == Cand2[iCand2].fEventNumber);
342 if (opt.find("mix") != string::npos && isSameEvent) continue; // only mixed events
343 if (opt.find("likesign") != string::npos && iCand2 <= iCand1) continue; // no double counting of like-sign pairs
344
345 //As in analysis, don't pair two electrons from different Plutos (check if you really want to skip these pairs!)
346 bool isPlutoPair = (Cand1[iCand1].fPluto != "" && Cand2[iCand2].fPluto != "");
347 bool isTruePlutoPair = (isPlutoPair && isSameEvent && Cand1[iCand1].fPluto == Cand2[iCand2].fPluto);
348 if (fPlutoMode == "one" && isPlutoPair && !isTruePlutoPair) continue;
349
350 TVector3 momP(Cand1[iCand1].fPx, Cand1[iCand1].fPy, Cand1[iCand1].fPz);
351 TVector3 momM(Cand2[iCand2].fPx, Cand2[iCand2].fPy, Cand2[iCand2].fPz);
352 double angle = momP.Angle(momM);
353 double minv = 2. * TMath::Sin(angle / 2.) * TMath::Sqrt(momP.Mag() * momM.Mag());
354 fHMean.H1(hName)->Fill(minv);
355 }
356 }
357}
358
359void LmvmFastSim::CheckHistogramNames() // Check if handed-over histogram names fulfill the conventions and set mode
360{
361 // It is assumed that the histo structure of all input files is the same, so only first file is checked.
362 vector<TH1*> hnames = fH[0]->fHM.H1Vector("hfsc_.*"); // TODO: |hfsp)
363 vector<TH1*> hmult = fH[0]->fHM.H1Vector("hfsc_mult_.*"); // TODO: |hfsp)
364 //vector<THnSparse*> hmom = fH[0]->fHM.HnSparseVector("(hfsc|hfsp)_mom_.*"); //TODO: has to be implemented in CbmHistManager
365 if (hnames.size() < 1) LOG(fatal) << "CheckHistogramNames: No histograms found for Fast Simulation procedure!";
366 if (hmult.size() < 1)
367 LOG(fatal) << "CheckHistogramNames: No multiplicity histograms found for Fast Simulation procedure!";
368 //if (hmom.size() < 1) LOG(fatal) << "No momentum histograms found for Fast Simulation procedure!";
369 for (size_t iH = 0; iH < hnames.size(); iH++) {
370 LOG(info) << "CheckHistogramNames(): HistName = " << hnames[iH]->GetName();
371 vector<string> split = Split(hnames[iH]->GetName(), '_');
372 if (split.size() != 6)
373 LOG(fatal) << "CheckHistogramNames: Histogram name '" << hnames[iH]
374 << "' is not valid. Check manual for name conventions!";
375 if (split[1] != "mult" && split[1] != "mom")
376 LOG(fatal) << "CheckHistogramNames: Histogram name section '" << split[1]
377 << "' (2nd section) is not valid. Check manual for name conventions!";
378 if (split[2] != "pluto" && split[2] != "urqmd")
379 LOG(fatal) << "CheckHistogramNames: Histogram name section '" << split[2]
380 << "' (3rd section) is not valid. Check manual for name conventions!";
381 if (split[4] != "plus" && split[4] != "minus" && split[4] != "plusEl" && split[4] != "minusEl")
382 LOG(fatal) << "CheckHistogramNames: Histogram name section '" << split[4]
383 << "' (5th section) is not valid. Check manual for name conventions!";
384 }
385 cout << endl;
386
387 vector<TH1*> hnamesC = fH[0]->fHM.H1Vector("hfsc_.*");
388 vector<TH1*> hnamesP = fH[0]->fHM.H1Vector("hfsp_.*");
389 if (hnamesC.size() >= 1) fDoChargeBased = true;
390 if (hnamesP.size() >= 1)
391 fDoParticleBased =
392 true; // TODO: some histograms are empty because of rare occurence of particles; hence 'GetRandom' (which calls TH1::ComputeIntegral) throws out error; add corresponding request before calling GetRandom in that cased!
393 if (!fDoChargeBased && !fDoParticleBased)
394 LOG(fatal)
395 << "CheckHistogramNames: No histograms found for Fast Simulation procedure!"; // TODO: redundant, should actually been covered by previous line (see top)
396 if (fDoChargeBased == true && fDoParticleBased == false)
397 LOG(info) << "Running charge-based mode only.";
398 else if (fDoChargeBased == false && fDoParticleBased == true)
399 LOG(info) << "Running particle-based mode only.";
400 else if (fDoChargeBased == true && fDoParticleBased == true)
401 LOG(info) << "Running charge-based and particle-based mode.";
402}
403
404void LmvmFastSim::GetChargeCategoryVector()
405{
406 fChargeCats = {"plus", "minus"};
407 vector<TH1*> hTrueElP = fH[0]->fHM.H1Vector("hfsc_mult_.*");
408 vector<TH1*> hTrueElM = fH[0]->fHM.H1Vector("hfsc_mult_.*_minusEl_.*");
409
410 if (hTrueElP.size() > 0 && hTrueElM.size() > 0) {
411 fDoTrueEl = true;
412 fChargeCats.push_back("plusEl");
413 fChargeCats.push_back("minusEl");
414 }
415 if ((hTrueElP.size() > 0 && hTrueElM.size() == 0) || (hTrueElP.size() == 0 && hTrueElM.size() > 0)) {
416 LOG(error) << "LmvmFastSim::GetChargeCategoryVector: Either only 'plusEl' or only 'minusEl' histos exist. Can not "
417 "procede on true electrons.";
418 }
419}
420
421void LmvmFastSim::GetIndividualCategories()
422{
423 // It is assumed that the individual categories are the same for charge-based and particle-based modes
424 string hname = (fDoChargeBased) ? "hfsc" : "hfsp";
425 vector<TH1*> hcats = fH[0]->fHM.H1Vector(hname + "_mult_urqmd_.*_plus_.*");
426 for (size_t iH = 0; iH < hcats.size(); iH++) {
427 vector<string> split = Split(hcats[iH]->GetName(), '_');
428 if (!ContainsString(fCatNames, split[5])) fCatNames.push_back(split[5]);
429 }
430 fNofCats = fCatNames.size();
431 LOG(info) << "Number of Categories: " << fNofCats;
432 cout << "[INFO] Category Names: ";
433 for (size_t iC = 0; iC < fCatNames.size(); iC++) {
434 cout << fCatNames[iC] << " ";
435 }
436 cout << endl << endl;
437}
438
439void LmvmFastSim::GetPlutoAndUrqmdNames(const vector<string>& files)
440{
441 // Pluto Names
442 if (fPlutoMode == "one") {
443 fNofPlutos = files.size();
444 for (int iPluto = 0; iPluto < fNofPlutos; iPluto++) {
445 vector<TH1*> histos = fH[iPluto]->fHM.H1Vector("hfsc_mult_pluto_.*_plus_" + fCatNames[0]);
446 if (histos.size() > 1) LOG(error) << "More than one Pluto particle included. Check!";
447 for (size_t iH = 0; iH < histos.size(); iH++) {
448 vector<string> split = Split(histos[iH]->GetName(), '_');
449 fPlutoNames.push_back(split[3]);
450 }
451 }
452 }
453 else {
454 vector<TH1*> histos = fH[0]->fHM.H1Vector("hfsc_mult_pluto_.*_plus_" + fCatNames[0]);
455 fNofPlutos = histos.size();
456 for (size_t iH = 0; iH < histos.size(); iH++) {
457 vector<string> split = Split(histos[iH]->GetName(), '_');
458 if (split[2] == "pluto") fPlutoNames.push_back(split[3]);
459 }
460 }
461 LOG(info) << "Number of Pluto Particles: " << fNofPlutos;
462 cout << "[INFO] Pluto Names: ";
463 for (size_t iPluto = 0; iPluto < fPlutoNames.size(); iPluto++) {
464 cout << fPlutoNames[iPluto] << " ";
465 }
466 cout << endl << endl;
467
468 // UrQMD Names (in case of particle-based mode)
469 vector<TH1*> uHistos = fH[0]->fHM.H1Vector("hfsp_mult_urqmd_.*_plus_" + fCatNames[0]);
470 for (size_t iH = 0; iH < uHistos.size(); iH++) {
471 vector<string> split = Split(uHistos[iH]->GetName(), '_');
472 //if (!ContainsString(fUrqmdNames, split[3]) && split[3] != "urqmd") fUrqmdNames.push_back(split[3]);
473 if (!ContainsString(fUrqmdNames, split[3])) fUrqmdNames.push_back(split[3]);
474 }
475 fNofUrqmds = fUrqmdNames.size();
476 LOG(info) << "Number of UrQMD Particles (particle-based): " << fNofUrqmds;
477 cout << "[INFO] UrQMD Names: ";
478 for (size_t iUrqmd = 0; iUrqmd < fUrqmdNames.size(); iUrqmd++) {
479 cout << fUrqmdNames[iUrqmd] << " ";
480 }
481 cout << endl << endl;
482}
483
484bool LmvmFastSim::ContainsString(const vector<string>& vName, const string& name)
485{
486 return std::find(vName.begin(), vName.end(), name) != vName.end();
487}
488
489void LmvmFastSim::CreateHistos()
490{
491 // Number of Fast Simulation Events
492 fHMean.CreateH1("hfs_EventNumber", "", "", 1, 0, 1.);
493 fHMean.CreateH1("hfs_nofJobs", "", "", 1, 0, 1.); // for correct scaling in draw class
494
495 CreateMeanHist<TH1D>(fBgHistName); // Template for background histograms (copy from bg histogram in analysis)
496
497 // Create blank background histos that will be filled in this procedure and returned
498 LOG(info) << "LmvmFastSim::CreateHistos(): fBgHistName = " << fBgHistName;
499 for (int iCat = 0; iCat < fNofCats; iCat++) {
500 // Charge-based
501 TH1D* hcBg = fHMean.CreateHByClone<TH1D>(fBgHistName,
502 "hfsc_MinvBg_" + fCatNames[iCat]); // standard histogram that is returned
503 TH1D* hcBgU = fHMean.CreateHByClone<TH1D>(fBgHistName, "hfsc_MinvBg_urqmdAll_" + fCatNames[iCat]); // only UrQMD
504 TH1D* hcBgPM = fHMean.CreateHByClone<TH1D>(fBgHistName, "hfsc_MinvPM_sameEv_"
505 + fCatNames[iCat]); // same as MinvBg ->TODO: unify this!
506 TH1D* hcBgPP = fHMean.CreateHByClone<TH1D>(fBgHistName, "hfsc_MinvPP_sameEv_" + fCatNames[iCat]);
507 TH1D* hcBgMM = fHMean.CreateHByClone<TH1D>(fBgHistName, "hfsc_MinvMM_sameEv_" + fCatNames[iCat]);
508 TH1D* hcBgpm = fHMean.CreateHByClone<TH1D>(fBgHistName, "hfsc_MinvPM_mixedEv_" + fCatNames[iCat]);
509 TH1D* hcBgpp = fHMean.CreateHByClone<TH1D>(fBgHistName, "hfsc_MinvPP_mixedEv_" + fCatNames[iCat]);
510 TH1D* hcBgmm = fHMean.CreateHByClone<TH1D>(fBgHistName, "hfsc_MinvMM_mixedEv_" + fCatNames[iCat]);
511 hcBg->Reset();
512 hcBgU->Reset();
513 hcBgPM->Reset();
514 hcBgPP->Reset();
515 hcBgMM->Reset();
516 hcBgpm->Reset();
517 hcBgpp->Reset();
518 hcBgmm->Reset();
519 if (fDoTrueEl) {
520 TH1D* hcBgEl = fHMean.CreateHByClone<TH1D>(
521 fBgHistName, "hfsc_MinvBg_trueEl_" + fCatNames[iCat]); // histogram with only true electrons
522 TH1D* hcBgUEl =
523 fHMean.CreateHByClone<TH1D>(fBgHistName, "hfsc_MinvBg_urqmdEl_" + fCatNames[iCat]); // only UrQMD electrons
524 hcBgEl->Reset();
525 hcBgUEl->Reset();
526 }
527
528 // Particle-based
529 if (fDoParticleBased) {
530 TH1D* hpBg = fHMean.CreateHByClone<TH1D>(fBgHistName, "hfsp_MinvBg_" + fCatNames[iCat]);
531 TH1D* hpBgEl = fHMean.CreateHByClone<TH1D>(fBgHistName, "hfsp_MinvBg_trueEl_" + fCatNames[iCat]);
532 TH1D* hpBgPM = fHMean.CreateHByClone<TH1D>(
533 fBgHistName, "hfsp_MinvPM_sameEv_" + fCatNames[iCat]); // same as MinvBg ->TODO: unify this!
534 TH1D* hpBgPP = fHMean.CreateHByClone<TH1D>(fBgHistName, "hfsp_MinvPP_sameEv_" + fCatNames[iCat]);
535 TH1D* hpBgMM = fHMean.CreateHByClone<TH1D>(fBgHistName, "hfsp_MinvMM_sameEv_" + fCatNames[iCat]);
536 TH1D* hpBgpm = fHMean.CreateHByClone<TH1D>(fBgHistName, "hfsp_MinvPM_mixedEv_" + fCatNames[iCat]);
537 TH1D* hpBgpp = fHMean.CreateHByClone<TH1D>(fBgHistName, "hfsp_MinvPP_mixedEv_" + fCatNames[iCat]);
538 TH1D* hpBgmm = fHMean.CreateHByClone<TH1D>(fBgHistName, "hfsp_MinvMM_mixedEv_" + fCatNames[iCat]);
539 hpBg->Reset();
540 hpBgEl->Reset();
541 hpBgPM->Reset();
542 hpBgPP->Reset();
543 hpBgMM->Reset();
544 hpBgpm->Reset();
545 hpBgpp->Reset();
546 hpBgmm->Reset();
547 }
548 }
549
550 // Create multiplicity and momentum histograms for UrQMD particles
551 // Charge-based
552 for (const string& charge : fChargeCats) {
553 for (int iCat = 0; iCat < fNofCats; iCat++) {
554 string hNameMult = "hfsc_mult_urqmd_urqmd_" + charge + "_" + fCatNames[iCat];
555 string hNameMom = "hfsc_mom_urqmd_urqmd_" + charge + "_" + fCatNames[iCat];
556 string hNameMultRndm = "hfsc_multRndm_urqmd_urqmd_" + charge + "_" + fCatNames[iCat];
557 string hNameMomRndm = "hfsc_momRndm_urqmd_urqmd_" + charge + "_" + fCatNames[iCat];
558
559 // Create mean hists from analysis file
560 CreateMeanHist<TH1D>(hNameMult);
561 CreateMeanHist<THnSparseD>(hNameMom);
562
563 // Create blank mean histos to be filled with random numbers (to check with true values)
564 TH1D* hMultRndm = fHMean.CreateHByClone<TH1D>(hNameMult, hNameMultRndm);
565 THnSparseD* hMomRndm = fHMean.CreateHByClone<THnSparseD>(hNameMom, hNameMomRndm);
566 hMultRndm->Reset();
567 hMomRndm->Reset();
568 }
569 }
570 // Particle-based
571 if (fDoParticleBased) {
572 for (const string& ptcl : fUrqmdNames) {
573 for (const string charge : {"plus", "minus"}) {
574 for (int iCat = 0; iCat < fNofCats; iCat++) {
575 string hNameMult = "hfsp_mult_urqmd_" + ptcl + "_" + charge + "_" + fCatNames[iCat];
576 string hNameMom = "hfsp_mom_urqmd_" + ptcl + "_" + charge + "_" + fCatNames[iCat];
577 string hNameMultRndm = "hfsp_multRndm_urqmd_" + ptcl + "_" + charge + "_" + fCatNames[iCat];
578 string hNameMomRndm = "hfsp_momRndm_urqmd_" + ptcl + "_" + charge + "_" + fCatNames[iCat];
579
580 // Create mean hists from analysis file
581 CreateMeanHist<TH1D>(hNameMult);
582 CreateMeanHist<THnSparseD>(hNameMom);
583
584 // Create blank mean histos to be filled with random numbers (to check with true values)
585 TH1D* hMultRndm = fHMean.CreateHByClone<TH1D>(hNameMult, hNameMultRndm);
586 THnSparseD* hMomRndm = fHMean.CreateHByClone<THnSparseD>(hNameMom, hNameMomRndm);
587 hMultRndm->Reset();
588 hMomRndm->Reset();
589 }
590 }
591 }
592 }
593
594 // Create multiplicity and momentum histos for Pluto particles (original, random, corrected bin='0/Ev')
595 for (int iPluto = 0; iPluto < fNofPlutos; iPluto++) {
596 int nofEventsSig = (int) fH[iPluto]->H1("hEventNumber")->GetEntries();
597 for (const string charge : {"plus", "minus"}) {
598 for (int iCat = 0; iCat < fNofCats; iCat++) {
599
600 // Multiplicity
601 string hNameMult = "hfsc_mult_pluto_" + fPlutoNames[iPluto] + "_" + charge + "_" + fCatNames[iCat];
602 string hNameMultCorr = "hfsc_multCorr_pluto_" + fPlutoNames[iPluto] + "_" + charge + "_" + fCatNames[iCat];
603 string hNameMultRndm = "hfsc_multRndm_pluto_" + fPlutoNames[iPluto] + "_" + charge + "_" + fCatNames[iCat];
604
605 fHMean.fHM.Add(hNameMult, (TH1D*) fH[iPluto]->H1Clone(hNameMult));
606
607 TH1D* hMultRndm = fH[iPluto]->H1Clone(hNameMult);
608 hMultRndm->Reset();
609 hMultRndm->SetNameTitle(hNameMultRndm.c_str(), hNameMultRndm.c_str());
610 fHMean.fHM.Add(hNameMultRndm, hMultRndm);
611
612 // The multiplicity yield for the bin '0/Event' are corrected for Pluto particles, since they are feeded
613 // with a weighting factor into the yield (integral over multiplicity histos should return nofEvents).
614 TH1D* hMultCorr = fH[iPluto]->H1Clone(hNameMult);
615 hMultCorr->SetNameTitle(hNameMultCorr.c_str(), hNameMultCorr.c_str());
616 double c1 = hMultCorr->GetBinContent(2); // Bin content of bin '1/Event'
617 double c0 = nofEventsSig - c1;
618 hMultCorr->SetBinContent(1, c0);
619 fHMean.fHM.Add(hNameMultCorr, hMultCorr);
620
621 // Momentum
622 string hNameMom = "hfsc_mom_pluto_" + fPlutoNames[iPluto] + "_" + charge + "_" + fCatNames[iCat];
623 string hNameMomRndm = "hfsc_momRndm_pluto_" + fPlutoNames[iPluto] + "_" + charge + "_" + fCatNames[iCat];
624
625 THnSparseD* hMomRndm = (THnSparseD*) fH[iPluto]->fHM.GetObject(hNameMom)->Clone();
626 hMomRndm->SetNameTitle(hNameMomRndm.c_str(), hNameMomRndm.c_str());
627 hMomRndm->Reset();
628 fHMean.fHM.Add(hNameMomRndm, hMomRndm);
629 } // categories
630 } // charge
631 } // Plutos
632
633 // Time per event
634 fHMean.CreateH1("hfs_time", "Time per Event [#mus]", "a.u.", 100, 0., 100.);
635
636 // Occupied memory
637 fHMean.CreateH1("hfs_memory_vir", "Memory [MB]", "a.u.", 100, 0., 6000.);
638 fHMean.CreateH1("hfs_memory_res", "Memory [MB]", "a.u.", 100, 0., 6000.);
639}
640
641template<class T>
642void LmvmFastSim::CreateMeanHist(const string& name)
643{
644 if (fPlutoMode == "one") {
645 for (int iFile = 0; iFile < fNofPlutos; iFile++) {
646 if (iFile == 0)
647 fHMean.fHM.Add(name, static_cast<T*>(fH[iFile]->GetObject(name)->Clone()));
648 else
649 static_cast<T*>(fHMean.GetObject(name))->Add(static_cast<T*>(fH[iFile]->GetObject(name)->Clone()));
650 }
651 }
652 else
653 fHMean.fHM.Add(name, static_cast<T*>(fH[0]->GetObject(name)->Clone()));
654}
655
656void LmvmFastSim::SaveHist()
657{
658 if (fOutputDir != "") {
659 gSystem->mkdir(fOutputDir.c_str(), true);
660 TFile* f = TFile::Open(string(fOutFileName).c_str(), "RECREATE");
661 fHMean.WriteToFile();
662 f->Close();
663 }
664 fHMean.fHM.SaveCanvasToImage(fOutputDir, "png");
665}
666
667void LmvmFastSim::CheckMemory(const string& text)
668{
669 static ProcInfo_t info;
670 const float toMB = 1.f / 1024.f;
671 gSystem->GetProcInfo(&info);
672 cout << Form((text + ": Memory res | vir: %g | %g Mbytes\n").c_str(), info.fMemResident * toMB,
673 info.fMemVirtual * toMB);
674 cout << endl;
675}
676
677ClassImp(LmvmFastSim);
CbmHistManager.h
Histogram manager.
ClassImp
ClassImp(LmvmFastSim)
LmvmFastSim::fPlutoMode
std::string fPlutoMode
Definition LmvmFastSim.h:33
LmvmFastSim::fNofFastSimEv
int fNofFastSimEv
Definition LmvmFastSim.h:35
LmvmFastSim::CalcMinvAndFillHist
void CalcMinvAndFillHist(const std::vector< LmvmDataFastSim > &Plus, const std::vector< LmvmDataFastSim > &Minus, const std::string &hName)
LmvmFastSim::fCatNames
std::vector< std::string > fCatNames
Definition LmvmFastSim.h:41
LmvmFastSim::fH
std::vector< LmvmHist * > fH
Definition LmvmFastSim.h:53
LmvmFastSim::CreateMeanHist
void CreateMeanHist(const std::string &name)
LmvmFastSim::fDoChargeBased
bool fDoChargeBased
Definition LmvmFastSim.h:45
LmvmFastSim::CheckHistogramNames
void CheckHistogramNames()
Definition LmvmFastSim.cxx:359
LmvmFastSim::GetChargeCategoryVector
void GetChargeCategoryVector()
Definition LmvmFastSim.cxx:404
LmvmFastSim::fOutputDir
std::string fOutputDir
Definition LmvmFastSim.h:31
LmvmFastSim::fDoParticleBased
bool fDoParticleBased
Definition LmvmFastSim.h:46
LmvmFastSim::CheckMemory
void CheckMemory(const std::string &text)
Definition LmvmFastSim.cxx:667
LmvmFastSim::fChargeCats
std::vector< std::string > fChargeCats
Definition LmvmFastSim.h:42
LmvmFastSim::GetIndividualCategories
void GetIndividualCategories()
Definition LmvmFastSim.cxx:421
LmvmFastSim::fDoMixEvents
bool fDoMixEvents
Definition LmvmFastSim.h:50
LmvmFastSim::fBgHistName
std::string fBgHistName
Definition LmvmFastSim.h:43
LmvmFastSim::fUrqmdNames
std::vector< std::string > fUrqmdNames
Definition LmvmFastSim.h:40
LmvmFastSim::fHMean
LmvmHist fHMean
Definition LmvmFastSim.h:54
LmvmFastSim::ContainsString
bool ContainsString(const std::vector< std::string > &vName, const std::string &name)
Definition LmvmFastSim.cxx:484
LmvmFastSim::fPlutoNames
std::vector< std::string > fPlutoNames
Definition LmvmFastSim.h:39
LmvmFastSim::GetPlutoAndUrqmdNames
void GetPlutoAndUrqmdNames(const std::vector< std::string > &files)
Definition LmvmFastSim.cxx:439
LmvmFastSim::CreateHistos
void CreateHistos()
Definition LmvmFastSim.cxx:489
LmvmFastSim::fOutFileName
std::string fOutFileName
Definition LmvmFastSim.h:32
LmvmFastSim::DoFastSim
void DoFastSim(std::vector< std::string > files, const std::string &dataDir, const std::string &outFile, int nofFastSimEv, const std::string &bgHistName)
Definition LmvmFastSim.cxx:29
Cbm::Split
vector< string > Split(const string &name, char delimiter)
Definition CbmUtils.cxx:66
std
Hash for CbmL1LinkKey.
Definition algo/base/Options.cxx:22