Heavy ion Analysis Libriares
Loading...
Searching...
No Matches
CorrFitDumpedPairAna.cxx
1/*
2 * FemtoNDimMap.cxx
3 *
4 * Created on: 24 paź 2019
5 * Author: Daniel Wielanek
6 * E-mail: daniel.wielanek@gmail.com
7 * Warsaw University of Technology, Faculty of Physics
8 */
9
10#include "CorrFitDumpedPairAna.h"
11
12
13#include <TApplication.h>
14#include <TBranch.h>
15#include <TChain.h>
16#include <TClonesArray.h>
17#include <TCollection.h>
18#include <TDirectory.h>
19#include <TDirectoryFile.h>
20#include <TFile.h>
21#include <TKey.h>
22#include <TList.h>
23#include <TNamed.h>
24#include <TObjArray.h>
25#include <TSystem.h>
26#include <TTree.h>
27#include <cstdlib>
28#include <fstream>
29#include <string>
30
31#include "Array.h"
32#include "CorrFitInfo.h"
33#include "CorrFitParamsSetup.h"
34#include "Cout.h"
35#include "Femto1DCF.h"
36#include "Femto3DCF.h"
37#include "FemtoConst.h"
38#include "FemtoCorrFunc.h"
39#include "FemtoCorrFuncSimple.h"
40#include "FemtoDumpPairAna.h"
41#include "FemtoFreezoutGenerator.h"
42#include "FemtoMiniPair.h"
43#include "FemtoPair.h"
44#include "FemtoSHCF.h"
45#include "FemtoSerializationInterface.h"
46#include "FemtoSourceModel.h"
47#include "FemtoWeightGenerator.h"
48#include "FemtoWeightGeneratorLednicky.h"
49#include "XMLNode.h"
50
51namespace Hal {
52 CorrFitDumpedPairAna::CorrFitDumpedPairAna(Int_t jobid, Int_t maps_perAna) :
53 TObject(), fPairFile(""), fJobId(jobid), fMultiplyJobs(maps_perAna), fMode(eDumpCalcMode::kSignalPairs) {}
54
55 void CorrFitDumpedPairAna::SetCorrFunc(const FemtoCorrFunc& func) {
56 if (fTempCF) delete fTempCF;
57 fTempCF = (FemtoCorrFunc*) func.Clone();
58 }
59
60 void CorrFitDumpedPairAna::SetWeightGenerator(const FemtoWeightGenerator& weight) {
61 if (fWeight) delete fWeight;
62 fWeight = weight.MakeCopy();
63 }
64
65 Bool_t CorrFitDumpedPairAna::Init() {
66 if (!ConfigureFromXML()) {
67 Hal::Cout::PrintInfo("Cannot configure XML input", EInfo::kError);
68 return kFALSE;
69 }
70 Bool_t fileOk = CorrFitParamsSetup::TestMapFile(fJobId * fMultiplyJobs);
71 if (fileOk) {
72 Hal::Cout::PrintInfo("Test file found no need to calculate map", EInfo::kError);
73 exit(0); // force exit to prevent crash
74 return kFALSE;
75 }
76 if (!ConfigureInput()) {
77 Hal::Cout::PrintInfo("Cannot configure input", EInfo::kError);
78 return kFALSE;
79 }
80 if (!fWeight) {
81 Hal::Cout::PrintInfo("Cannot find weight", EInfo::kError);
82 return kFALSE;
83 }
84 if (!fTempCF) {
85 Hal::Cout::PrintInfo("Cannot find CF", EInfo::kError);
86 return kFALSE;
87 }
88 if (fTempCF->GetEntries() > 1) {
89 Hal::Cout::PrintInfo("To many CF", EInfo::kError);
90 return kFALSE;
91 }
92 if (!fTree) {
93 Hal::Cout::PrintInfo("Cannot find tree", EInfo::kError);
94 return kFALSE;
95 }
96 if (fJobId == -1) {
97 Hal::Cout::PrintInfo("Wrong job iD", EInfo::kError);
98 return kFALSE;
99 }
100 DividedHisto1D* dummy = fTempCF->GetCF(0);
101 Femto::EKinematics kinem = Femto::EKinematics::kLCMS;
102 if (dummy->GetLabelsNo() > 0) {
103 TString label = dummy->GetLabel(0);
104 kinem = Femto::LabelToKinematics(label);
105 }
106 fPair = Femto::MakePair(kinem, fImgMom);
107
108 if (fIgnoreSing) {
109 if (dummy->InheritsFrom("Hal::FemtoSHCF")) {
110 Cout::PrintInfo("Cannot ignore sign with SH correlation function!", EInfo::kError);
111 } else {
112 fPair->UseAbs();
113 }
114 }
115 return kTRUE;
116 }
117
118 void CorrFitDumpedPairAna::Run(Int_t events) {
119 Int_t nEvents = 0;
120 if (events > 0) {
121 nEvents = events;
122 } else {
123 nEvents = fTree->GetEntries();
124 }
125 switch (fMode) {
126 case eDumpCalcMode::kSignalBackgroundPairs: {
127 Int_t step = nEvents / 100;
128 for (int iEvent = 0; iEvent < nEvents; iEvent++) {
129 fTree->GetEntry(iEvent);
130 if (iEvent % step == 0) { Cout::ProgressBar(iEvent, nEvents); }
131 RunSignalBackgroundPair();
132 if (fPairThreshold > 0 && fPairThreshold < fPairsProcessed) break;
133 }
134
135 } break;
136 case eDumpCalcMode::kSignalPairs: {
137 RunSignalPair();
138 Int_t step = nEvents / 100;
139 for (int iEvent = 0; iEvent < nEvents; iEvent++) {
140 fTree->GetEntry(iEvent);
141 if (iEvent % step == 0) { Cout::ProgressBar(iEvent, nEvents); }
142 RunSignalPair();
143 if (fPairThreshold > 0 && fPairThreshold < fPairsProcessed) break;
144 }
145 } break;
146 case eDumpCalcMode::kBackgroundPairsOnly: {
147 Int_t step = nEvents / 100;
148 for (int iEvent = 0; iEvent < nEvents; iEvent++) {
149 fTree->GetEntry(iEvent);
150 if (iEvent % step == 0) { Cout::ProgressBar(iEvent, nEvents); }
151 RunBackgroundPair();
152 if (fPairThreshold > 0 && fPairThreshold < fPairsProcessed) break;
153 }
154 } break;
155 }
156 }
157
158 Bool_t CorrFitDumpedPairAna::SaveAsRawArray(TObject* cf, Int_t step) {
159 TFile* file = new TFile(Form("files/corrfit_map_%i.root", GetSimStepNo() + step), "recreate");
160 auto cfx = (Hal::DividedHisto1D*) cf;
161 auto serialization = (FemtoSerializationInterface*) cfx->GetSpecial("serialization");
162 serialization->BindCFs(cf);
163 Array_1<Float_t>* Data = new Array_1<Float_t>();
164 serialization->BindArray(Data);
165 serialization->Init();
166 serialization->Serialize();
167 if (!Data) return kFALSE;
168 TTree* tree = new TTree("map", "map");
169 tree->Branch("data", &Data);
170 tree->Fill();
171 tree->Write();
172 file->Close();
173 delete serialization;
174 return kTRUE;
175 }
176
177 Bool_t CorrFitDumpedPairAna::ConfigureInput() {
178 if (fPairFile.EndsWith(".root")) {
179 return ConfigureRootInput();
180 } else {
181 ConfigureListInput();
182 return kFALSE;
183 }
184 }
185
186 Bool_t CorrFitDumpedPairAna::ConfigureFromXML() {
187 auto printErr = [](TString str) {
188 Hal::Cout::PrintInfo(Form("Cannot find node: %s in XML", str.Data()), EInfo::kError);
189 return kFALSE;
190 };
191 XMLFile file("corrfit_conf.xml");
192 CorrFitParamsSetup setup("corrfit_conf.xml");
193 fTotalNumberOfPoints = 1;
194 int parNo = setup.GetNParams();
195 for (int i = 0; i < parNo; i++) {
196 fTotalNumberOfPoints = fTotalNumberOfPoints * setup.GetNPoints(i);
197 }
198 XMLNode* root = file.GetRootNode();
199 fPairFile = root->GetChild("PairFile")->GetValue();
200 XMLNode* parameters = root->GetChild("Parameters");
201 XMLNode* dumpAna = root->GetChild("DumpAnalysisConf");
202 if (!dumpAna) return printErr("DumpAnalysisConf");
203
204 if (!parameters) return printErr("Parameters");
205
206 XMLNode* freezXml = dumpAna->GetChild("FreezoutGenerator");
207 XMLNode* sourceXml = dumpAna->GetChild("SourceModel");
208 if (!freezXml) return printErr("FreezoutGenerator");
209
210 if (!sourceXml) return printErr("SourceModel");
211 TClass* freezoutClass = TClass::GetClass(freezXml->GetValue(), 1, 0);
212 TClass* sourceClass = TClass::GetClass(sourceXml->GetValue(), 1, 0);
213 FemtoFreezoutGenerator* generator = static_cast<FemtoFreezoutGenerator*>(freezoutClass->New());
214 FemtoSourceModel* source = static_cast<FemtoSourceModel*>(sourceClass->New());
215 if (source && generator) {
216 generator->SetSourceModel(*source);
217 SetGenerator(*generator);
218 } else {
219 Hal::Cout::PrintInfo("Cannot create source or generator class", EInfo::kError);
220 return kFALSE;
221 }
222 if (source) delete source;
223 if (generator) delete generator;
224 XMLNode* calcOpts = dumpAna->GetChild("CalcOptions");
225 if (calcOpts) {
226 XMLNode* pairCut = calcOpts->GetChild("NoPairCut");
227 if (pairCut) { fPairThreshold = pairCut->GetValue().Atoi(); }
228 Hal::Cout::PrintInfo("Pair threshold detected", EInfo::kInfo);
229 }
230 XMLNode* multiJobs = calcOpts->GetChild("JobMultiplyFactor");
231 if (multiJobs) {
232 if (multiJobs->GetValue().Length() > 0 && fMultiplyJobs <= 0) { fMultiplyJobs = multiJobs->GetValue().Atoi(); }
233 }
234 if (fMultiplyJobs <= 0) fMultiplyJobs = 1;
235 std::vector<int> dims = setup.GetDimensions();
236 if (!InitGenerators(dims, parameters, setup)) {
237 Cout::PrintInfo("Cannot init generators !", EInfo::kError);
238 return kFALSE;
239 }
240 if (calcOpts) {
241 XMLNode* multiplyXmlWeight = calcOpts->GetChild("WeightMultiplyFactor");
242 if (multiplyXmlWeight) SetMultiplyFactorWeight(multiplyXmlWeight->GetValue().Atoi());
243
244 XMLNode* multiplyXmlPreproc = calcOpts->GetChild("PreprocessMultiplyFactor");
245 if (multiplyXmlPreproc) SetMultiplyFactorPreprocess(multiplyXmlPreproc->GetValue().Atoi());
246
247 XMLNode* calcXml = calcOpts->GetChild("CalcMode");
248 if (calcXml) {
249 TString opt = calcXml->GetValue();
250 if (opt == "S") fMode = eDumpCalcMode::kSignalPairs;
251 if (opt == "S+B") fMode = eDumpCalcMode::kSignalBackgroundPairs; // UseMixed
252 if (opt == "B") fMode = eDumpCalcMode::kBackgroundPairsOnly; // UseMixedOnly
253 }
254 XMLNode* ignore = calcOpts->GetChild("IgnoreSign");
255 if (ignore) {
256 if (ignore->GetValue() == "kTRUE") this->IgnoreSign();
257 }
258 }
259
260 XMLNode* weight = dumpAna->GetChild("WeightConf");
261 if (weight) {
262 auto w = Hal::Femto::GetWeightGeneratorFromXLM(weight);
263 if (!w) {
264 Cout::PrintInfo("Problem with weight configuration !", EInfo::kError);
265 return kFALSE;
266 }
267 this->SetWeightGenerator(*w);
268 delete w;
269 }
270 XMLNode* cfXml = dumpAna->GetChild("CorrelationFunction");
271 DividedHisto1D* cf = Femto::GetHistoFromXML(cfXml);
272 if (cf) {
273 FemtoCorrFuncSimple corrFunc(*cf);
274 delete cf;
275 this->SetCorrFunc(corrFunc);
276 if (!InitCFs()) { Cout::PrintInfo("Cannot init CF !", EInfo::kError); }
277 } else {
278 Cout::PrintInfo("Cannot find CF in XML !", EInfo::kError);
279 return kFALSE;
280 }
281 return kTRUE;
282 }
283
284 CorrFitDumpedPairAna::~CorrFitDumpedPairAna() {
285 if (fTempCF) delete fTempCF;
286 if (fTempGenerator) delete fTempGenerator;
287 for (auto i : fGenerator) {
288 delete i;
289 }
290 }
291
292 void CorrFitDumpedPairAna::Print(Option_t* /*option*/) const {
293 Cout::Text("CorrFitDumpedPairAna Info", "M");
294 Cout::Text(Form("JobID: %i", fJobId), "L");
295 Cout::Text(Form("Multiply factor jobs: %i", fMultiplyJobs), "L");
296 Cout::Text(Form("Multiply factor weight: %i", fMultiplyWeight), "L");
297 Cout::Text(Form("Multiply factor preprocess: %i", fMultiplyPreprocess), "L");
298 Cout::Text(Form("IgnoreSign: %i", fIgnoreSing), "L");
299 Cout::Text(Form("Pair file: %s", fPairFile.Data()), "L");
300 switch (fMode) {
301 case eDumpCalcMode::kBackgroundPairsOnly: {
302 Cout::Text("Use only background branch", "L");
303 } break;
304 case eDumpCalcMode::kSignalBackgroundPairs: {
305 Cout::Text("Use signal & background branches", "L");
306 } break;
307 case eDumpCalcMode::kSignalPairs: {
308 Cout::Text("Use only signal branch", "L");
309 } break;
310 }
311 Cout::Text("CF Info", "M");
312 if (fTempCF) fTempCF->Print();
313 Cout::Text("Generator Info", "M");
314
315 for (auto& gen : fGenerator) {
316 if (gen) gen->Print();
317 }
318 Cout::Text("Weight Info", "M");
319 if (fWeight) fWeight->Print();
320 }
321
322 void CorrFitDumpedPairAna::ConnectToSignal(const std::vector<TString>& branches) {
323 for (auto name : branches) {
324 fUsedBranches.push_back(name);
325 TClonesArray* array = new TClonesArray("Hal::FemtoMicroPair");
326 fTree->SetBranchAddress(name, &array);
327 fSignalClones.push_back(array);
328 }
329 }
330
331 void CorrFitDumpedPairAna::ConnectToBackground(const std::vector<TString>& branches) {
332 for (auto name : branches) {
333 fUsedBranches.push_back(name);
334 TClonesArray* array = new TClonesArray("Hal::FemtoMicroPair");
335 fTree->SetBranchAddress(name, &array);
336 fBackgroundClones.push_back(array);
337 }
338 }
339
340 void CorrFitDumpedPairAna::LockUnusedBranches() {
341 auto list = fTree->GetListOfBranches();
342 for (int i = 0; i < list->GetEntries(); i++) {
343 auto branch = (TBranch*) list->At(i);
344 Bool_t lock = kTRUE;
345 TString brName = branch->GetName();
346 for (auto name : fUsedBranches)
347 if (name == brName) {
348 lock = kFALSE;
349 continue;
350 }
351 if (lock) fTree->SetBranchStatus(brName, 0);
352 }
353 }
354
355 Bool_t CorrFitDumpedPairAna::ConfigureRootInput() {
356 fFile = new TFile(fPairFile);
357 fGrouping = dynamic_cast<CorrFitMapGroupConfig*>(fFile->Get("HalInfo/CorrFitMapGroup"));
358 if (!fGrouping) {
359 Hal::Cout::PrintInfo("Cannot find grouping config", EInfo::kError);
360 return kFALSE;
361 }
362 TString treeName = FindTreeName(fPairFile);
363 fTree = new TChain(treeName);
364 fTree->AddFile(fPairFile);
365
366 if (treeName.Length() == 0) return kFALSE;
367 Bool_t res = ConnectToData();
368 LockUnusedBranches();
369 return res;
370 }
371
372 Bool_t CorrFitDumpedPairAna::ConfigureListInput() {
373 std::vector<TString> filelist;
374 if (fPairFile.EndsWith("/")) {
375 filelist = Hal::Std::GetListOfFiles(fPairFile, "root", kTRUE, kFALSE);
376 } else {
377 filelist = Hal::Std::GetLinesFromFile(fPairFile, kTRUE);
378 }
379 if (filelist.size() < 1) return kFALSE;
380 fFile = new TFile(filelist[0]);
381 fGrouping = dynamic_cast<CorrFitMapGroupConfig*>(fFile->Get("HalInfo/CorrFitMapGroup"));
382 if (!fGrouping) {
383 Hal::Cout::PrintInfo("Cannot find grouping config", EInfo::kError);
384 return kFALSE;
385 }
386 TString treeName = FindTreeName(filelist[0]);
387 fTree = new TChain(treeName);
388 for (auto file : filelist) {
389 fTree->AddFile(file);
390 }
391 return kTRUE;
392 }
393
394 TString CorrFitDumpedPairAna::FindTreeName(TString name) const {
395 TFile* f = new TFile(name);
396 TList* keys = f->GetListOfKeys();
397 TString treeName;
398 for (int i = 0; i < keys->GetEntries(); i++) {
399 TKey* key = (TKey*) keys->At(i);
400 TObject* obj = f->Get(key->GetName());
401 if (obj->InheritsFrom("TTree")) {
402 treeName = key->GetName();
403 break;
404 }
405 }
406 f->Close();
407 delete f;
408 return treeName;
409 }
410
411} // namespace Hal
Hal::Array_1
Definition Array.h:19
Hal::Cout::PrintInfo
static void PrintInfo(TString text, Hal::EInfo status)
Definition Cout.cxx:370
Hal::FemtoWeightGenerator::MakeCopy
virtual FemtoWeightGenerator * MakeCopy() const
Definition FemtoWeightGenerator.h:87