Belle II Software  release-05-02-19
CDCDedxMomentumAlgorithm.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2016 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: jvbennett *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <reconstruction/calibration/CDCDedxMomentumAlgorithm.h>
12 
13 #include <reconstruction/dbobjects/CDCDedxMomentumCor.h>
14 
15 #include <TF1.h>
16 #include <TH1F.h>
17 #include <TCanvas.h>
18 
19 using namespace Belle2;
20 
21 
22 //-----------------------------------------------------------------
23 // Implementation
24 //-----------------------------------------------------------------
25 
27 {
28  // Set module properties
29  setDescription("A calibration algorithm for CDC dE/dx electron cos(theta) dependence");
30 }
31 
32 //-----------------------------------------------------------------
33 // Run the calibration
34 //-----------------------------------------------------------------
35 
37 {
38  // Get data objects
39  auto ttree = getObjectPtr<TTree>("tree");
40 
41  // require at least 100 tracks (arbitrary for now)
42  if (ttree->GetEntries() < 100)
43  return c_NotEnoughData;
44 
45  double dedx, p;
46  ttree->SetBranchAddress("dedx", &dedx);
47  ttree->SetBranchAddress("p", &p);
48 
49  // make histograms to store dE/dx values in bins of cos(theta)
50  // bin size can be arbitrary, for now just make uniform bins
51  const int nbins = 100;
52  TH1F dedxp[nbins];
53  for (unsigned int i = 0; i < nbins; ++i) {
54  dedxp[i] = TH1F(TString::Format("dedxp%d", i), "dE/dx in bins of momentum", 200, 0, 4);
55  }
56 
57  // fill histograms, bin size may be arbitrary
58  for (int i = 0; i < ttree->GetEntries(); ++i) {
59  ttree->GetEvent(i);
60  if (dedx == 0) continue;
61  if (p <= 0.0 || p > 10.0) continue;
62  int bin = std::floor(p / (10.0 / nbins));
63  if (bin < 0 || bin >= nbins) continue;
64  dedxp[bin].Fill(dedx);
65  }
66 
67  // Print the histograms for quality control
68  TCanvas* ctmp = new TCanvas("tmp", "tmp", 900, 900);
69  ctmp->Divide(3, 3);
70  std::stringstream psname; psname << "dedx_momentum.ps[";
71  ctmp->Print(psname.str().c_str());
72  psname.str(""); psname << "dedx_momentum.ps";
73 
74  // fit histograms to get gains in bins of cos(theta)
75  std::vector<double> momentum;
76  for (unsigned int i = 0; i < nbins; ++i) {
77  ctmp->cd(i % 9 + 1); // each canvas is 9x9
78  dedxp[i].DrawCopy("hist");
79 
80  double mean = 1.0;
81  if (dedxp[i].Integral() < 10)
82  momentum.push_back(mean); // FIXME! --> should return not enough data
83  else {
84  if (dedxp[i].Fit("gaus")) {
85  B2WARNING("Fit failed...");
86  momentum.push_back(mean); // FIXME! --> should return not enough data
87  } else {
88  B2WARNING("Fit succeeded: " << mean << "\t" << dedxp[i].GetFunction("gaus")->GetParameter(1));
89  mean *= dedxp[i].GetFunction("gaus")->GetParameter(1);
90  momentum.push_back(mean);
91  }
92  }
93  if ((i + 1) % 9 == 0)
94  ctmp->Print(psname.str().c_str());
95  }
96 
97  psname.str(""); psname << "dedx_momentum.ps]";
98  ctmp->Print(psname.str().c_str());
99 
100  B2INFO("dE/dx Calibration done for CDC dE/dx momentum correction");
101  std::cout << "dE/dx Calibration done for CDC dE/dx momentum correction" << std::endl;
102 
103  CDCDedxMomentumCor* gain = new CDCDedxMomentumCor(momentum);
104  saveCalibration(gain, "CDCDedxMomentumCor");
105 
106  delete ctmp;
107 
108  return c_OK;
109 }
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::CDCDedxMomentumAlgorithm::calibrate
virtual EResult calibrate() override
Scale Momentum algorithm.
Definition: CDCDedxMomentumAlgorithm.cc:36
Belle2::CalibrationAlgorithm::c_OK
@ c_OK
Finished successfuly =0 in Python.
Definition: CalibrationAlgorithm.h:51
Belle2::CDCDedxMomentumAlgorithm::CDCDedxMomentumAlgorithm
CDCDedxMomentumAlgorithm()
Constructor: Sets the description, the properties and the parameters of the algorithm.
Definition: CDCDedxMomentumAlgorithm.cc:26
Belle2::CalibrationAlgorithm::setDescription
void setDescription(const std::string &description)
Set algorithm description (in constructor)
Definition: CalibrationAlgorithm.h:331
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::CalibrationAlgorithm::EResult
EResult
The result of calibration.
Definition: CalibrationAlgorithm.h:50
Belle2::CalibrationAlgorithm::c_NotEnoughData
@ c_NotEnoughData
Needs more data =2 in Python.
Definition: CalibrationAlgorithm.h:53
Belle2::CDCDedxMomentumCor
dE/dx wire gain calibration constants
Definition: CDCDedxMomentumCor.h:36
Belle2::CalibrationAlgorithm
Base class for calibration algorithms.
Definition: CalibrationAlgorithm.h:47