Belle II Software  release-05-02-19
eclCosmicEAlgorithm Class Reference

class eclCosmiEAlgorithm. More...

#include <eclCosmicEAlgorithm.h>

Inheritance diagram for eclCosmicEAlgorithm:
Collaboration diagram for eclCosmicEAlgorithm:

Public Types

enum  EResult {
  c_OK,
  c_Iterate,
  c_NotEnoughData,
  c_Failure,
  c_Undefined
}
 The result of calibration. More...
 

Public Member Functions

 eclCosmicEAlgorithm ()
 Constructor.
 
virtual ~eclCosmicEAlgorithm ()
 Destructor.
 
std::string getPrefix () const
 Get the prefix used for getting calibration data.
 
bool checkPyExpRun (PyObject *pyObj)
 Checks that a PyObject can be successfully converted to an ExpRun type. More...
 
Calibration::ExpRun convertPyExpRun (PyObject *pyObj)
 Performs the conversion of PyObject to ExpRun. More...
 
std::string getCollectorName () const
 Alias for prefix. More...
 
void setPrefix (const std::string &prefix)
 Set the prefix used to identify datastore objects.
 
void setInputFileNames (PyObject *inputFileNames)
 Set the input file names used for this algorithm from a Python list. More...
 
PyObject * getInputFileNames ()
 Get the input file names used for this algorithm and pass them out as a Python list of unicode strings.
 
std::vector< Calibration::ExpRun > getRunListFromAllData () const
 Get the complete list of runs from inspection of collected data.
 
RunRange getRunRangeFromAllData () const
 Get the complete RunRange from inspection of collected data.
 
IntervalOfValidity getIovFromAllData () const
 Get the complete IoV from inspection of collected data.
 
void fillRunToInputFilesMap ()
 Fill the mapping of ExpRun -> Files.
 
std::string getGranularity () const
 Get the granularity of collected data.
 
EResult execute (std::vector< Calibration::ExpRun > runs={}, int iteration=0, IntervalOfValidity iov=IntervalOfValidity())
 Runs calibration over vector of runs for a given iteration. More...
 
EResult execute (PyObject *runs, int iteration=0, IntervalOfValidity iov=IntervalOfValidity())
 Runs calibration over Python list of runs. Converts to C++ and then calls the other execute() function.
 
std::list< Database::DBImportQuery > & getPayloads ()
 Get constants (in TObjects) for database update from last execution.
 
std::list< Database::DBImportQuerygetPayloadValues ()
 Get constants (in TObjects) for database update from last execution but passed by VALUE.
 
bool commit ()
 Submit constants from last calibration into database.
 
bool commit (std::list< Database::DBImportQuery > payloads)
 Submit constants from a (potentially previous) set of payloads.
 
const std::string & getDescription () const
 Get the description of the algoithm (set by developers in constructor)
 
bool loadInputJson (const std::string &jsonString)
 Load the m_inputJson variable from a string (useful from Python interface). The rturn bool indicates success or failure.
 
const std::string dumpOutputJson () const
 Dump the JSON string of the output JSON object.
 
const std::vector< Calibration::ExpRun > findPayloadBoundaries (std::vector< Calibration::ExpRun > runs, int iteration=0)
 Used to discover the ExpRun boundaries that you want the Python CAF to execute on. This is optional and only used in some.
 
template<>
std::shared_ptr< TTree > getObjectPtr (const std::string &name, const std::vector< Calibration::ExpRun > &requestedRuns)
 Specialization of getObjectPtr<TTree>.
 

Public Attributes

int cellIDLo
 Parameters to control Novosibirsk fit to signal measured in each crystal. More...
 
int cellIDHi
 Last cellID to be fit.
 
int minEntries
 All crystals to be fit must have at least minEntries events in the fit range.
 
int maxIterations
 no more than maxIteration iterations
 
double tRatioMin
 Fit range is adjusted so that fit at upper endpoint is between tRatioMin and tRatioMax of peak.
 
double tRatioMax
 Fit range is adjusted so that fit at upper endpoint is between tRatioMin and tRatioMax of peak.
 
bool performFits
 if false, input histograms are copied to output, but no fits are done.
 
bool findExpValues
 if true, fits are used to find expected energy deposit for each crystal instead of the calibration constant
 
int storeConst
 controls which values are written to the database. More...
 

Protected Member Functions

virtual EResult calibrate () override
 Run algorithm on events. More...
 
void setInputFileNames (std::vector< std::string > inputFileNames)
 Set the input file names used for this algorithm. More...
 
virtual bool isBoundaryRequired (const Calibration::ExpRun &)
 Given the current collector data, make a decision about whether or not this run should be the start of a payload boundary.
 
virtual void boundaryFindingSetup (std::vector< Calibration::ExpRun >, int)
 If you need to make some changes to your algorithm class before 'findPayloadBoundaries' is run, make them in this function.
 
virtual void boundaryFindingTearDown ()
 Put your algorithm back into a state ready for normal execution if you need to.
 
const std::vector< Calibration::ExpRun > & getRunList () const
 Get the list of runs for which calibration is called.
 
int getIteration () const
 Get current iteration.
 
std::vector< std::string > getVecInputFileNames () const
 Get the input file names used for this algorithm as a STL vector.
 
template<class T >
std::shared_ptr< T > getObjectPtr (const std::string &name, const std::vector< Calibration::ExpRun > &requestedRuns)
 Get calibration data object by name and list of runs, the Merge function will be called to generate the overall object.
 
template<class T >
std::shared_ptr< T > getObjectPtr (std::string name)
 Get calibration data object (for all runs the calibration is requested for) This function will only work during or after execute() has been called once.
 
template<>
shared_ptr< TTree > getObjectPtr (const string &name, const vector< ExpRun > &requestedRuns)
 We cheekily cast the TChain to TTree for the returned pointer so that the user never knows Hopefully this doesn't cause issues if people do low level stuff to the tree...
 
std::string getGranularityFromData () const
 Get the granularity of collected data.
 
void saveCalibration (TClonesArray *data, const std::string &name)
 Store DBArray payload with given name with default IOV.
 
void saveCalibration (TClonesArray *data, const std::string &name, const IntervalOfValidity &iov)
 Store DBArray with given name and custom IOV.
 
void saveCalibration (TObject *data)
 Store DB payload with default name and default IOV.
 
void saveCalibration (TObject *data, const IntervalOfValidity &iov)
 Store DB payload with default name and custom IOV.
 
void saveCalibration (TObject *data, const std::string &name)
 Store DB payload with given name with default IOV.
 
void saveCalibration (TObject *data, const std::string &name, const IntervalOfValidity &iov)
 Store DB payload with given name and custom IOV.
 
void updateDBObjPtrs (const unsigned int event, const int run, const int experiment)
 Updates any DBObjPtrs by calling update(event) for DBStore.
 
void setDescription (const std::string &description)
 Set algorithm description (in constructor)
 
void clearCalibrationData ()
 Clear calibration data.
 
Calibration::ExpRun getAllGranularityExpRun () const
 Returns the Exp,Run pair that means 'Everything'. Currently unused.
 
void resetInputJson ()
 Clears the m_inputJson member variable.
 
void resetOutputJson ()
 Clears the m_outputJson member variable.
 
template<class T >
void setOutputJsonValue (const std::string &key, const T &value)
 Set a key:value pair for the outputJson object, expected to used interally during calibrate()
 
template<class T >
const T getOutputJsonValue (const std::string &key) const
 Get a value using a key from the JSON output object, not sure why you would want to do this.
 
template<class T >
const T getInputJsonValue (const std::string &key) const
 Get an input JSON value using a key. The normal exceptions are raised when the key doesn't exist.
 
const nlohmann::json & getInputJsonObject () const
 Get the entire top level JSON object. We explicitly say this must be of object type so that we might pick.
 
bool inputJsonKeyExists (const std::string &key) const
 Test for a key in the input JSON object.
 

Protected Attributes

std::vector< Calibration::ExpRun > m_boundaries
 

Private Member Functions

std::string getExpRunString (Calibration::ExpRun &expRun) const
 Gets the "exp.run" string repr. of (exp,run)
 
std::string getFullObjectPath (std::string name, Calibration::ExpRun expRun) const
 constructs the full TDirectory + Key name of an object in a TFile based on its name and exprun
 

Private Attributes

int fitOK = 16
 fit is OK
 
int iterations = 8
 fit reached max number of iterations, but is useable
 
int atLimit = 4
 a parameter is at the limit; fit not useable
 
int poorFit = 3
 low chi square; fit not useable
 
int noPeak = 2
 Novosibirsk component of fit is negligible; fit not useable.
 
int notFit = -1
 no fit performed
 
std::vector< std::string > m_inputFileNames
 List of input files to the Algorithm, will initially be user defined but then gets the wildcards expanded during execute()
 
std::map< Calibration::ExpRun, std::vector< std::string > > m_runsToInputFiles
 Map of Runs to input files. Gets filled when you call getRunRangeFromAllData, gets cleared when setting input files again.
 
std::string m_granularityOfData
 Granularity of input data. This only changes when the input files change so it isn't specific to an execution.
 
ExecutionData m_data
 Data specific to a SINGLE execution of the algorithm. Gets reset at the beginning of execution.
 
std::string m_description {""}
 Description of the algorithm.
 
std::string m_prefix {""}
 The name of the TDirectory the collector objects are contained within.
 
nlohmann::json m_jsonExecutionInput = nlohmann::json::object()
 Optional input JSON object used to make decisions about how to execute the algorithm code.
 
nlohmann::json m_jsonExecutionOutput = nlohmann::json::object()
 Optional output JSON object that can be set during the execution by the underlying algorithm code.
 

Static Private Attributes

static const Calibration::ExpRun m_allExpRun = make_pair(-1, -1)
 

Detailed Description

class eclCosmiEAlgorithm.

Analyze histograms of normalized energy for each ECL crystal from cosmic ray events. Code can either find most-likely energy deposit for each crystal using CRY MC or calibration constant for each crystal (data)

Definition at line 33 of file eclCosmicEAlgorithm.h.

Member Enumeration Documentation

◆ EResult

enum EResult
inherited

The result of calibration.

Enumerator
c_OK 

Finished successfuly =0 in Python.

c_Iterate 

Needs iteration =1 in Python.

c_NotEnoughData 

Needs more data =2 in Python.

c_Failure 

Failed =3 in Python.

c_Undefined 

Not yet known (before execution) =4 in Python.

Definition at line 50 of file CalibrationAlgorithm.h.

Member Function Documentation

◆ calibrate()

CalibrationAlgorithm::EResult calibrate ( )
overrideprotectedvirtual

Run algorithm on events.


..ranges of various fit parameters, and tolerance to determine that fit is at the limit

..Put root into batch mode so that we don't try to open a graphics window


..Clean up existing histograms if necessary


..Histograms containing the data collected by eclCosmicECollectorModule


..Record the number of entries per crystal in each of the two normalized energy histograms and average the constants obtained from DB


..Write out the basic histograms in all cases


..If we have not been asked to do fits, we can quit now


..Check that every crystal has enough entries. If we are finding calibration constants (normal data mode), at least 1 histogram must have sufficient statistics. If we are finding expected values (used with MC), both must have sufficient statistics.


Insufficient data. Quit if we are required to have a successful fit for every crystal


..Some prep for the many fits about to follow

..1D summary histograms

..Histograms to store results for DB


..Loop over specified crystals and performs fits to the two normalized energy distributions

..Extract the 1D normalized energy distribution from the appropriate 2D histogram

..Fit function (xmin, xmax, nparameters) for this histogram

..Estimate initial parameters. For peak, use maximum bin in the allowed range

..Fit range is histogram low edge plus a few bins to peak + 2.5*effective sigma

..Constant from lower edge of plot

..Eta is nominal

..parameters to control iterations. dIter checks if we are stuck in a loop


..Iterate from this point if needed

..Set the initial parameters

..Perform the fit and note the resulting parameters

..The upper fit range should correspond to 20-25% of the peak. Iterate if necessary.

..Check if we are oscillating between two end points

..Many iterations may mean we are stuck in a loop. Try a different end point.

..Set the constant term to 0 if we are close to the limit

..No more than specified number of iterations


..Calculate fit probability. Same as P option in fit, which cannot be used with L


..Fit status

No peak; normalization of Novo component is too small

..poor fit, or relatively poor fit with too many iterations

..parameter at limit

..Store the fit results

..Write out the fit distribution


..Find expected energies from MC, if requested

..Write out expected energies if status is adequate. Check that every crystal has at least one good fit


..Otherwise, find calibration constants

..Find calibration constant separately for the two normalized energy distributions for each crystal

..Peak and uncertainty; assume uncertainties on expected energy and elec calib are negligible

..Find the weighted average of the two constants and store in the histogram

..If both fits failed, use the negative of the initial "same" calibration constant


..Write output to DB if requested and successful

..Store expected energy for each crystal and neighbour type from CRY MC

..Store calibration constant for each crystal (nominally real data)

..Write out some diagnostic histograms

..Histograms containing values written to DB


..Clean up histograms in case Algorithm is called again


..Set the return code appropriately

Implements CalibrationAlgorithm.

Definition at line 49 of file eclCosmicEAlgorithm.cc.

50 {
51 
54  double limitTol = 0.0005; /*< tolerance for checking if a parameter is at the limit */
55  double minFitLimit = 1e-25; /*< cut off for labeling a fit as poor */
56  double minFitProbIter = 1e-8; /*< cut off for labeling a fit as poor if it also has many iterations */
57  double constRatio = 0.5; /*< Novosibirsk normalization must be greater than constRatio x constant term */
58  double peakMin(0.5), peakMax(1.75); /*< range for peak of measured energy distribution */
59  double peakTol = limitTol * (peakMax - peakMin); /*< fit is at limit if it is within peakTol of min or max */
60  double effSigMin(0.08), effSigMax(0.5); /*< range for effective sigma of measured energy distribution */
61  double effSigTol = limitTol * (effSigMax - effSigMin);
62  double etaMin(-3.), etaMax(1.); /*< Novosibirsk tail parameter range */
63  double etaNom(-0.41); /*< Nominal tail parameter */
64  double etaTol = limitTol * (etaMax - etaMin);
65  double constTol = 0.1; /*< if constant is less than constTol, it will be fixed to 0 */
66 
68  gROOT->SetBatch();
69 
72  TH1F* dummy;
73  dummy = (TH1F*)gROOT->FindObject("EnvsCrysSameRing");
74  if (dummy) {delete dummy;}
75  dummy = (TH1F*)gROOT->FindObject("EnvsCrysDifferentRing");
76  if (dummy) {delete dummy;}
77  dummy = (TH1F*)gROOT->FindObject("IntegralVsCrysSame");
78  if (dummy) {delete dummy;}
79  dummy = (TH1F*)gROOT->FindObject("IntegralVsCrysDifferent");
80  if (dummy) {delete dummy;}
81  dummy = (TH1F*)gROOT->FindObject("AverageExpECrysSame");
82  if (dummy) {delete dummy;}
83  dummy = (TH1F*)gROOT->FindObject("AverageExpECrysDifferent");
84  if (dummy) {delete dummy;}
85  dummy = (TH1F*)gROOT->FindObject("AverageElecCalibSame");
86  if (dummy) {delete dummy;}
87  dummy = (TH1F*)gROOT->FindObject("AverageElecCalibDifferent");
88  if (dummy) {delete dummy;}
89  dummy = (TH1F*)gROOT->FindObject("AverageInitialCalibSame");
90  if (dummy) {delete dummy;}
91  dummy = (TH1F*)gROOT->FindObject("AverageInitialCalibDifferent");
92  if (dummy) {delete dummy;}
93 
96  std::vector<std::shared_ptr<TH2F>> EnvsCrys;
97  EnvsCrys.push_back(getObjectPtr<TH2F>("EnvsCrysSameRing"));
98  EnvsCrys.push_back(getObjectPtr<TH2F>("EnvsCrysDifferentRing"));
99 
100  std::vector<std::shared_ptr<TH1F>> ExpEvsCrys;
101  ExpEvsCrys.push_back(getObjectPtr<TH1F>("ExpEvsCrysSameRing"));
102  ExpEvsCrys.push_back(getObjectPtr<TH1F>("ExpEvsCrysDifferentRing"));
103 
104  std::vector<std::shared_ptr<TH1F>> ElecCalibvsCrys;
105  ElecCalibvsCrys.push_back(getObjectPtr<TH1F>("ElecCalibvsCrysSameRing"));
106  ElecCalibvsCrys.push_back(getObjectPtr<TH1F>("ElecCalibvsCrysDifferentRing"));
107 
108  std::vector<std::shared_ptr<TH1F>> InitialCalibvsCrys;
109  InitialCalibvsCrys.push_back(getObjectPtr<TH1F>("InitialCalibvsCrysSameRing"));
110  InitialCalibvsCrys.push_back(getObjectPtr<TH1F>("InitialCalibvsCrysDifferentRing"));
111 
112  std::vector<std::shared_ptr<TH1F>> CalibEntriesvsCrys;
113  CalibEntriesvsCrys.push_back(getObjectPtr<TH1F>("CalibEntriesvsCrysSameRing"));
114  CalibEntriesvsCrys.push_back(getObjectPtr<TH1F>("CalibEntriesvsCrysDifferentRing"));
115 
116  auto RawDigitAmpvsCrys = getObjectPtr<TH2F>("RawDigitAmpvsCrys");
117 
121  TH1F* IntegralVsCrys[2];
122  IntegralVsCrys[0] = new TH1F("IntegralVsCrysSame", "Integral of EnVsCrys for each crystal, same theta ring;Crystal ID", 8736, 0,
123  8736);
124  IntegralVsCrys[1] = new TH1F("IntegralVsCrysDifferent", "Integral of EnVsCrys for each crystal, different theta rings;Crystal ID",
125  8736, 0, 8736);
126 
127  TH1F* AverageExpECrys[2];
128  AverageExpECrys[0] = new TH1F("AverageExpECrysSame",
129  "Average expected E per crys from collector, same theta ring;Crystal ID;Energy (GeV)", 8736, 0, 8736);
130  AverageExpECrys[1] = new TH1F("AverageExpECrysDifferent",
131  "Average expected E per crys from collector, different theta ring;Crystal ID;Energy (GeV)", 8736, 0, 8736);
132 
133  TH1F* AverageElecCalib[2];
134  AverageElecCalib[0] = new TH1F("AverageElecCalibSame",
135  "Average electronics calib const vs crys, same theta ring;Crystal ID;Calibration constant", 8736, 0, 8736);
136  AverageElecCalib[1] = new TH1F("AverageElecCalibDifferent",
137  "Average electronics calib const vs crys, different theta rings;Crystal ID;Calibration constant", 8736, 0, 8736);
138 
139  TH1F* AverageInitialCalib[2];
140  AverageInitialCalib[0] = new TH1F("AverageInitialCalibSame",
141  "Average initial cosmic calib const vs crys, same theta ring;Crystal ID;Calibration constant", 8736, 0, 8736);
142  AverageInitialCalib[1] = new TH1F("AverageInitialCalibDifferent",
143  "Average initial cosmic calib const vs crys, different theta rings;Crystal ID;Calibration constant", 8736, 0, 8736);
144 
145  for (int crysID = 0; crysID < 8736; crysID++) {
146  int histbin = crysID + 1;
147  for (int idir = 0; idir < 2; idir++) {
148  TH1D* hEnergy = EnvsCrys[idir]->ProjectionY("hEnergy", histbin, histbin);
149  int Integral = hEnergy->Integral();
150  IntegralVsCrys[idir]->SetBinContent(histbin, Integral);
151 
152  double TotEntries = CalibEntriesvsCrys[idir]->GetBinContent(histbin);
153 
154  double expectedE = 0.;
155  if (TotEntries > 0.) {expectedE = ExpEvsCrys[idir]->GetBinContent(histbin) / TotEntries;}
156  AverageExpECrys[idir]->SetBinContent(histbin, expectedE);
157  AverageExpECrys[idir]->SetBinError(histbin, 0.);
158 
159  double calibconst = 0.;
160  if (TotEntries > 0.) {calibconst = ElecCalibvsCrys[idir]->GetBinContent(histbin) / TotEntries;}
161  AverageElecCalib[idir]->SetBinContent(histbin, calibconst);
162  AverageElecCalib[idir]->SetBinError(histbin, 0);
163 
164  calibconst = 0.;
165  if (TotEntries > 0.) {calibconst = InitialCalibvsCrys[idir]->GetBinContent(histbin) / TotEntries;}
166  AverageInitialCalib[idir]->SetBinContent(histbin, calibconst);
167  AverageInitialCalib[idir]->SetBinError(histbin, 0);
168  }
169  }
170 
173  TFile* histfile = new TFile("eclCosmicEAlgorithm.root", "recreate");
174  for (int idir = 0; idir < 2; idir++) {
175  EnvsCrys[idir]->Write();
176  IntegralVsCrys[idir]->Write();
177  AverageExpECrys[idir]->Write();
178  AverageElecCalib[idir]->Write();
179  AverageInitialCalib[idir]->Write();
180  }
181  RawDigitAmpvsCrys->Write();
182 
185  if (!performFits) {
186  B2INFO("eclCosmicEAlgorithm has not been asked to perform fits; copying input histograms and quitting");
187  histfile->Close();
188  return c_NotEnoughData;
189  }
190 
195  bool sufficientData = true;
196  for (int crysID = cellIDLo - 1; crysID < cellIDHi; crysID++) {
197  int histbin = crysID + 1;
198  bool SameLow = IntegralVsCrys[0]->GetBinContent(histbin) < minEntries;
199  bool DifferentLow = IntegralVsCrys[1]->GetBinContent(histbin) < minEntries;
200  if ((SameLow && DifferentLow) || (findExpValues && (SameLow || DifferentLow))) {
201  if (storeConst == 1) {B2INFO("eclCosmicEAlgorithm: cellID " << histbin << " has insufficient statistics: " << IntegralVsCrys[0]->GetBinContent(histbin) << " and " << IntegralVsCrys[1]->GetBinContent(histbin) << ". Requirement is " << minEntries);}
202  sufficientData = false;
203  break;
204  }
205  }
206 
209  if (!sufficientData && storeConst == 1) {
210  histfile->Close();
211  return c_NotEnoughData;
212  }
213 
216  TString preName[2] = {"SameRing", "DifferentRing"};
217 
218  TH1F* PeakperCrys[2];
219  PeakperCrys[0] = new TH1F("PeakperCrysSame", "Fit peak per crystal, same theta ring;Crystal ID;Peak normalized energy", 8736, 0,
220  8736);
221  PeakperCrys[1] = new TH1F("PeakperCrysDifferent", "Fit peak per crystal, different theta ring;Crystal ID;Peak normalized energy",
222  8736, 0, 8736);
223 
224  TH1F* SigmaperCrys[2];
225  SigmaperCrys[0] = new TH1F("SigmaperCrysSame", "Fit sigma per crysal, same theta ring;Crystal ID;Sigma (ADC)", 8736, 0, 8736);
226  SigmaperCrys[1] = new TH1F("SigmaperCrysDifferent", "Fit sigma per crysal, different theta ring;Crystal ID;Sigma (ADC)", 8736, 0,
227  8736);
228 
229  TH1F* EtaperCrys[2];
230  EtaperCrys[0] = new TH1F("EtaperCrysSame", "Fit eta per crysal, same theta ring;Crystal ID;Eta", 8736, 0, 8736);
231  EtaperCrys[1] = new TH1F("EtaperCrysDifferent", "Fit eta per crysal, different theta ring;Crystal ID;Eta", 8736, 0, 8736);
232 
233  TH1F* ConstperCrys[2];
234  ConstperCrys[0] = new TH1F("ConstperCrysSame", "Fit constant per crystal, same theta ring;Crystal ID;Constant", 8736, 0, 8736);
235  ConstperCrys[1] = new TH1F("ConstperCrysDifferent", "Fit constant per crystal, different theta ring;Crystal ID;Constant", 8736, 0,
236  8736);
237 
238  TH1F* StatusperCrys[2];
239  StatusperCrys[0] = new TH1F("StatusperCrysSame", "Fit status, same theta ring;Crystal ID;Status", 8736, 0, 8736);
240  StatusperCrys[1] = new TH1F("StatusperCrysDifferent", "Fit status, different theta ring;Crystal ID;Status", 8736, 0, 8736);
241 
243  TH1F* hStatus[2];
244  hStatus[0] = new TH1F("StatusSame", "Fit status, same theta ring", 25, -5, 20);
245  hStatus[1] = new TH1F("StatusDifferent", "Fit status, different theta ring", 25, -5, 20);
246 
247  TH1F* fracPeakUnc[2];
248  fracPeakUnc[0] = new TH1F("fracPeakUncSame", "Fractional uncertainty on peak location, same theta ring", 100, 0, 0.1);
249  fracPeakUnc[1] = new TH1F("fracPeakUncDifferent", "Fractional uncertainty on peak location, different theta ring", 100, 0, 0.1);
250 
251  TH1F* hfitProb[2];
252  hfitProb[0] = new TH1F("fitProbSame", "Probability of fit, same theta ring", 200, 0, 0.02);
253  hfitProb[1] = new TH1F("fitProbDifferent", "Probability of fit, different theta ring", 200, 0, 0.02);
254 
256  TH1F* ExpEnergyperCrys[2];
257  ExpEnergyperCrys[0] = new TH1F("ExpEnergyperCrysSame", "Expected energy per crystal, same theta ring;Crystal ID;Peak energy (GeV)",
258  8736, 0, 8736);
259  ExpEnergyperCrys[1] = new TH1F("ExpEnergyperCrysDifferent",
260  "Expected energy per crystal, different theta ring;Crystal ID;Peak energy (GeV)", 8736, 0, 8736);
261 
262  TH1F* CalibvsCrys = new TH1F("CalibvsCrys", "Energy calibration constant from cosmics;Crystal ID;Calibration constant", 8736, 0,
263  8736);
264 
267  bool allFitsOK = true;
268  for (int crysID = cellIDLo - 1; crysID < cellIDHi; crysID++) {
269  int histbin = crysID + 1;
270  for (int idir = 0; idir < 2; idir++) {
271 
273  TString hname = preName[idir];
274  hname += "Enormalized";
275  hname += crysID;
276  TH1D* hEnergy = EnvsCrys[idir]->ProjectionY(hname, histbin, histbin);
277 
279  double histMin = hEnergy->GetXaxis()->GetXmin();
280  double histMax = hEnergy->GetXaxis()->GetXmax();
281  TF1* func = new TF1("eclCosmicNovoConst", eclCosmicNovoConst, histMin, histMax, 5);
282  func->SetParNames("normalization", "peak", "effSigma", "eta", "const");
283  func->SetParLimits(1, peakMin, peakMax);
284  func->SetParLimits(2, effSigMin, effSigMax);
285  func->SetParLimits(3, etaMin, etaMax);
286 
288  hEnergy->GetXaxis()->SetRangeUser(peakMin, peakMax);
289  int maxBin = hEnergy->GetMaximumBin();
290  double peakE = hEnergy->GetBinLowEdge(maxBin);
291  double peakEUnc = 0.;
292  double normalization = hEnergy->GetMaximum();
293  double effSigma = hEnergy->GetRMS();
294  double sigmaUnc = 0.;
295  hEnergy->GetXaxis()->SetRangeUser(histMin, histMax);
296 
298  double fitlow = 0.25;
299  double fithigh = peakE + 2.5 * effSigma;
300 
302  int il0 = hEnergy->GetXaxis()->FindBin(fitlow);
303  int il1 = hEnergy->GetXaxis()->FindBin(fitlow + 0.1);
304  double constant = hEnergy->Integral(il0, il1) / (1 + il1 - il0);
305  double constUnc = 0.;
306 
308  double eta = etaNom;
309  double etaUnc = 0.;
310 
312  double dIter = 0.1 * (histMax - histMin) / hEnergy->GetNbinsX();
313  double highold(0.), higholdold(0.);
314  double fitProb(0.);
315  double fitProbDefault(0.);
316  bool fitHist = IntegralVsCrys[idir]->GetBinContent(histbin) >= minEntries; /* fit only if enough events */
317  bool fixConst = false;
318  int nIter = 0;
319 
322  while (fitHist) {
323  nIter++;
324 
326  func->SetParameters(normalization, peakE, effSigma, eta, constant);
327  if (fixConst) { func->FixParameter(4, 0); }
328 
330  hEnergy->Fit(func, "LIQ", "", fitlow, fithigh);
331  normalization = func->GetParameter(0);
332  peakE = func->GetParameter(1);
333  peakEUnc = func->GetParError(1);
334  effSigma = func->GetParameter(2);
335  sigmaUnc = func->GetParError(2);
336  eta = func->GetParameter(3);
337  etaUnc = func->GetParError(3);
338  constant = func->GetParameter(4);
339  constUnc = func->GetParError(4);
340  fitProbDefault = func->GetProb();
341 
343  fitHist = false;
344  double peak = func->Eval(peakE) - constant;
345  double tRatio = (func->Eval(fithigh) - constant) / peak;
346  if (tRatio < tRatioMin || tRatio > tRatioMax) {
347  double targetY = constant + 0.5 * (tRatioMin + tRatioMax) * peak;
348  higholdold = highold;
349  highold = fithigh;
350  fithigh = func->GetX(targetY, peakE, histMax);
351  fitHist = true;
352 
354  if (abs(fithigh - higholdold) < dIter) {fithigh = 0.5 * (highold + higholdold); }
355 
357  if (nIter > maxIterations - 3) {fithigh = 0.33333 * (fithigh + highold + higholdold); }
358  }
359 
361  if (constant < constTol && !fixConst) {
362  constant = 0;
363  fixConst = true;
364  fitHist = true;
365  }
366 
368  if (nIter == maxIterations) {fitHist = false;}
369  B2DEBUG(200, "cellID = " << histbin << " " << nIter << " " << preName[idir] << " " << peakE << " " << constant << " " << tRatio <<
370  " " << fithigh);
371 
372  }
373 
376  fitProb = 0.;
377  if (nIter > 0) {
378  int lowbin = hEnergy->GetXaxis()->FindBin(fitlow);
379  int highbin = hEnergy->GetXaxis()->FindBin(fithigh);
380  int npar = 5;
381  if (fixConst) {npar = 4;}
382  int ndeg = (highbin - lowbin) + 1 - npar;
383  double chisq = 0.;
384  double halfbinwidth = 0.5 * hEnergy->GetBinWidth(1);
385  for (int ib = lowbin; ib <= highbin; ib++) {
386  double yexp = func->Eval(hEnergy->GetBinLowEdge(ib) + halfbinwidth);
387  double yobs = hEnergy->GetBinContent(ib);
388  double dchi2 = (yexp - yobs) * (yexp - yobs) / yexp;
389  chisq += dchi2;
390  }
391  fitProb = 0.5 * (TMath::Prob(chisq, ndeg) + fitProbDefault);
392  }
393 
396  int iStatus = fitOK; // success
397  if (nIter == maxIterations) {iStatus = iterations;} // too many iterations
398 
400  if (normalization < constRatio * constant) {iStatus = noPeak;}
401 
403  if (fitProb <= minFitLimit || (fitProb < minFitProbIter && iStatus == iterations)) {iStatus = poorFit;}
404 
406  if ((peakE < peakMin + peakTol) || (peakE > peakMax - peakTol)) {iStatus = atLimit;}
407  if ((effSigma < effSigMin + effSigTol) || (effSigma > effSigMax - effSigTol)) {iStatus = atLimit;}
408  if ((eta < etaMin + etaTol) || (eta > etaMax - etaTol)) {iStatus = atLimit;}
409 
410  //**..No fit
411  if (nIter == 0) {iStatus = notFit;} // not fit
412 
414  PeakperCrys[idir]->SetBinContent(histbin, peakE);
415  PeakperCrys[idir]->SetBinError(histbin, peakEUnc);
416  SigmaperCrys[idir]->SetBinContent(histbin, effSigma);
417  SigmaperCrys[idir]->SetBinError(histbin, sigmaUnc);
418  EtaperCrys[idir]->SetBinContent(histbin, eta);
419  EtaperCrys[idir]->SetBinError(histbin, etaUnc);
420  ConstperCrys[idir]->SetBinContent(histbin, constant);
421  ConstperCrys[idir]->SetBinError(histbin, constUnc);
422  StatusperCrys[idir]->SetBinContent(histbin, iStatus);
423  hStatus[idir]->Fill(iStatus);
424  fracPeakUnc[idir]->Fill(peakEUnc / peakE);
425  hfitProb[idir]->Fill(fitProb);
426 
428  B2INFO("cellID " << histbin << " " << preName[idir] << " status = " << iStatus << " fit probability = " << fitProb);
429  histfile->cd();
430  hEnergy->Write();
431  }
432  }
433 
436  if (findExpValues) {
437 
439  for (int crysID = 0; crysID < 8736; crysID++) {
440  int histbin = crysID + 1;
441  bool atLeastOneOK = false;
442  for (int idir = 0; idir < 2; idir++) {
443  double fitstatus = StatusperCrys[idir]->GetBinContent(histbin);
444  double peakE = PeakperCrys[idir]->GetBinContent(histbin);
445  double peakEUnc = PeakperCrys[idir]->GetBinError(histbin);
446 
447  //**..For failed fits, store the negative of the input expected energy */
448  if (fitstatus < iterations) {
449  if (histbin >= cellIDLo && histbin <= cellIDHi) {
450  B2INFO("eclCosmicEAlgorithm: crystal " << crysID << " " << preName[idir] << " is not a successful fit. Status = " << fitstatus);
451  }
452  peakE = -1.;
453  peakEUnc = 0.;
454  } else {
455  atLeastOneOK = true;
456  }
457  double inputExpE = abs(AverageExpECrys[idir]->GetBinContent(histbin));
458  ExpEnergyperCrys[idir]->SetBinContent(histbin, inputExpE * peakE);
459  ExpEnergyperCrys[idir]->SetBinError(histbin, inputExpE * peakEUnc / peakE);
460  }
461  if (!atLeastOneOK) {allFitsOK = false;}
462  }
463 
466  } else {
467 
469  for (int crysID = 0; crysID < 8736; crysID++) {
470  int histbin = crysID + 1;
471  double calibConst[2] = {};
472  double calibConstUnc[2] = {999999., 999999.};
473  double weight[2] = {};
474  bool bothFitsBad = true;
475  for (int idir = 0; idir < 2; idir++) {
476 
478  double peakE = PeakperCrys[idir]->GetBinContent(histbin);
479  double fracPeakEUnc = PeakperCrys[idir]->GetBinError(histbin) / peakE;
480  double inputConst = AverageInitialCalib[idir]->GetBinContent(histbin);
481  double fitstatus = StatusperCrys[idir]->GetBinContent(histbin);
482  double inputExpE = AverageExpECrys[idir]->GetBinContent(histbin);
483  if (fitstatus >= iterations && inputConst == 0) {B2FATAL("eclCosmicEAlgorithm: input calibration = 0 for idir = " << idir << " and crysID = " << crysID);}
484 
485  //** Find constant only if fit was successful and we have a value for the expected energy */
486  if (fitstatus >= iterations && inputExpE > 0.) {
487  calibConst[idir] = abs(inputConst) / peakE;
488  calibConstUnc[idir] = calibConst[idir] * fracPeakEUnc / peakE;
489  weight[idir] = 1. / (calibConstUnc[idir] * calibConstUnc[idir]);
490  bothFitsBad = false;
491  }
492  if (fitstatus < iterations && histbin >= cellIDLo && histbin <= cellIDHi) {
493  B2INFO("eclCosmicEAlgorithm: cellID " << histbin << " " << preName[idir] << " is not a successful fit. Status = " << fitstatus);
494  } else if (inputExpE < 0. && histbin >= cellIDLo && histbin <= cellIDHi) {
495  B2WARNING("eclCosmicEAlgorithm: cellID " << histbin << " " << preName[idir] << " has no expected energy. Status = " << fitstatus);
496  }
497  }
498 
499 
501  double averageConst;
502  double averageConstUnc;
503 
505  if (bothFitsBad) {
506  if (histbin >= cellIDLo && histbin <= cellIDHi) {B2INFO("eclCosmicEAlgorithm: no constant found for cellID = " << histbin << " status = " << StatusperCrys[0]->GetBinContent(histbin) << " and " << StatusperCrys[1]->GetBinContent(histbin));}
507  averageConst = -1.*abs(AverageInitialCalib[0]->GetBinContent(histbin));
508  averageConstUnc = 0.;
509  } else {
510  averageConst = (calibConst[0] * weight[0] + calibConst[1] * weight[1]) / (weight[0] + weight[1]);
511  averageConstUnc = 1. / sqrt(weight[0] + weight[1]);
512  }
513  CalibvsCrys->SetBinContent(histbin, averageConst);
514  CalibvsCrys->SetBinError(histbin, averageConstUnc);
515  }
516  }
517 
520  bool DBsuccess = false;
521  if (storeConst == 0 || (storeConst == 1 && allFitsOK)) {
522  DBsuccess = true;
523 
525  if (findExpValues) {
526  std::vector<std::string> DBname = {"ECLExpCosmicESame", "ECLExpCosmicEDifferent"};
527  for (int idir = 0; idir < 2; idir++) {
528  std::vector<float> tempE;
529  std::vector<float> tempUnc;
530  for (int crysID = 0; crysID < 8736; crysID++) {
531  int histbin = crysID + 1;
532  tempE.push_back(ExpEnergyperCrys[idir]->GetBinContent(histbin));
533  tempUnc.push_back(ExpEnergyperCrys[idir]->GetBinError(histbin));
534  }
535  ECLCrystalCalib* ExpectedE = new ECLCrystalCalib();
536  ExpectedE->setCalibVector(tempE, tempUnc);
537  saveCalibration(ExpectedE, DBname[idir]);
538  B2INFO("eclCosmicEAlgorithm: successfully stored expected values for " << DBname[idir]);
539  }
540 
542  } else {
543  std::vector<float> tempCalib;
544  std::vector<float> tempCalibUnc;
545  for (int crysID = 0; crysID < 8736; crysID++) {
546  int histbin = crysID + 1;
547  tempCalib.push_back(CalibvsCrys->GetBinContent(histbin));
548  tempCalibUnc.push_back(CalibvsCrys->GetBinError(histbin));
549  }
550  ECLCrystalCalib* CosmicECalib = new ECLCrystalCalib();
551  CosmicECalib->setCalibVector(tempCalib, tempCalibUnc);
552  saveCalibration(CosmicECalib, "ECLCrystalEnergyCosmic");
553  B2INFO("eclCosmicEAlgorithm: successfully stored calibration constants");
554  }
555  }
556 
558  for (int idir = 0; idir < 2; idir++) {
559  PeakperCrys[idir]->Write();
560  SigmaperCrys[idir]->Write();
561  EtaperCrys[idir]->Write();
562  ConstperCrys[idir]->Write();
563  StatusperCrys[idir]->Write();
564  hStatus[idir]->Write();
565  fracPeakUnc[idir]->Write();
566  hfitProb[idir]->Write();
567  }
568 
570  if (findExpValues) {
571  ExpEnergyperCrys[0]->Write();
572  ExpEnergyperCrys[1]->Write();
573  } else {
574  CalibvsCrys->Write();
575  }
576  histfile->Close();
577 
580  dummy = (TH1F*)gROOT->FindObject("PeakperCrysSame"); delete dummy;
581  dummy = (TH1F*)gROOT->FindObject("SigmaperCrysSame"); delete dummy;
582  dummy = (TH1F*)gROOT->FindObject("EtaperCrysSame"); delete dummy;
583  dummy = (TH1F*)gROOT->FindObject("ConstperCrysSame"); delete dummy;
584  dummy = (TH1F*)gROOT->FindObject("StatusperCrysSame"); delete dummy;
585  dummy = (TH1F*)gROOT->FindObject("StatusSame"); delete dummy;
586  dummy = (TH1F*)gROOT->FindObject("fracPeakUncSame"); delete dummy;
587  dummy = (TH1F*)gROOT->FindObject("fitProbSame"); delete dummy;
588  dummy = (TH1F*)gROOT->FindObject("ExpEnergyperCrysSame"); delete dummy;
589  dummy = (TH1F*)gROOT->FindObject("PeakperCrysDifferent"); delete dummy;
590  dummy = (TH1F*)gROOT->FindObject("SigmaperCrysDifferent"); delete dummy;
591  dummy = (TH1F*)gROOT->FindObject("EtaperCrysDifferent"); delete dummy;
592  dummy = (TH1F*)gROOT->FindObject("ConstperCrysDifferent"); delete dummy;
593  dummy = (TH1F*)gROOT->FindObject("StatusperCrysDifferent"); delete dummy;
594  dummy = (TH1F*)gROOT->FindObject("StatusDifferent"); delete dummy;
595  dummy = (TH1F*)gROOT->FindObject("fracPeakUncDifferent"); delete dummy;
596  dummy = (TH1F*)gROOT->FindObject("fitProbDifferent"); delete dummy;
597  dummy = (TH1F*)gROOT->FindObject("ExpEnergyperCrysDifferent"); delete dummy;
598  dummy = (TH1F*)gROOT->FindObject("CalibvsCrys"); delete dummy;
599 
602  if (storeConst == -1) {
603  B2INFO("eclCosmicEAlgorithm performed fits but was not asked to store contants");
604  return c_Failure;
605  } else if (!DBsuccess) {
606  if (findExpValues) { B2INFO("eclCosmicEAlgorithm: failed to store expected values"); }
607  else { B2INFO("eclCosmicEAlgorithm: failed to store calibration constants"); }
608  return c_Failure;
609  }
610  return c_OK;
611 }

◆ checkPyExpRun()

bool checkPyExpRun ( PyObject *  pyObj)
inherited

Checks that a PyObject can be successfully converted to an ExpRun type.

Checks if the PyObject can be converted to ExpRun.

Definition at line 21 of file CalibrationAlgorithm.cc.

◆ convertPyExpRun()

ExpRun convertPyExpRun ( PyObject *  pyObj)
inherited

Performs the conversion of PyObject to ExpRun.

Converts the PyObject to an ExpRun. We've preoviously checked the object so this assumes a lot about the PyObject.

Definition at line 63 of file CalibrationAlgorithm.cc.

◆ execute()

EResult execute ( std::vector< Calibration::ExpRun >  runs = {},
int  iteration = 0,
IntervalOfValidity  iov = IntervalOfValidity() 
)
inherited

Runs calibration over vector of runs for a given iteration.

You can also specify the IoV to save the database payload as. By default the Algorithm will create an IoV from your requested ExpRuns, or from the overall ExpRuns of the input data if you haven't specified ExpRuns in this function.

No checks are performed to make sure that a IoV you specify matches the data you ran over, it simply labels the IoV to commit to the database later.

◆ getCollectorName()

std::string getCollectorName ( ) const
inlineinherited

Alias for prefix.

For convenience and less writing, we say developers to set this to default collector module name in constructor of base class. One can however use the dublets of collector+algorithm multiple times with different settings. To bind these together correctly, the prefix has to be set the same for algo and collector. So we call the setter setPrefix rather than setModuleName or whatever. This getter will work out of the box for default cases -> return the name of module you have to add to your path to collect data for this algorihtm.

Definition at line 174 of file CalibrationAlgorithm.h.

◆ setInputFileNames() [1/2]

void setInputFileNames ( PyObject *  inputFileNames)
inherited

Set the input file names used for this algorithm from a Python list.

Set the input file names used for this algorithm and resolve the wildcards.

Definition at line 159 of file CalibrationAlgorithm.cc.

◆ setInputFileNames() [2/2]

void setInputFileNames ( std::vector< std::string >  inputFileNames)
protectedinherited

Set the input file names used for this algorithm.

Set the input file names used for this algorithm and resolve the wildcards.

Definition at line 187 of file CalibrationAlgorithm.cc.

Member Data Documentation

◆ cellIDLo

int cellIDLo

Parameters to control Novosibirsk fit to signal measured in each crystal.

First cellID to be fit

Definition at line 43 of file eclCosmicEAlgorithm.h.

◆ storeConst

int storeConst

controls which values are written to the database.

0 (default): store value found by successful fits, or -|input value| otherwise; -1 : do not store values 1 : store values if every fit for [cellIDLo,cellIDHi] was successful

Definition at line 51 of file eclCosmicEAlgorithm.h.


The documentation for this class was generated from the following files:
Belle2::ECL::eclCosmicEAlgorithm::cellIDLo
int cellIDLo
Parameters to control Novosibirsk fit to signal measured in each crystal.
Definition: eclCosmicEAlgorithm.h:43
Belle2::ECL::eclCosmicEAlgorithm::iterations
int iterations
fit reached max number of iterations, but is useable
Definition: eclCosmicEAlgorithm.h:63
prepareAsicCrosstalkSimDB.e
e
aux.
Definition: prepareAsicCrosstalkSimDB.py:53
Belle2::ECL::eclCosmicEAlgorithm::atLimit
int atLimit
a parameter is at the limit; fit not useable
Definition: eclCosmicEAlgorithm.h:64
Belle2::CalibrationAlgorithm::saveCalibration
void saveCalibration(TClonesArray *data, const std::string &name)
Store DBArray payload with given name with default IOV.
Definition: CalibrationAlgorithm.cc:290
Belle2::ECL::eclCosmicEAlgorithm::poorFit
int poorFit
low chi square; fit not useable
Definition: eclCosmicEAlgorithm.h:65
Belle2::CalibrationAlgorithm::c_OK
@ c_OK
Finished successfuly =0 in Python.
Definition: CalibrationAlgorithm.h:51
Belle2::ECLCrystalCalib
General DB object to store one calibration number per ECL crystal.
Definition: ECLCrystalCalib.h:34
Belle2::ECL::eclCosmicEAlgorithm::noPeak
int noPeak
Novosibirsk component of fit is negligible; fit not useable.
Definition: eclCosmicEAlgorithm.h:66
Belle2::ECL::eclCosmicEAlgorithm::cellIDHi
int cellIDHi
Last cellID to be fit.
Definition: eclCosmicEAlgorithm.h:44
Belle2::ECL::eclCosmicEAlgorithm::storeConst
int storeConst
controls which values are written to the database.
Definition: eclCosmicEAlgorithm.h:51
Belle2::ECL::eclCosmicEAlgorithm::tRatioMax
double tRatioMax
Fit range is adjusted so that fit at upper endpoint is between tRatioMin and tRatioMax of peak.
Definition: eclCosmicEAlgorithm.h:48
Belle2::CalibrationAlgorithm::c_Failure
@ c_Failure
Failed =3 in Python.
Definition: CalibrationAlgorithm.h:54
Belle2::ECL::eclCosmicEAlgorithm::fitOK
int fitOK
fit is OK
Definition: eclCosmicEAlgorithm.h:62
Belle2::ECL::eclCosmicEAlgorithm::maxIterations
int maxIterations
no more than maxIteration iterations
Definition: eclCosmicEAlgorithm.h:46
Belle2::ECL::eclCosmicEAlgorithm::minEntries
int minEntries
All crystals to be fit must have at least minEntries events in the fit range.
Definition: eclCosmicEAlgorithm.h:45
Belle2::CalibrationAlgorithm::c_NotEnoughData
@ c_NotEnoughData
Needs more data =2 in Python.
Definition: CalibrationAlgorithm.h:53
Belle2::ECL::eclCosmicEAlgorithm::performFits
bool performFits
if false, input histograms are copied to output, but no fits are done.
Definition: eclCosmicEAlgorithm.h:49
Belle2::ECLCrystalCalib::setCalibVector
void setCalibVector(const std::vector< float > &CalibConst, const std::vector< float > &CalibConstUnc)
Set vector of constants with uncertainties.
Definition: ECLCrystalCalib.h:48
Belle2::ECL::eclCosmicEAlgorithm::findExpValues
bool findExpValues
if true, fits are used to find expected energy deposit for each crystal instead of the calibration co...
Definition: eclCosmicEAlgorithm.h:50
Belle2::ECL::eclCosmicEAlgorithm::tRatioMin
double tRatioMin
Fit range is adjusted so that fit at upper endpoint is between tRatioMin and tRatioMax of peak.
Definition: eclCosmicEAlgorithm.h:47
Belle2::ECL::eclCosmicEAlgorithm::notFit
int notFit
no fit performed
Definition: eclCosmicEAlgorithm.h:67