35 #include <framework/database/DBImportObjPtr.h>
36 #include <framework/database/DBObjPtr.h>
37 #include <framework/database/DBStore.h>
38 #include <framework/datastore/StoreObjPtr.h>
39 #include <framework/datastore/DataStore.h>
40 #include <framework/dataobjects/EventMetaData.h>
41 #include <framework/database/Configuration.h>
42 #include <framework/logging/LogSystem.h>
43 #include <ecl/dbobjects/ECLCrystalCalib.h>
53 void setupDatabase(
int exp,
int run,
int eventNr = 1)
58 evtPtr.registerInDataStore();
60 std::cout <<
"about to construct EventMetaData, exp = " << exp <<
" run = " << run <<
" eventNr = " << eventNr << std::endl;
69 int main(
int argc,
char** argv)
71 if (argc < 4 || argc > 5) {
72 std::cout <<
"insufficient arguments for eclElectronicsPayloads" << std::endl;
75 std::string payloadName = argv[1];
76 if (payloadName !=
"ECLCrystalElectronics" and payloadName !=
"ECLCrystalElectronicsTime" and payloadName !=
"ECLRefAmplNom"
77 and payloadName !=
"ECLRefTimeNom") {
78 std::cout <<
"First argument must be ECLCrystalElectronics, ECLCrystalElectronicsTime, ECLRefAmplNom, or ECLRefTimeNom" <<
82 int experiment = std::stoi(argv[2]);
83 int run = std::stoi(argv[3]);
84 bool writeOutput =
true;
85 if (argc == 5) writeOutput = std::stoi(argv[4]);
86 std::cout <<
"eclElectronicsPayloads called with arguments " << payloadName <<
" " << experiment <<
" " << run <<
" " << writeOutput
92 conf.prependGlobalTag(
"ECL_localrun_data");
93 conf.prependTestingPayloadLocation(
"localdb/database.txt");
98 logging->setDebugLevel(10);
101 std::cout <<
"calling setupDatabase " << std::endl;
102 setupDatabase(experiment, run);
115 std::cout <<
"Reading ECLRefAmpl, ECLRefAmplNom, ECLRefTime, ECLRefTimeNom, ECLCrystalElectronics, and ECLCrystalElectronicsTime" <<
117 std::cout <<
"Dumping " << payloadName << std::endl;
118 existingObject->Dump();
121 std::vector<float> currentValues = existingObject->getCalibVector();
122 std::vector<float> currentUnc = existingObject->getCalibUncVector();
124 std::vector<float> refAmpl = InputAmpl->getCalibVector();
125 std::vector<float> refAmplUnc = InputAmpl->getCalibUncVector();
127 std::vector<float> refAmplNom = InputAmplNom->getCalibVector();
129 std::vector<float> refTime = InputTime->getCalibVector();
130 std::vector<float> refTimeUnc = InputTime->getCalibUncVector();
132 std::vector<float> refTimeNom = InputTimeNom->getCalibVector();
134 std::vector<float> crysElec = CurrentElec->getCalibVector();
136 std::vector<float> crysTime = CurrentTime->getCalibVector();
140 std::cout << std::endl <<
"Reference amplitudes and times read from database " << std::endl;
141 for (
int ic = 0; ic < 9000; ic += 1000) {
142 std::cout <<
"cellID " << ic + 1 <<
" ref amplitude = " << refAmpl[ic] <<
" +/- " << refAmplUnc[ic] <<
" nom = " << refAmplNom[ic]
143 <<
" ref time = " << refTime[ic]
144 <<
" +/- " << refTimeUnc[ic] <<
" nom = " << refTimeNom[ic] << std::endl;
149 std::vector<float> newValues;
150 std::vector<float> newUnc;
151 for (
int ic = 0; ic < 8736; ic++) {
152 if (payloadName ==
"ECLCrystalElectronics") {
153 newValues.push_back(refAmplNom[ic] / refAmpl[ic]);
154 newUnc.push_back(newValues[ic]*refAmplUnc[ic] / refAmpl[ic]);
155 }
else if (payloadName ==
"ECLCrystalElectronicsTime") {
156 newValues.push_back(refTime[ic] - refTimeNom[ic]);
157 newUnc.push_back(refTimeUnc[ic]);
158 }
else if (payloadName ==
"ECLRefAmplNom") {
159 newValues.push_back(crysElec[ic]*refAmpl[ic]);
160 newUnc.push_back(0.);
161 }
else if (payloadName ==
"ECLRefTimeNom") {
162 newValues.push_back(refTime[ic] - crysTime[ic]);
163 newUnc.push_back(0.);
170 std::cout << std::endl <<
"Comparison of existing and new values for " << payloadName << std::endl;
171 for (
int ic = 0; ic < 9000; ic += 1000) {
172 std::cout <<
"cellID " << ic + 1 <<
" existing = " << currentValues[ic] <<
" +/- " << currentUnc[ic] <<
" new = " << newValues[ic]
173 <<
" +/- " << newUnc[ic] << std::endl;
175 std::cout << std::endl;
177 TString payloadTitle = payloadName;
179 payloadTitle += experiment;
182 TString fname = payloadTitle;
184 TFile hfile(fname,
"recreate");
185 TString htitle = payloadTitle;
186 htitle +=
" existing calibration values;cellID";
187 TH1F* existingCalib =
new TH1F(
"existingCalib", htitle, 8736, 1, 8737);
189 htitle = payloadTitle;
190 htitle +=
" new calibration values;cellID";
191 TH1F* newCalib =
new TH1F(
"newCalib", htitle, 8736, 1, 8737);
193 htitle = payloadTitle;
195 TH1F* calibRatio =
new TH1F(
"calibRatio", htitle, 200, 0.9, 1.1);
197 htitle = payloadTitle;
198 htitle +=
" difference";
199 TH1F* calibDiff =
new TH1F(
"calibDiff", htitle, 200, -100, 100);
201 htitle = payloadTitle;
202 htitle +=
" reference";
203 TH1F* refValues =
new TH1F(
"refValues", htitle, 8736, 1, 8737);
205 htitle = payloadTitle;
206 htitle +=
" ratio vs cellID;cellID;new/old";
207 TH1F* ratioVsCellID =
new TH1F(
"ratioVsCellID", htitle, 8736, 1, 8737);
209 htitle = payloadTitle;
210 htitle +=
" diff vs cellID;cellID;new - old";
211 TH1F* diffVsCellID =
new TH1F(
"diffVsCellID", htitle, 8736, 1, 8737);
213 for (
int cellID = 1; cellID <= 8736; cellID++) {
214 float oldValue = currentValues[cellID - 1];
215 float newValue = newValues[cellID - 1];
217 if (oldValue != 0.) {
218 ratio = newValue / oldValue;
219 }
else if (newValue != 0.) {
220 ratio = ratio * newValue / fabs(newValue);
223 existingCalib->SetBinContent(cellID, oldValue);
224 existingCalib->SetBinError(cellID, currentUnc[cellID - 1]);
225 newCalib->SetBinContent(cellID, newValue);
226 newCalib->SetBinError(cellID, newUnc[cellID - 1]);
227 calibRatio->Fill(ratio);
228 ratioVsCellID->SetBinContent(cellID, ratio);
229 ratioVsCellID->SetBinError(cellID, 0);
230 calibDiff->Fill(newValue - oldValue);
231 diffVsCellID->SetBinContent(cellID, newValue - oldValue);
232 diffVsCellID->SetBinError(cellID, 0);
233 if (payloadName ==
"ECLCrystalElectronics" or payloadName ==
"ECLRefAmplNom") {
234 refValues->SetBinContent(cellID, refAmpl[cellID - 1]);
235 refValues->SetBinError(cellID, refAmplUnc[cellID - 1]);
237 refValues->SetBinContent(cellID, refTime[cellID - 1]);
238 refValues->SetBinError(cellID, refTimeUnc[cellID - 1]);
242 if ((payloadName ==
"ECLCrystalElectronics" or payloadName ==
"ECLRefAmplNom") and (ratio<0.99 or ratio>1.01)) {
243 std::cout <<
"Ratio = " << ratio <<
" for cellID = " << cellID <<
" refAmpl = " << refAmpl[cellID - 1] <<
" refAmplNom = " <<
244 refAmplNom[cellID - 1] << std::endl;
245 }
else if (abs(newValue - oldValue) > 20.) {
246 std::cout <<
"Difference = " << newValue - oldValue <<
" for cellID = " << cellID <<
" refTime = " << refTime[cellID - 1] <<
247 " refTimeNom = " << refTimeNom[cellID - 1] << std::endl;
254 std::cout << std::endl <<
"Comparison of existing and new calibration values written to " << fname << std::endl;
259 std::cout <<
"Creating importer" << std::endl;
261 importer.construct();
262 importer->setCalibVector(newValues, newUnc);
264 std::cout <<
"Successfully wrote payload " << payloadName <<
" with iov " << experiment <<
"," << run <<
",-1,-1" << std::endl;