Heavy ion Analysis Libriares
Loading...
Searching...
No Matches
hal_corrfit.cpp
1/*
2 * nica_corrfit_db.cpp
3 *
4 * Created on: 4 kwi 2020
5 * Author: Daniel Wielanek
6 * E-mail: daniel.wielanek@gmail.com
7 * Warsaw University of Technology, Faculty of Physics
8 */
9#include "AnaFile.h"
10#include "CorrFitConst.h"
11#include "CorrFitInfo.h"
12#include "CorrFitMapRebin.h"
13#include "CorrFitParamsSetup.h"
14#include "Cout.h"
15#include "Femto1DCF.h"
16#include "Femto3DCF.h"
17#include "FemtoSHCF.h"
18#include "FemtoSerializationInterface.h"
19#include "Package.h"
20#include "PackageSql.h"
21
22#include <TChain.h>
23#include <TError.h>
24#include <TFile.h>
25#include <TKey.h>
26#include <TString.h>
27#include <TSystem.h>
28#include <TTree.h>
29
30#include <fstream>
31
32#include <fstream>
33#include <iostream>
34
35void merge();
36void mergeVertical();
37void mergeHorizontal();
38void printjobs();
39void printjobsbad();
40void preparetemplate();
41void compress();
42int getNFiles();
43void rebin(TString opt);
44int main(int argc, char* argv[]) {
45 if (argc < 2) {
46 Hal::Cout::PrintInfo("No arguments! run: hal-corrfit --help to get help", Hal::EInfo::kCriticalError);
47 return 0;
48 }
49 TString arg1 = argv[1];
50 if (arg1 == "--help") {
51 std::cout << " Usage:" << std::endl;
52 std::cout << "--help print help" << std::endl;
53 std::cout << "0 - prepare template files and database" << std::endl;
54 std::cout << "1 - merge maps" << std::endl;
55 std::cout << "2 - compress maps" << std::endl;
56 std::cout << "--print - check/print info about executed jobs " << std::endl;
57 std::cout << "--print-bad - check/print info about executed jobs list of bad jobs is set to bad_jobs.txt" << std::endl;
58 std::cout << "--rebin=n rebin CF new_map.root fille will be created" << std::endl;
59 return 1;
60 } else if (arg1 == "0") { // prepare macro & xml files
61 preparetemplate();
62 return 1;
63 } else if (arg1 == "--print") {
64 printjobs();
65 } else if (arg1 == "--print-bad") {
66 printjobsbad();
67 } else if (arg1 == "1") {
68 printjobs();
69 merge();
70 } else if (arg1 == "2") {
71 compress();
72 } else if (arg1.BeginsWith("--rebin=")) {
73 rebin(arg1);
74 } else {
75 std::cout << "Unknown parameter type --help" << std::endl;
76 return 1;
77 }
78 return 0;
79}
80
81Bool_t AllJobsCompleted(Int_t jobs) {
82 for (int i = 0; i < jobs; i++) {
83 if (Hal::CorrFitParamsSetup::TestMapFile(i) == kFALSE) return kFALSE;
84 }
85 return kTRUE;
86}
87
88inline void merge() {
89 TFile* info_file = new TFile("files/config.root");
90 Hal::CorrFitInfo* obj = new Hal::CorrFitInfo(*(Hal::CorrFitInfo*) info_file->Get("Info"));
91 bool vert = obj->IsVertical();
92 info_file->Close();
93 vert = true;
94 if (vert) {
95 mergeVertical();
96 } else {
97 mergeHorizontal();
98 }
99}
100
101void printjobs() {
102 Int_t nJobs = getNFiles();
103 Bool_t all_jobs = AllJobsCompleted(nJobs);
104
105 Hal::Cout::Text(Form("CorrFitDB %i jobs", nJobs), "M", kWhite);
106 if (all_jobs) {
107 Hal::Cout::Text("All jobs completed", "", kWhite);
108 } else {
109 Hal::Cout::Text("Not all jobs are completed", "", kWhite);
110 }
111 for (int i = 0; i < nJobs; i++) {
112 Bool_t fileOk = Hal::CorrFitParamsSetup::TestMapFile(i);
113 if (fileOk) {
114 Hal::Cout::FailSucced(Form("job %i", i), "GOOD", kCyan);
115 } else {
116 Hal::Cout::FailSucced(Form("job %i", i), "BAD ", kRed);
117 }
118 }
119}
120
121void printjobsbad() {
122 Hal::CorrFitParamsSetup setup("corrfit_conf.xml");
123 Int_t nJobs = setup.GetNJobs();
124 Bool_t all_jobs = AllJobsCompleted(nJobs);
125 Hal::Cout::Text(Form("CorrFitDB %i jobs", nJobs), "M", kWhite);
126 if (all_jobs) {
127 Hal::Cout::Text("All jobs completed", "", kWhite);
128 } else {
129 Hal::Cout::Text("Not all jobs are completed", "", kWhite);
130 }
131 std::ofstream bad_jobs;
132 bad_jobs.open("bad_jobs.txt");
133 bad_jobs << "LIST_OF_BAD_JOBS" << std::endl;
134 int bad_job_id = 0;
135 for (int i = 0; i < nJobs; i++) {
136 Bool_t fileOk = Hal::CorrFitParamsSetup::TestMapFile(i);
137 if (!fileOk) {
138 Hal::Cout::FailSucced(Form("job %i", i), "BAD ", kRed);
139 bad_jobs << "BAD_JOB jobs/job_" << i << std::endl;
140 ++bad_job_id;
141 }
142 }
143 bad_jobs.close();
144}
145
146void preparetemplate() {
147
148 std::ofstream xml_file;
149 xml_file.open("corrfit_conf.xml");
150 xml_file << "<CorrfitConfig>" << std::endl;
151 xml_file << "<Parameters>" << std::endl;
152 xml_file << "\t<Param name=\"R_{out}\" min=\"1\" max=\"10\" step=\"1\"></Param>" << std::endl;
153 xml_file << "\t<Param name=\"R_{side}\" min=\"1\" max=\"10\" step=\"1\"></Param>" << std::endl;
154 xml_file << "\t<Param name=\"R_{long}\" min=\"1\" max=\"10\" step=\"1\"></Param>" << std::endl;
155 xml_file << "</Parameters>" << std::endl;
156 xml_file << "<!-- full path to file with pairs-->" << std::endl;
157 xml_file << "<PairFile>zz.root</PairFile>" << std::endl;
158 xml_file << "<!-- optional part, use to configure dump pair analysis-->" << std::endl;
159 xml_file << "<DumpAnalysisConf>" << std::endl;
160 xml_file << "\t<CorrelationFunction>" << std::endl;
161 xml_file << "\t\t<Name>CF</Name>" << std::endl;
162 xml_file << "\t\t<Frame>EKinematics::kLCMS</Frame>" << std::endl;
163 xml_file << "\t\t<Type>Femto3DCF</Type>" << std::endl;
164 xml_file << "<!-- optional part, used only for spherical harmonics-->" << std::endl;
165 xml_file << "\t\t<L>3</L>" << std::endl;
166 xml_file << "\t\t<Xaxis bins=\"100\" min=\"0.0\" max=\"1.0\"></Xaxis>" << std::endl;
167 xml_file << "\t\t<Yaxis bins=\"100\" min=\"0.0\" max=\"1.0\"></Yaxis>" << std::endl;
168 xml_file << "\t\t<Zaxis bins=\"100\" min=\"0.0\" max=\"1.0\"></Zaxis>" << std::endl;
169 xml_file << "\t</CorrelationFunction>" << std::endl;
170 xml_file << "\t<FreezoutGenerator>Hal::FemtoFreezoutGeneratorLCMS</FreezoutGenerator>" << std::endl;
171 xml_file << "\t<SourceModel>Hal::FemtoSourceModelGauss3D</SourceModel>" << std::endl;
172 xml_file << "\t<CalcOptions>" << std::endl;
173 xml_file << "\t\t<JobMultiplyFactor>1</JobMultiplyFactor>" << std::endl;
174 xml_file << "\t\t<WeightMultiplyFactor>1</WeightMultiplyFactor>" << std::endl;
175 xml_file << "\t\t<PreprocessMultiplyFactor>1</PreprocessMultiplyFactor>" << std::endl;
176 xml_file << "\t\t<!-- S/B/S+B for S(signal) B (background) B+S (both)-->" << std::endl;
177 xml_file << "\t\t<CalcMode>S</CalcMode>" << std::endl;
178 xml_file << "\t\t<IgnoreSign>kTRUE</IgnoreSign>" << std::endl;
179 xml_file << "\t</CalcOptions>" << std::endl;
180 xml_file << "\t<WeightConf>" << std::endl;
181 xml_file << "\t\t<Type>Hal::FemtoWeightGeneratorLednicky</Type>" << std::endl;
182 xml_file << "\t\t<QuantumOn>kTRUE</QuantumOn>" << std::endl;
183 xml_file << "\t\t<StrongOn>kFALSE</StrongOn>" << std::endl;
184 xml_file << "\t\t<CoulombOn>kFALSE</CoulombOn>" << std::endl;
185 xml_file << "\t\t<PairType>211;211</PairType>" << std::endl;
186 xml_file << "\t</WeightConf>" << std::endl;
187 xml_file << "</DumpAnalysisConf>" << std::endl;
188 xml_file << "</CorrfitConfig>" << std::endl;
189 xml_file.close();
190
191 std::ofstream macro;
192 macro.open("ana.C");
193 macro << "void ana(){" << std::endl;
194 macro << "\tTString jobEnv = gSystem->Getenv(\"JOB_ID_HAL\");" << std::endl;
195 macro << "\tInt_t job_index = jobEnv.Atoi();" << std::endl;
196 macro << "\tHal::CorrFitDumpedPairAna *task = new Hal::CorrFitDumpedPairAna(job_index);" << std::endl;
197 macro << "\tif(task->Init()==kFALSE) return;" << std::endl;
198 macro << "\ttask->Run();" << std::endl;
199 macro << "\ttask->Finish();" << std::endl;
200 macro << "}" << std::endl;
201 macro.close();
202
203 gSystem->mkdir("jobs");
204 gSystem->mkdir("files");
205}
206
207void rebin(TString rebin) {
208
209 TString opt = Hal::Std::RemoveString(rebin, "--rebin=");
210 Int_t Nrebin = opt.Atoi();
211 if (Nrebin <= 0) {
212 std::cout << "invalid number of bins " << std::endl;
213 return;
214 }
215 Hal::CorrFitMapRebin* rebinner = new Hal::CorrFitMapRebin("merged_map.root", Nrebin);
216 rebinner->Compress("merged_rebined_map.root");
217 delete rebinner;
218}
219int getNFiles() {
220 TFile* info_file = new TFile("files/config.root");
221 Hal::CorrFitInfo* obj = new Hal::CorrFitInfo(*(Hal::CorrFitInfo*) info_file->Get("Info"));
222 Hal::CorrFitParamsSetup setup("corrfit_conf.xml");
223 int nFiles = setup.GetNJobs();
224 if (obj->IsVertical()) { nFiles = obj->GetMainBins(); }
225 return nFiles;
226}
227void mergeHorizontal() {
228 TFile* info_file = new TFile("files/config.root");
229 Hal::CorrFitInfo* obj = new Hal::CorrFitInfo(*(Hal::CorrFitInfo*) info_file->Get("Info"));
230 Hal::CorrFitParamsSetup setup("corrfit_conf.xml");
231 int nJobs = setup.GetNJobs();
232 std::cout << "Merging vertical CF's" << nJobs << std::endl;
233 if (!AllJobsCompleted(nJobs)) {
234 Hal::Cout::PrintInfo("Not all jobs are completed", Hal::EInfo::kError);
235 Hal::Cout::Text("Do you really want to continue y/n?", "L");
236 TString ans;
237 std::cin >> ans;
238 if (!(ans == "y" || ans == "Y")) { return; }
239 }
240 TFile* file = new TFile("merged_map.root", "recreate");
241 TChain* chain = new TChain("map", "");
242
243 for (int i = 0; i < nJobs; i++) {
244 chain->AddFile(Form("files/corrfit_map_%i.root", i));
245 }
246
247 TTree* tree = new TTree("map", "");
248 Hal::Array_1<Float_t>* Data_out = new Hal::Array_1<Float_t>();
249 Hal::Array_1<Float_t>* Data_in = nullptr;
250 chain->SetBranchAddress("data", &Data_in);
251 tree->Branch("data", &Data_out);
252 int entries = chain->GetEntries() / nJobs; // n_files *n_bins /nJobs(nbins) = n_files
253 chain->GetEntry(0);
254 Data_out->MakeBigger(Data_in->GetSize() * nJobs);
255 int single_size = Data_in->GetSize();
256 for (int i = 0; i < nJobs; i++) { // nsamples
257 chain->GetEntry(i);
258 for (int k = 0; k < single_size; k++) {
259 Data_out->Set(k, Data_in->Get(k));
260 }
261 tree->Fill();
262 }
263
264 tree->Write();
265 obj->Write();
266 file->Close();
267 info_file->Close();
268}
269void mergeVertical() {
270 TFile* info_file = new TFile("files/config.root");
271 Hal::CorrFitInfo* obj = new Hal::CorrFitInfo(*(Hal::CorrFitInfo*) info_file->Get("Info"));
272 bool vert = obj->IsVertical();
273 auto CF = obj->GetCf()->Clone();
274 std::cout << "CF " << CF->ClassName() << std::endl;
275 auto cfType = Hal::Femto::GetCFType(CF);
276 int nSlices = Hal::Femto::GetNSlices(CF);
277
278
279 Hal::CorrFitParamsSetup setup("corrfit_conf.xml");
280 int n_params = setup.GetNJobs();
281 if (!AllJobsCompleted(nSlices)) {
282 Hal::Cout::PrintInfo("Not all nFiles are completed", Hal::EInfo::kError);
283 Hal::Cout::Text("Do you really want to continue y/n?", "L");
284 TString ans;
285 std::cin >> ans;
286 if (!(ans == "y" || ans == "Y")) { return; }
287 }
288
289 TFile* file = new TFile("merged_map.root", "recreate");
290 TChain* chain = new TChain("map", "");
291
292 for (int i = 0; i < nSlices; i++) {
293 chain->AddFile(Form("files/corrfit_map_%i.root", i));
294 }
295 std::cout << "Merging vertical CF's" << nSlices << " " << n_params << std::endl;
296 auto* tree = new TTree("map", "");
297 auto* Data_out = new Hal::Array_1<Float_t>();
298 Hal::Array_1<Float_t>* Data_in = nullptr;
299 chain->SetBranchAddress("data", &Data_in);
300 tree->Branch("data", &Data_out);
301 chain->GetEntry(0);
302 const int slice_size = Data_in->GetSize();
303 std::cout << "CHECK SIZE " << nSlices * n_params << " " << chain->GetEntries() << std::endl;
304 std::cout << nSlices << " " << n_params << std::endl;
305 Data_out->MakeBigger(slice_size * nSlices);
306 for (int iPar = 0; iPar < n_params; iPar++) {
307 for (int iBin = 0; iBin < nSlices; iBin++) {
308 chain->GetEntry(iPar + n_params * iBin);
309 for (int iS = 0; iS < slice_size; iS++) {
310 Data_out->Set(iBin * slice_size + iS, Data_in->Get(iS));
311 }
312 }
313 tree->Fill();
314 }
315 delete CF;
316 tree->Write();
317 obj->Write();
318 file->Close();
319 info_file->Close();
320}
321
322void compress() {
323 TFile* info_file = new TFile("files/config.root");
324 Hal::CorrFitInfo* obj = new Hal::CorrFitInfo(*(Hal::CorrFitInfo*) info_file->Get("Info"));
325 auto CF = (Hal::DividedHisto1D*) obj->GetCf()->Clone();
326 auto cfType = Hal::Femto::GetCFType(CF);
327 TFile* fileIn = new TFile("merged_map.root");
328 TTree* treeIn = (TTree*) fileIn->Get("map");
329
330 TFile* fileOut = new TFile("compressed_map.root", "recreate");
331 TTree* treeOut = new TTree("map", "");
332 Hal::Array_1<Float_t>* Data_out = new Hal::Array_1<Float_t>();
333 Hal::Array_1<Float_t>* Data_in = nullptr;
334 treeIn->SetBranchAddress("data", &Data_in);
335 treeOut->Branch("data", &Data_out);
336 auto interfaceFrom = (Hal::FemtoSerializationInterface*) CF->GetSpecial("serialization");
337 std::cout << interfaceFrom << std::endl;
338 std::cout << interfaceFrom->ClassName() << std::endl;
339 interfaceFrom->BindCFs(CF);
340 interfaceFrom->BindArray(Data_in);
341 interfaceFrom->SetOption(Hal::FemtoSerializationInterface::EOption::kFull);
342 interfaceFrom->Init();
343 auto interfaceTo = (Hal::FemtoSerializationInterface*) CF->GetSpecial("serialization");
344 interfaceTo->BindCFs(CF);
345 interfaceTo->BindArray(Data_out);
346 interfaceTo->SetOption(Hal::FemtoSerializationInterface::EOption::kSimple);
347 interfaceTo->Init();
348 const int entries = treeIn->GetEntries();
349 obj->Write();
350 auto sh = dynamic_cast<Hal::FemtoSHCF*>(CF);
351 for (int i = 0; i < entries; i++) {
352 treeIn->GetEntry(i);
353 interfaceFrom->Deserialize();
354 if (i == 0) std::cout << "ARSIZE " << Data_in->GetSize() << " " << Data_out->GetSize() << std::endl;
355 if (sh) {
356 sh->RecalculateCF();
357 sh->RecalculateCF(-1, kTRUE);
358 /*std::cout << "NUM " << sh->GetNum()->GetBinContent(1) << std::endl;
359 std::cout << "NUM " << sh->GetNum()->GetBinContent(2) << std::endl;
360 std::cout << "NUM " << sh->GetNum()->GetBinContent(3) << std::endl;
361 TFile* f = new TFile(Form("dump_%i.root", i), "recreate");
362 auto copy = sh->Clone();
363 f->cd();
364 copy->Write();
365
366 f->Close();*/
367 }
368 interfaceTo->Serialize();
369 treeOut->Fill();
370 }
371 treeOut->Write();
372 fileIn->Close();
373 fileOut->Close();
374 delete interfaceFrom;
375 delete interfaceTo;
376}
Hal::Array_1
Definition Array.h:19
Hal::Array_1::Set
void Set(Int_t index, T val)
Definition Array.h:103
Hal::Array_1::Get
T Get(Int_t i) const
Definition Array.h:97
Hal::Array_1::MakeBigger
void MakeBigger(Int_t new_dim)
Definition Array.cxx:5
Hal::Array_1::GetSize
Int_t GetSize() const
Definition Array.h:50
Hal::CorrFitMapRebin::Compress
Bool_t Compress(TString outFile)
Definition CorrFitMapRebin.cxx:34
Hal::Cout::Text
static void Text(TString text, TString option="L", Color_t color=-1)
Definition Cout.cxx:92
Hal::Cout::PrintInfo
static void PrintInfo(TString text, Hal::EInfo status)
Definition Cout.cxx:370
Hal::Cout::FailSucced
static void FailSucced(TString value, TString flag, Color_t color)
Definition Cout.cxx:382