Belle II Software  release-05-01-25
eclCrystalEnergy.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2019 - 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  * - Read in payloads energy calibration payloads from specified *
11  * global tag (default ECL_energy_calibrations) or localdb, if present. *
12  * - Combine to create new ECLCrystalEnergy *
13  * - Write new ECLCrystalEnergy to localdb, with iov exp,run,-1,-1 *
14  * - Compare new version to existing version and store histograms to *
15  * ECLCrystalEnergy_exp_run.root *
16  * *
17  * Usage: eclCrystalEnergy exp run [globaltag [writeToDB]] *
18  * Set option argument writeToDB = 0 to not write output to database *
19  * *
20  **************************************************************************/
21 
22 #include <framework/database/DBImportObjPtr.h>
23 #include <framework/database/DBObjPtr.h>
24 #include <framework/database/DBStore.h>
25 #include <framework/datastore/StoreObjPtr.h>
26 #include <framework/datastore/DataStore.h>
27 #include <framework/dataobjects/EventMetaData.h>
28 #include <framework/database/Configuration.h>
29 #include <ecl/dbobjects/ECLCrystalCalib.h>
30 #include <iostream>
31 #include <TFile.h>
32 #include <TH1F.h>
33 
34 using namespace Belle2;
35 
37 namespace {
38  void setupDatabase(int exp, int run, int eventNr = 1)
39  {
41  // simulate the initialize() phase where we can register objects in the DataStore
43  evtPtr.registerInDataStore();
45  // now construct the event metadata
46  evtPtr.construct(eventNr, run, exp);
47  // and update the database contents
48  DBStore& dbstore = DBStore::Instance();
49  dbstore.update();
50  // this is only needed it the payload might be intra-run dependent,
51  // that is if it might change during one run as well
52  dbstore.updateEvent();
53  }
54 }
55 
56 
57 int main(int argc, char** argv)
58 {
59  //------------------------------------------------------------------------
60  //..Check the arguments
61  if (argc < 3 || argc > 5) {
62  std::cout << "incorrect number of arguments for eclCrystalEnergy" << std::endl;
63  return -1;
64  }
65  int experiment = std::stoi(argv[1]);
66  int run = std::stoi(argv[2]);
67  std::string gtName = "ECL_energy_calibrations";
68  if (argc >= 4) {gtName = argv[3];}
69  bool writeOutput = true;
70  if (argc == 5) {writeOutput = std::stoi(argv[4]);}
71 
72  std::cout << "eclCrystalEnergy called with arguments exp = " << experiment << " run = " << run << " GT = " << gtName <<
73  " writeOutput = " << writeOutput << std::endl;
74 
75 
76  //------------------------------------------------------------------------
79  conf.prependGlobalTag(gtName);
80  conf.prependTestingPayloadLocation("localdb/database.txt");
81 
83  DBObjPtr<Belle2::ECLCrystalCalib> GammaGamma("ECLCrystalEnergyGammaGamma");
84  DBObjPtr<Belle2::ECLCrystalCalib> Existing("ECLCrystalEnergy");
85 
87  setupDatabase(experiment, run);
88 
89  //------------------------------------------------------------------------
91  std::vector<float> GammaGammaCalib;
92  std::vector<float> GammaGammaCalibUnc;
93  GammaGammaCalib = GammaGamma->getCalibVector();
94  GammaGammaCalibUnc = GammaGamma->getCalibUncVector();
95 
96  std::vector<float> ExistingCalib;
97  std::vector<float> ExistingCalibUnc;
98  ExistingCalib = Existing->getCalibVector();
99  ExistingCalibUnc = Existing->getCalibUncVector();
100 
101  //------------------------------------------------------------------------
103  //..For now, use Gamma Gamma if available; otherwise, use existing value
104  std::vector<float> NewCalib;
105  std::vector<float> NewCalibUnc;
106  NewCalib.resize(8736);
107  NewCalibUnc.resize(8736);
108  for (int ic = 0; ic < 8736; ic++) {
109  if (GammaGammaCalib[ic] > 0.) {
110  NewCalib[ic] = GammaGammaCalib[ic];
111  NewCalibUnc[ic] = GammaGammaCalibUnc[ic];
112  } else {
113  NewCalib[ic] = ExistingCalib[ic];
114  NewCalibUnc[ic] = ExistingCalibUnc[ic];
115  }
116  }
117 
118  //------------------------------------------------------------------------
120  for (int ic = 0; ic < 9000; ic += 1000) {
121  std::cout << ic + 1 << " " << ExistingCalib[ic] << " " << ExistingCalibUnc[ic] << " " << GammaGammaCalib[ic] << " " <<
122  GammaGammaCalibUnc[ic] << " " << NewCalib[ic] << " " << NewCalibUnc[ic] << std::endl;
123  }
124 
125  //------------------------------------------------------------------------
126  //..Histograms of existing calibration, new calibration, and ratio new/old
127  TString payloadTitle = "ECLCrystalEnergy";
128  payloadTitle += "_";
129  payloadTitle += experiment;
130  payloadTitle += "_";
131  payloadTitle += run;
132  TString fname = payloadTitle;
133  fname += ".root";
134  TFile hfile(fname, "recreate");
135 
136  TString htitle = payloadTitle;
137  htitle += " existing values;cellID";
138  TH1F* existingPayload = new TH1F("existingPayload", htitle, 8736, 1, 8737);
139 
140  htitle = payloadTitle;
141  htitle += " new values;cellID";
142  TH1F* newPayload = new TH1F("newPayload", htitle, 8736, 1, 8737);
143 
144  htitle = payloadTitle;
145  htitle += " ratio new/old;cellID";
146  TH1F* payloadRatioVsCellID = new TH1F("payloadRatioVsCellID", htitle, 8736, 1, 8737);
147 
148  htitle = payloadTitle;
149  htitle += " ratio new/old";
150  TH1F* payloadRatio = new TH1F("payloadRatio", htitle, 200, 0.95, 1.05);
151 
152  for (int cellID = 1; cellID <= 8736; cellID++) {
153  existingPayload->SetBinContent(cellID, ExistingCalib[cellID - 1]);
154  existingPayload->SetBinError(cellID, ExistingCalibUnc[cellID - 1]);
155 
156  newPayload->SetBinContent(cellID, NewCalib[cellID - 1]);
157  newPayload->SetBinError(cellID, NewCalibUnc[cellID - 1]);
158 
159  float ratio = 1.;
160  float ratioUnc = 0.;
161  if (abs(ExistingCalib[cellID - 1]) > 1.0e-12) {
162  ratio = NewCalib[cellID - 1] / ExistingCalib[cellID - 1];
163  float rUnc0 = ExistingCalibUnc[cellID - 1] / ExistingCalib[cellID - 1];
164  float rUnc1 = 0.;
165  if (abs(NewCalib[cellID - 1]) > 1.0e-12) {rUnc1 = NewCalibUnc[cellID - 1] / NewCalib[cellID - 1];}
166  ratioUnc = ratio * sqrt(rUnc0 * rUnc0 + rUnc1 * rUnc1);
167  }
168 
169  payloadRatioVsCellID->SetBinContent(cellID, ratio);
170  payloadRatioVsCellID->SetBinError(cellID, ratioUnc);
171 
172  payloadRatio->Fill(ratio);
173  }
174 
175  hfile.cd();
176  hfile.Write();
177  hfile.Close();
178  std::cout << std::endl << "Values written to " << fname << std::endl;
179 
180  //------------------------------------------------------------------------
181  //..Write new payload to localdb if requested
182  if (writeOutput) {
183  std::cout << "Creating importer" << std::endl;
184  Belle2::DBImportObjPtr<Belle2::ECLCrystalCalib> importer("ECLCrystalEnergy");
185  importer.construct();
186  importer->setCalibVector(NewCalib, NewCalibUnc);
187  importer.import(Belle2::IntervalOfValidity(experiment, run, -1, -1));
188  std::cout << "Successfully wrote payload ECLCrystalEnergy with iov " << experiment << "," << run << ",-1,-1" << std::endl;
189  } else {
190  std::cout << "Payload not written to database" << std::endl;
191  }
192 }
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 >
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::DBStore::update
void update()
Updates all objects that are outside their interval of validity.
Definition: DBStore.cc:87