\ #include #include #include #include "TFile.h" #include "TH1D.h" #include "TH2D.h" #include "TF1.h" #include "TMath.h" #include "JDetector/JDetector.hh" #include "JDetector/JDetectorToolkit.hh" #include "JDetector/JStringRouter.hh" #include "Jeep/JPrint.hh" #include "Jeep/JParser.hh" /** * \file * * Auxiliary program to extract string-string correlation information from the output created with JMonitorL1dt * * \author dsamtleben */ int main(int argc, char **argv) { using namespace std; using namespace JPP; string inputFile; string outputFile; string detectorFile; int neighbour, setmax; int debug; try { JParser<> zap("Program to extract time offsets of DOM-DOM correlations"); zap['f'] = make_field(inputFile, "input file") = "monitor.root"; zap['o'] = make_field(outputFile, "output file") = "corr_max.root"; zap['a'] = make_field(detectorFile, "detector file"); zap['n'] = make_field(neighbour, "neighbour level") = 5; zap['m'] = make_field(setmax, "minimal entries") = 5; zap['d'] = make_field(debug) = 0; if (zap.read(argc, argv) != 0) { return 1; } } catch(const exception &error) { FATAL(error.what() << endl); } JDetector detector; try { load(detectorFile, detector); } catch(const JException& error) { FATAL(error); } const JStringRouter string(detector); map > zmap; for (JDetector::const_iterator module = detector.begin(); module != detector.end(); ++module) { zmap[module->getString()][module->getFloor()] = module->getID(); } struct h1_t { h1_t() : p(NULL) // pointer must be initialised {} TH1D* p; }; TF1 f1("f1", "[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2])) + [3]"); const int number_of_strings = getNumberOfStrings(detector); // cout << "number of strings " << number_of_strings << endl; double maxarr[number_of_strings][number_of_strings][3]; int duarr[number_of_strings]; map > H1; TFile* in = TFile::Open(inputFile.c_str(), "exist"); if (in == NULL || !in->IsOpen()) { FATAL("File: " << inputFile << " not opened." << endl); } for (const auto& string_1 : zmap) { duarr[string.getIndex(string_1.first)]=string_1.first; for (const auto& floor_1 : string_1.second) { const int module_1 = floor_1.second; // parent module TH2D* h2 = (TH2D*) in->Get(MAKE_CSTRING(module_1 << ".2S")); if (h2 == NULL) { continue; } for (const auto& string_2 : zmap) { duarr[string.getIndex(string_2.first)]=string_2.first; if (string_1.first != string_2.first) { TH1D* h1 = H1[string_1.first][string_2.first].p; if (h1 == NULL) { h1 = new TH1D(MAKE_CSTRING(string_1.first << "_" << string_2.first << "_" << neighbour << ".2T"), NULL, h2->GetYaxis()->GetNbins(), h2->GetYaxis()->GetXmin(), h2->GetYaxis()->GetXmax()); H1[string_1.first][string_2.first].p = h1; // book keeping } for (const auto& floor_2 : string_2.second) { if (floor_1.first > floor_2.first && (floor_1.first - floor_2.first) == neighbour) { // floor condition const int module_2 = floor_2.second; // daughter module TH1D* py = h2->ProjectionY("__py", h2->GetXaxis()->FindBin(TString(Form("%i", module_2))), h2->GetXaxis()->FindBin(TString(Form("%i", module_2))), "e"); h1->Add(py); // add data delete py; } } //------------ if (floor_1.first == 18){ double mm=h1->GetXaxis()->GetBinCenter(h1->GetMaximumBin()); f1.SetParameter(0, double(h1->GetMaximum())); f1.SetParameter(1, double(mm)); f1.SetParameter(2, 50.); // fit if (h1->GetMaximum()>0){ // cout << "we are fitting " << string_1.first << " " << string_2.first << endl; h1->Fit(&f1,"LQ", "same"); // maxarr[ string.getIndex(string_1.first) ][ string.getIndex(string_2.first) ][0] = h1->GetXaxis()->GetBinCenter(h1->GetMaximumBin()); // maxarr[ string.getIndex(string_1.first) ][ string.getIndex(string_2.first) ][1] = h1->GetMaximum(); maxarr[ string.getIndex(string_1.first) ][ string.getIndex(string_2.first) ][0] = h1->GetFunction("f1")->GetParameter(1); maxarr[ string.getIndex(string_1.first) ][ string.getIndex(string_2.first) ][1] = h1->GetFunction("f1")->GetParameter(0); maxarr[ string.getIndex(string_1.first) ][ string.getIndex(string_2.first) ][2] = h1->GetFunction("f1")->GetParameter(2); } } } // if string 1 != string 2 } // loop over string 2 } // loop over modules of string 1 } // loop over string 1 for (int i=0;isetmax && maxarr[j][i][1]>setmax && maxarr[i][j][2]>20. && maxarr[j][i][2]>20.){ cout << i << " " << j << " " << maxarr[i][j][0] << " " << maxarr[j][i][0] << " " << (maxarr[i][j][0]-maxarr[j][i][0])/2. << " " << maxarr[i][j][1] << " " << maxarr[j][i][1] << " " << " " << duarr[i] << " " << duarr[j] << " " << neighbour << endl; } } } TFile out(outputFile.c_str(), "recreate"); out.cd(); for (auto& h1 : H1) { for (auto& h2 : h1.second) { h2.second.p->Write(); } } out.Close(); }