18 #include "TObjString.h"
20 #include "TObjArray.h"
21 #include "TObjString.h"
143 unsigned int nNull = 0;
145 bool ok = (
size > 0);
151 bool ok = (
size > 0);
157 cout << i1 <<
"-" << i2 << endl;
179 std::vector<unsigned int> hashV, indexV;
180 unsigned int nTot = 0;
184 if (!history)
continue;
185 for (
unsigned int j = 0; j < history->
nData(); j++) {
188 if (
data.isDisconnected())
continue;
189 if (nTot % 10000 == 0) cout << nTot << endl;
190 if (
data.energy() > eCut) {
191 cout <<
"E = " <<
data.energy() <<
" " <<
i <<
" " << j << endl;
201 for (
unsigned int i = 0;
i < hashV.size();
i++) {
206 cout << hashV.size() <<
"/" << nTot << endl;
215 if (!history)
continue;
217 cout <<
"Invalid LArSamplesHistory at hash = " <<
i << endl;
228 std::vector<const Interface*> interfaces;
229 interfaces.push_back(
this);
230 interfaces.push_back(&
other);
238 for (
unsigned int i = 0;
i < interfaces.size();
i++)
239 accessors.push_back(&interfaces[
i]->accessor());
246 std::vector<const Interface*> interfaces;
247 interfaces.push_back(
this);
248 interfaces.push_back(&
other);
256 for (
unsigned int i = 0;
i < interfaces.size();
i++)
257 accessors.push_back(&interfaces[
i]->accessor());
266 if (!
multi)
return nullptr;
267 std::vector<const Interface*> justOne;
268 justOne.push_back(
multi);
275 if (!
multi)
return nullptr;
276 std::vector<const Interface*> justOne;
277 justOne.push_back(
multi);
286 TObjArray*
list = filters.Tokenize(
",;");
287 if (
list->GetEntries() == 0) {
288 cout <<
"No filtering specified, exiting.";
293 for (
int k = 0;
k <
list->GetEntries();
k++) {
294 TObjString* tobs = (TObjString*)(
list->At(
k));
295 TObjArray*
items = tobs->String().Tokenize(
":");
296 if (
items->GetEntries() != 2) {
297 cout <<
"Invalid filter entry " << tobs->String() <<
", exiting." << endl;
302 TString
params = ((TObjString*)(
items->At(0)))->String();
303 TString
suffix = ((TObjString*)(
items->At(1)))->String();
306 cout <<
"---" << endl;
314 if (!tweak.
set(tweaks))
return 0;
317 if (!
multi)
return 0;
320 std::vector<MultiTreeAccessor*> filtered_mts = mt->
filterComponents(filterList, tweak);
321 if (filtered_mts.size() != filterList.
size())
return 0;
323 cout <<
"Component filtering done!" << endl;
327 for (
unsigned int f = 0;
f < filtered_mts.size();
f++) {
328 std::vector<TString>
files;
329 for (
unsigned int i = 0;
i < filtered_mts[
f]->nAccessors();
i++) {
331 cout <<
"Added " <<
files.back() << endl;
333 delete filtered_mts[
f];
336 std::vector<const Interface*> justOne;
337 justOne.push_back(filtered_multi);
349 if (!
f.set(
sel))
return nullptr;
352 if (!tweak.
set(tweaks))
return nullptr;
357 if (newFN !=
"") thisFN = newFN;
359 return filter(
f, tweak, thisFN);
368 return rootName +
"_" +
suffix +
".root";
392 return filter(
f, tw, newFileName);
405 if (
info->region() != region)
continue;
419 if (
info->feb() != feb)
continue;
434 if (
ft >= 0 &&
info->feedThrough() !=
ft)
continue;
435 if (slot >= 0 &&
info->slot() != slot)
continue;
448 if (!
f.set(
sel))
return nullptr;
450 std::vector<TString> vars;
451 std::vector<DataFuncSet> funcs;
452 std::vector<DataFuncArgs>
args;
454 cout <<
"Invalid variable specification " <<
var << endl;
458 TString
title = vars[0];
461 TH1D*
h =
m.dist(funcs[0],
args[0], vars[0],
nBins, xMin, xMax,
462 title, vars[0],
"digits",
f);
463 if (!
h)
return nullptr;
470 int nBinsY,
double yMin,
double yMax,
471 const TString&
sel,
const TString&
opt)
const
475 if (!
f.set(
sel))
return nullptr;
477 std::vector<TString> vars;
478 std::vector<DataFuncSet> funcs;
479 std::vector<DataFuncArgs>
args;
481 cout <<
"Invalid variable specification " <<
varList << endl;
485 TString
title = vars[1] +
" vs. " + vars[0];
488 TH2D*
h =
m.dist(funcs[0],
args[0], funcs[1],
args[1], vars[0] +
"_" + vars[1],
489 nBinsX, xMin, xMax, nBinsY, yMin, yMax,
490 title, vars[0], vars[1],
f);
491 if (!
h)
return nullptr;
498 const TString&
sel,
const TString&
opt,
503 if (!
f.set(
sel))
return nullptr;
505 std::vector<TString> vars;
506 std::vector<DataFuncSet> funcs;
507 std::vector<DataFuncArgs>
args;
509 cout <<
"Invalid variable specification " <<
var << endl;
518 if (!
h)
return nullptr;
525 const TString&
sel,
const TString&
opt,
530 if (!
f.set(
sel))
return nullptr;
532 std::vector<TString> vars;
533 std::vector<DataFuncSet> funcs;
534 std::vector<DataFuncArgs>
args;
536 cout <<
"Invalid variable specification " <<
var << endl;
544 if (!
h)
return nullptr;
554 if (!
f.set(
sel))
return 0;
563 if (!
f.set(
sel))
return 0;
571 if (!history)
return false;
579 if (!
f.set(
sel))
return false;
582 if (!history)
continue;
586 if (hDesc ==
"")
continue;
587 cout << Form(
"Hash = %-5d : ",
i) << hDesc
588 <<
"-----------------------------------------------------------------------------"
598 if (!
f.set(
sel))
return false;
599 std::map< std::pair<unsigned int, unsigned int>,
unsigned int > eventCells;
600 std::map< std::pair<unsigned int, unsigned int>,
double > eventEnergy;
604 if ((
i+1) % 50000 == 0) cout <<
"Cell # " <<
i+1 <<
"/" <<
nChannels() << endl;
606 if (!history)
continue;
607 for (
unsigned int j = 0; j < history->
nData(); j++) {
608 std::pair<unsigned int, unsigned int>
ev = std::make_pair(history->
data(j)->
run(), history->
data(j)->
event());
614 for (
unsigned int i = 0;
i <
nEvents();
i++) {
616 std::pair<unsigned int, unsigned int>
id(evtData->
run(), evtData->
event());
619 printout += Form(
" : %6d LAr hits, %7.1f MeV", eventCells[
id], eventEnergy[
id]);
620 cout << printout << endl;
629 std::map<unsigned int, unsigned int>
events;
636 for (
unsigned int i = 0;
i <
nRuns();
i++) {
640 printout += Form(
" : %6d events",
events[rData->
run()]);
641 cout << printout << endl;
652 if (!
f.set(
sel))
return 0;
654 std::vector<TString> vars;
655 std::vector<DataFuncSet> funcs;
656 std::vector<DataFuncArgs>
args;
658 cout <<
"Invalid variable specification " <<
varList << endl;
662 TVectorD
mean(vars.size()), meanErr(vars.size());
663 TMatrixD
covMatrix(vars.size(), vars.size()), covMatrixErr(vars.size(), vars.size());
667 cout <<
"---------------------------" << endl;
668 for (
unsigned int i = 0;
i < vars.size();
i++)
669 cout << Form(
"| %10s | %-9.4g |", vars[
i].
Data(),
mean(
i)) << endl;
670 cout <<
"---------------------------" << endl << endl;
673 for (
unsigned int i = 0;
i < vars.size();
i++) cout <<
" |";
674 cout << endl <<
"--------------";
675 for (
unsigned int i = 0;
i < vars.size();
i++) cout <<
"-------------";
677 for (
unsigned int i1 = 0; i1 < vars.size(); i1++) {
678 cout << Form(
"| %10s |", vars[i1].
Data());
679 for (
unsigned int i2 = 0; i2 < vars.size(); i2++)
680 cout << Form(
" %-9.4g |",
covMatrix(i1, i2));
683 cout <<
"--------------";
684 for (
unsigned int i = 0;
i < vars.size();
i++) cout <<
"-------------";
688 cout <<
"---------------------------" << endl;
689 for (
unsigned int i = 0;
i < vars.size();
i++)
690 cout << Form(
"| %10s | %-9.4g +/- %-9.4g |", vars[
i].
Data(),
mean(
i), meanErr(
i)) << endl;
691 cout <<
"---------------------------" << endl << endl;
694 for (
unsigned int i = 0;
i < vars.size();
i++) cout <<
" |";
695 cout << endl <<
"--------------";
696 for (
unsigned int i = 0;
i < vars.size();
i++) cout <<
"---------------------------";
698 for (
unsigned int i1 = 0; i1 < vars.size(); i1++) {
699 cout << Form(
"| %10s |", vars[i1].
Data());
700 for (
unsigned int i2 = 0; i2 < vars.size(); i2++)
701 cout << Form(
" %-9.4g +/- %-9.4g |",
covMatrix(i1, i2), covMatrixErr(i1, i2));
704 cout <<
"--------------";
705 for (
unsigned int i = 0;
i < vars.size();
i++) cout <<
"---------------------------";
717 if (!otherCell)
continue;
718 if (
cell.position().DeltaR(otherCell->
position()) > dRCut) {
delete otherCell;
continue; }
730 if (!
cell)
return true;
732 if (
layer < 0)
return true;
733 std::vector<unsigned int> allHashes;
741 if (
layer == -2) {
hashes = allHashes;
return true; }
742 for (
unsigned int h : allHashes) {
767 if (!history)
continue;
768 const Data* dataForEvent = history->data_for_event(
event);
769 if (dataForEvent)
data.push_back(
new Data(*dataForEvent));
777 std::vector<float*> floatVars;
778 std::vector<int*> intVars;
779 std::vector<std::vector<float>*> floatVects;
780 std::vector<std::vector<int>*> intVects;
781 std::map<TString, unsigned int> varIndex;
783 std::vector<TString> vars;
784 std::vector<DataFuncSet> funcs;
785 std::vector<DataFuncArgs>
args;
788 cout <<
"Making trees..." << endl;
790 TFile* flatFile = TFile::Open(
fileName +
"_tmpFlatFile.root",
"RECREATE");
791 TTree* flatTree =
new TTree(
"flatTree",
"Flat tree");
793 TFile* eventFile = TFile::Open(
fileName,
"RECREATE");
794 TTree* eventTree =
new TTree(
"eventTree",
"Event tree");
796 for (
unsigned int j = 0; j < vars.size(); j++) {
797 unsigned int index = 0;
798 if (funcs[j].isNull())
return false;
799 if (funcs[j].isInt()) {
800 int* varCont =
new int(0);
801 std::vector<int>* vectCont =
new std::vector<int>();
802 index = intVars.size();
803 intVars.push_back(varCont);
804 intVects.push_back(vectCont);
805 flatTree->Branch(vars[j], varCont);
806 eventTree->Branch(vars[j], vectCont);
809 float* varCont =
new float(0);
810 std::vector<float>* vectCont =
new std::vector<float>();
811 index = floatVars.size();
812 floatVars.push_back(varCont);
813 floatVects.push_back(vectCont);
814 flatTree->Branch(vars[j], varCont);
815 eventTree->Branch(vars[j], vectCont);
817 varIndex[vars[j]] =
index;
818 cout << vars[j] <<
" -> " <<
index << endl;
820 std::map< unsigned int, std::map< unsigned int, std::vector<long long> > > runEventIndices;
821 cout <<
"Making flat ntuple" << endl;
822 unsigned int count = 0;
826 if (
count % 100 == 0) cout <<
"Processing entry " <<
count <<
" (hash = " << iter.pos() <<
"), size = " <<
hist->nData() << endl;
827 for (
unsigned int k = 0;
k <
hist->nData();
k++) {
828 for (
unsigned int j = 0; j < vars.size(); j++) {
829 if (funcs[j].isInt())
830 *intVars[varIndex[vars[j]]] =
int(funcs[j].intVal(*
hist->data(
k),
args[j]));
832 *floatVars[varIndex[vars[j]]] = funcs[j].doubleVal(*
hist->data(
k),
args[j]);
834 runEventIndices[
hist->data(
k)->run()][
hist->data(
k)->event()].push_back(flatTree->GetEntries());
839 cout <<
"Making event tuple" << endl;
841 for (
const auto&
run : runEventIndices) {
843 cout <<
"Processing run " <<
run.first <<
" (" <<
runCount <<
" of " << runEventIndices.size() <<
")" << endl;
844 unsigned int eventCount = 0;
845 for (
const auto&
event :
run.second) {
847 if (eventCount % 1000 == 0)
848 cout <<
" processing event " <<
event.first <<
" (" << eventCount <<
" of " <<
run.second.size() <<
"), size = " <<
event.second.size() << endl;
849 for (
unsigned int j = 0; j < intVects.size(); j++) intVects[j]->
clear();
850 for (
unsigned int j = 0; j < floatVects.size(); j++) floatVects[j]->
clear();
852 flatTree->GetEntry(
index);
853 for (
unsigned int j = 0; j < vars.size(); j++) {
854 if (funcs[j].isInt())
855 intVects[varIndex[vars[j]]]->push_back(*intVars[varIndex[vars[j]]]);
857 floatVects[varIndex[vars[j]]]->push_back(*floatVars[varIndex[vars[j]]]);
864 cout <<
"Writing data..." << endl;
870 cout <<
"Cleaning up..." << endl;
876 for (
unsigned int j = 0; j < intVects.size(); j++)
delete intVects[j];
877 for (
unsigned int j = 0; j < floatVects.size(); j++)
delete floatVects[j];
878 for (
unsigned int j = 0; j < intVars.size(); j++)
delete intVars[j];
879 for (
unsigned int j = 0; j < floatVars.size(); j++)
delete floatVars[j];
881 cout <<
"Done!" << endl;