Belle II Software  release-05-01-25
eclElectronicsPayloads.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2018 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Christopher Hearty hearty@physics.ubc.ca *
7  * *
8  * This software is provided "as is" without any warranty. *
9  * *
10  * Standard usage is to read payloads ECLRefAmpl, ECLRefAmplNom, *
11  * ECLRefTime and ECLRefTimeNom and use them to derive payloads *
12  * ECLCrystalElectronics or ECLCrystalElectronicsTime *
13  * *
14  * Alternatively, find new payloads ECLRefAmplNom or ECLRefTimeNom to *
15  * keep ECLCrystalElectronics or ECLCrystalElectronicsTime constant *
16  * *
17  * Also performs a comparison of new and existing calibration values and *
18  * writes these to a root file. *
19  * *
20  * Payloads are read from localdb if present, otherwise from *
21  * ECL_localrun_data *
22  * They are written to localdb with iov = exp,run,-1,-1 *
23  * *
24  * Usage: *
25  * eclElectronicsPayloads payloadName exp run [writeToDB] *
26  * where payloadName = ECLCrystalElectronics, ECLCrystalElectronicsName, *
27  * ECLRefAmplNom, or ECLRefTimeNom *
28  * *
29  * exp and run specify the start of the iov, and are used to read *
30  * the reference amplitudes and times *
31  * Option argument writeToDB = 0 to not write output to database *
32  * *
33  **************************************************************************/
34 
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>
44 #include <iostream>
45 #include <TFile.h>
46 #include <TH1F.h>
47 
48 using namespace Belle2;
49 
50 //------------------------------------------------------------------------
51 //..Set experiment, run, and event numbers before reading a payload from the DB
52 namespace {
53  void setupDatabase(int exp, int run, int eventNr = 1)
54  {
56  // simulate the initialize() phase where we can register objects in the DataStore
58  evtPtr.registerInDataStore();
60  std::cout << "about to construct EventMetaData, exp = " << exp << " run = " << run << " eventNr = " << eventNr << std::endl;
61  evtPtr.construct(eventNr, run, exp);
62  DBStore& dbstore = DBStore::Instance();
63  dbstore.update();
64  dbstore.updateEvent();
65  }
66 }
67 
68 //------------------------------------------------------------------------
69 int main(int argc, char** argv)
70 {
71  if (argc < 4 || argc > 5) {
72  std::cout << "insufficient arguments for eclElectronicsPayloads" << std::endl;
73  return -1;
74  }
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" <<
79  std::endl;
80  return -1;
81  }
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
87  << std::endl;
88 
89  //------------------------------------------------------------------------
90  //..Specify database
92  conf.prependGlobalTag("ECL_localrun_data");
93  conf.prependTestingPayloadLocation("localdb/database.txt");
94 
95  //..set debug level
97  logging->setLogLevel(LogConfig::c_Debug);
98  logging->setDebugLevel(10);
99 
100  //..Populate database contents
101  std::cout << "calling setupDatabase " << std::endl;
102  setupDatabase(experiment, run);
103 
104  //------------------------------------------------------------------------
105  //..Read input and existing output payloads from database
106  DBObjPtr<Belle2::ECLCrystalCalib> existingObject(payloadName);
107  DBObjPtr<Belle2::ECLCrystalCalib> InputAmpl("ECLRefAmpl");
108  DBObjPtr<Belle2::ECLCrystalCalib> InputAmplNom("ECLRefAmplNom");
109  DBObjPtr<Belle2::ECLCrystalCalib> InputTime("ECLRefTime");
110  DBObjPtr<Belle2::ECLCrystalCalib> InputTimeNom("ECLRefTimeNom");
111  DBObjPtr<Belle2::ECLCrystalCalib> CurrentElec("ECLCrystalElectronics");
112  DBObjPtr<Belle2::ECLCrystalCalib> CurrentTime("ECLCrystalElectronicsTime");
113 
114  //..Print out some information about the existing payload
115  std::cout << "Reading ECLRefAmpl, ECLRefAmplNom, ECLRefTime, ECLRefTimeNom, ECLCrystalElectronics, and ECLCrystalElectronicsTime" <<
116  std::endl;
117  std::cout << "Dumping " << payloadName << std::endl;
118  existingObject->Dump();
119 
120  //..Get vectors of values from the payloads
121  std::vector<float> currentValues = existingObject->getCalibVector();
122  std::vector<float> currentUnc = existingObject->getCalibUncVector();
123 
124  std::vector<float> refAmpl = InputAmpl->getCalibVector();
125  std::vector<float> refAmplUnc = InputAmpl->getCalibUncVector();
126 
127  std::vector<float> refAmplNom = InputAmplNom->getCalibVector();
128 
129  std::vector<float> refTime = InputTime->getCalibVector();
130  std::vector<float> refTimeUnc = InputTime->getCalibUncVector();
131 
132  std::vector<float> refTimeNom = InputTimeNom->getCalibVector();
133 
134  std::vector<float> crysElec = CurrentElec->getCalibVector();
135 
136  std::vector<float> crysTime = CurrentTime->getCalibVector();
137 
138 
139  //..Print out a few values for quality control
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;
145  }
146 
147  //------------------------------------------------------------------------
148  //..Calculate the new values for requested payload
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.);
164  }
165 
166  }
167 
168  //------------------------------------------------------------------------
169  //..Compare current values to new ones
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;
174  }
175  std::cout << std::endl;
176 
177  TString payloadTitle = payloadName;
178  payloadTitle += "_";
179  payloadTitle += experiment;
180  payloadTitle += "_";
181  payloadTitle += run;
182  TString fname = payloadTitle;
183  fname += ".root";
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);
188 
189  htitle = payloadTitle;
190  htitle += " new calibration values;cellID";
191  TH1F* newCalib = new TH1F("newCalib", htitle, 8736, 1, 8737);
192 
193  htitle = payloadTitle;
194  htitle += " ratio";
195  TH1F* calibRatio = new TH1F("calibRatio", htitle, 200, 0.9, 1.1);
196 
197  htitle = payloadTitle;
198  htitle += " difference";
199  TH1F* calibDiff = new TH1F("calibDiff", htitle, 200, -100, 100);
200 
201  htitle = payloadTitle;
202  htitle += " reference";
203  TH1F* refValues = new TH1F("refValues", htitle, 8736, 1, 8737);
204 
205  htitle = payloadTitle;
206  htitle += " ratio vs cellID;cellID;new/old";
207  TH1F* ratioVsCellID = new TH1F("ratioVsCellID", htitle, 8736, 1, 8737);
208 
209  htitle = payloadTitle;
210  htitle += " diff vs cellID;cellID;new - old";
211  TH1F* diffVsCellID = new TH1F("diffVsCellID", htitle, 8736, 1, 8737);
212 
213  for (int cellID = 1; cellID <= 8736; cellID++) {
214  float oldValue = currentValues[cellID - 1];
215  float newValue = newValues[cellID - 1];
216  float ratio = 9999.;
217  if (oldValue != 0.) {
218  ratio = newValue / oldValue;
219  } else if (newValue != 0.) {
220  ratio = ratio * newValue / fabs(newValue);
221  }
222 
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]);
236  } else {
237  refValues->SetBinContent(cellID, refTime[cellID - 1]);
238  refValues->SetBinError(cellID, refTimeUnc[cellID - 1]);
239  }
240 
241  //..Note any large changes
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;
248  }
249  }
250 
251  hfile.cd();
252  hfile.Write();
253  hfile.Close();
254  std::cout << std::endl << "Comparison of existing and new calibration values written to " << fname << std::endl;
255 
256  //------------------------------------------------------------------------
257  //..Write out to localdb if requested
258  if (writeOutput) {
259  std::cout << "Creating importer" << std::endl;
261  importer.construct();
262  importer->setCalibVector(newValues, newUnc);
263  importer.import(Belle2::IntervalOfValidity(experiment, run, -1, -1));
264  std::cout << "Successfully wrote payload " << payloadName << " with iov " << experiment << "," << run << ",-1,-1" << std::endl;
265  }
266 }
Belle2::IntervalOfValidity
A class that describes the interval of experiments/runs for which an object in the database is valid.
Definition: IntervalOfValidity.h:35
Belle2::DataStore::Instance
static DataStore & Instance()
Instance of singleton Store.
Definition: DataStore.cc:54
Belle2::DataStore::setInitializeActive
void setInitializeActive(bool active)
Setter for m_initializeActive.
Definition: DataStore.cc:94
Belle2::DBObjPtr< Belle2::ECLCrystalCalib >
Belle2::LogSystem::getLogConfig
LogConfig * getLogConfig()
Returns global log system configuration.
Definition: LogSystem.h:88
main
int main(int argc, char **argv)
Run all tests.
Definition: test_main.cc:77
Belle2::DBStore
Singleton class to cache database objects.
Definition: DBStore.h:42
Belle2::StoreObjPtr::construct
bool construct(Args &&... params)
Construct an object of type T in this StoreObjPtr, using the provided constructor arguments.
Definition: StoreObjPtr.h:128
Belle2::Conditions::Configuration::getInstance
static Configuration & getInstance()
Get a reference to the instance which will be used when the Database is initialized.
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::StoreObjPtr
Type-safe access to single objects in the data store.
Definition: ParticleList.h:33
Belle2::DBImportObjPtr
Class for importing a single object to the database.
Definition: DBImportObjPtr.h:33
Belle2::DBStore::Instance
static DBStore & Instance()
Instance of a singleton DBStore.
Definition: DBStore.cc:36
Belle2::DBStore::updateEvent
void updateEvent()
Updates all intra-run dependent objects.
Definition: DBStore.cc:150
Belle2::LogSystem::Instance
static LogSystem & Instance()
Static method to get a reference to the LogSystem instance.
Definition: LogSystem.cc:33
Belle2::LogConfig::c_Debug
@ c_Debug
Debug: for code development.
Definition: LogConfig.h:36
Belle2::LogConfig
The LogConfig class.
Definition: LogConfig.h:32
Belle2::DBStore::update
void update()
Updates all objects that are outside their interval of validity.
Definition: DBStore.cc:87