Belle II Software development
CDCInitialT0Determination.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8#include "cdc/modules/cdcInitialT0Determination/CDCInitialT0Determination.h"
9#include <cdc/geometry/CDCGeometryPar.h>
10#include <framework/gearbox/Const.h>
11#include "TF1.h"
12#include "TDirectory.h"
13#include "TFile.h"
14#include "TGraphErrors.h"
15#include "TROOT.h"
16
17using namespace std;
18using namespace Belle2;
19using namespace CDC;
20REG_MODULE(CDCInitialT0Determination);
21
23{
24 setDescription("Module to determine crude t0");
25 setPropertyFlags(c_ParallelProcessingCertified); // specify this flag if you need parallel processing
26 addParam("OutputFileName", m_outputFileName, "t0 file output file", std::string("t0.dat"));
27 addParam("LowTDC", m_tdcMin, "lower boundary of tdc histogram", m_tdcMin);
28 addParam("UpTDC", m_tdcMax, "Upper boundary of tdc histogram", m_tdcMax);
29 addParam("InitialT0", m_initT0, "initial t0 for fitting", 3579.);
30 addParam("Cosmic", m_cosmic, "true; tof negative for upper part of cdc", true);
31 addParam("Zoffset", m_zOffset, "z offset", 0.);
32 addParam("ADCCut", m_adcMin, "threshold of ADC", m_adcMin);
33 addParam("StoreFittedHisto", m_storeFittedHisto, "Store fitted histogram for each channel or not", true);
34 addParam("HistoFileName", m_histoFileName, "file contain TDC histo", std::string("TDC.root"));
35 addParam("MinEntries", m_minEntries, "minimum entries per channel", m_minEntries);
36}
37
41
42
44{
46 for (int il = 0; il < 56; ++il) {
47 for (unsigned short w = 0; w < cdcgeo.nWiresInLayer(il); ++w) {
48 m_hTDC[il][w] = new TH1D(Form("hLay%d_ch%d", il, w), "tdc", m_tdcMax - m_tdcMin, m_tdcMin, m_tdcMax);
49 }
50 }
51 for (int ib = 0; ib < 300; ++ib) {
52 m_hTDCBoard[ib] = new TH1D(Form("hTDCBoard%d", ib), "", m_tdcMax - m_tdcMin, m_tdcMin, m_tdcMax);
53 }
54 m_hT0All = new TH1D("hT0All", "", 8500, 0, 8500);
55 m_CDCHits.isRequired();
56}
58{
60 for (const auto& hit : m_CDCHits) {
61 WireID wireid(hit.getID());
62 unsigned short lay = wireid.getICLayer();
63 unsigned short w = wireid.getIWire();
64 if (hit.getADCCount() > m_adcMin) {
65 m_hTDC[lay][w]->Fill(hit.getTDCCount());
66 m_hTDCBoard[cdcgeo.getBoardID(WireID(lay, w))]->Fill(hit.getTDCCount());
67 }
68 }
69}
71{
72 gROOT->SetBatch(1);
73 std::vector<double> sb;
74 std::vector<double> dsb;
75 std::vector<double> t0b;
76 std::vector<double> dt0b;
77 std::vector<double> b;
78 std::vector<double> db;
79 TH1D* hs = new TH1D("hs", "sigma", 100, 0, 20);
80
82 TF1* f1 = new TF1("f1", "[0]+[1]*(exp([2]*(x-[3]))/(1+exp(-([4]-x)/[5])))", m_tdcMin, m_tdcMax);
83 f1->SetParLimits(0, 0., 1000.);
84 f1->SetLineColor(kRed);
85 double tdcBinWidth = cdcgeo.getTdcBinWidth();
86 double bflag[300];
87 for (int ib = 1; ib < 300; ++ib) {
88 if (m_hTDCBoard[ib]->GetEntries() < m_minEntries) {
89 B2DEBUG(199, "Warning: this board low statistic: " << m_hTDCBoard[ib]->GetEntries());
90 bflag[ib] = 0;
91 m_t0b[ib] = m_initT0;
92 continue;
93 }
94 double p3 = m_hTDCBoard[ib]->GetXaxis()->GetBinCenter(m_hTDCBoard[ib]->GetMaximumBin());
95 f1->SetParameters(0, m_hTDCBoard[ib]->GetMaximum(), -0.001, p3, m_initT0, 2.5);
96 m_hTDCBoard[ib]->Fit("f1", "QM", "", m_initT0 - 60, m_initT0 + 60);
97
98 if ((fabs(f1->GetParameter(4) - m_initT0) > 100)
99 || (fabs(f1->GetParameter(5)) < 0.01)
100 || (fabs(f1->GetParameter(5)) > 16)) {
101
102 bflag[ib] = 0;
103 m_t0b[ib] = m_initT0;
104 continue;
105 }
106
107 bflag[ib] = 1;
108 m_t0b[ib] = f1->GetParameter(4) * tdcBinWidth;
109
110 sb.push_back(f1->GetParameter(5));
111 dsb.push_back(f1->GetParError(5));
112 t0b.push_back(f1->GetParameter(4));
113 dt0b.push_back(f1->GetParError(4));
114 b.push_back(ib);
115 db.push_back(0);
116 }
117
118 for (int il = 0; il < 56; ++il) {
119 for (unsigned short w = 0; w < cdcgeo.nWiresInLayer(il); ++w) {
120 B2DEBUG(99, "fitting for channel: " << il << " - " << w);
121 B2DEBUG(99, "number of entries" << m_hTDC[il][w]->GetEntries());
122 m_t0[il][w] = m_initT0 * tdcBinWidth;
123 int bid = cdcgeo.getBoardID(WireID(il, w));
124 if (m_hTDC[il][w]->GetEntries() < m_minEntries) {
125 B2DEBUG(99, "Warning: low statistic channel: " << m_hTDC[il][w]->GetEntries());
126 if (bflag[bid] != 0) {
127 m_t0[il][w] = m_t0b[bid];
128 m_flag[il][w] = true;
129 } else {m_flag[il][w] = false;}
130 } else {
131 double p3 = m_hTDC[il][w]->GetXaxis()->GetBinCenter(m_hTDC[il][w]->GetMaximumBin());
132 f1->SetParameters(0, m_hTDC[il][w]->GetMaximum(), -0.001, p3, m_initT0, 2.5);
133 m_hTDC[il][w]->Fit("f1", "QM", "", m_initT0 - 60, m_initT0 + 60);
134 B2DEBUG(99, "prob of fit : " << f1->GetProb());
135 if ((f1->GetProb() < 1E-150) || (fabs(f1->GetParameter(4) - m_initT0) > 100) || (f1->GetParameter(5) < 0.1)
136 || (f1->GetParameter(5) > 20)) {
137 if (bflag[bid] != 0) {
138 m_t0[il][w] = m_t0b[bid];
139 m_flag[il][w] = true;
140 } else {m_flag[il][w] = false;}
141 } else {
142 m_t0[il][w] = f1->GetParameter(4) * tdcBinWidth;
143 hs->Fill(f1->GetParameter(5));
144 m_flag[il][w] = true;
145 }
146 }
147 B2DEBUG(99, "P4 = " << m_t0[il][w]);
148 if (m_cosmic && cdcgeo.wireBackwardPosition(il, w).Y() > 0) {
149 m_t0[il][w] -= cdcgeo.senseWireR(il) / Const::speedOfLight;
150 } else {
151 m_t0[il][w] += cdcgeo.senseWireR(il) / Const::speedOfLight;
152 }
153 m_t0[il][w] += (m_zOffset - cdcgeo.wireBackwardPosition(il, w).Z()) / 27.25;
154 m_t0[il][w] += 6.122;
155 m_hT0All->Fill(m_t0[il][w]);
156 // m_hT0b[cdcgeo.getBoardID(WireID(il, w))]->Fill(m_t0[il][w]);
157 }
158 }
159
160 //check t0, and add t0 for low static channel
161 ofstream ofs(m_outputFileName.c_str());
162 for (int il = 0; il < 56; ++il) {
163 for (unsigned short w = 0; w < cdcgeo.nWiresInLayer(il); ++w) {
164 if (m_flag[il][w] != true) {
165 m_t0[il][w] = m_hT0All->GetMean();
166 }
167 ofs << il << "\t" << w << "\t" << m_t0[il][w] << endl;
168 }
169 }
170 ofs.close();
171 if (m_storeFittedHisto) {
172 TFile* fhist = new TFile(m_histoFileName.c_str(), "recreate");
173 fhist->cd();
174 TDirectory* top = gDirectory;
175 TDirectory* Direct[56];
176 for (int il = 0; il < 56; ++il) {
177 top->cd();
178 Direct[il] = gDirectory->mkdir(Form("lay_%d", il));
179 Direct[il]->cd();
180 for (unsigned short w = 0; w < cdcgeo.nWiresInLayer(il); ++w) {
181 if (m_flag[il][w] == true) {
182 m_hTDC[il][w]->Write();
183 }
184 }
185 }
186 top->cd();
187 TDirectory* board = gDirectory->mkdir("board");
188 board->cd();
189 for (int ib = 0; ib < 300; ++ib) {
190 if (m_hTDCBoard[ib]) {
191 m_hTDCBoard[ib]->Write();
192 }
193 }
194 top->cd();
195 m_hT0All->Write();
196 hs->Write();
197 if (b.size() > 20) {
198 TGraphErrors* gr = new TGraphErrors(b.size(), &b.at(0), &sb.at(0), &db.at(0), &dsb.at(0));
199 gr->SetName("reso");
200 gr->Write();
201 TGraphErrors* grT0b = new TGraphErrors(b.size(), &b.at(0), &t0b.at(0), &db.at(0), &dt0b.at(0));
202 grT0b->SetName("T0Board");
203 grT0b->Write();
204 }
205 fhist->Close();
206 }
207}
208
R E
internal precision of FFTW codelets
DataType Y() const
access variable Y (= .at(1) without boundary check)
Definition B2Vector3.h:433
TH1D * m_hT0All
T0 distribution of all channel.
unsigned short m_tdcMax
Upper boundary TDC histogram.
void initialize() override
Initializes the Module.
TH1D * m_hTDCBoard[300]
T0 distribution of each board.
bool m_storeFittedHisto
Store fitted histogram or not.
void event() override
Event action (main routine).
bool m_flag[56][400]
flag =1 for good, =0 for low statistic or bad fit
bool m_cosmic
for cosmic case, tof of upper sector will be negative
unsigned short m_adcMin
ADC cut to reject noise.
void terminate() override
Termination action, fit t0 and store histograms.
double m_t0[56][400]
T0 of each channel.
unsigned short m_tdcMin
Lower boundary TDC histogram.
double m_zOffset
z offset for calculate prop time, it is position of trigger counter,
double m_initT0
initial t0, use int fitting
std::string m_histoFileName
output file to store TDC histo after fit
unsigned short m_minEntries
min entries per histo.
TH1D * m_hTDC[56][400]
TDC distribution histo.
std::string m_outputFileName
output file name of t0 file.
The Class for CDC Geometry Parameters.
unsigned short getBoardID(const WireID &wID) const
Returns frontend board id. corresponding to the wire id.
const B2Vector3D wireBackwardPosition(uint layerId, int cellId, EWirePosition set=c_Base) const
Returns the backward position of the input sense wire.
double getTdcBinWidth() const
Return TDC bin width (nsec).
unsigned nWiresInLayer(int layerId) const
Returns wire numbers in a layer.
static CDCGeometryPar & Instance(const CDCGeometry *=nullptr)
Static method to get a reference to the CDCGeometryPar instance.
double senseWireR(int layerId) const
Returns radius of sense wire in each layer.
static const double speedOfLight
[cm/ns]
Definition Const.h:695
void setDescription(const std::string &description)
Sets the description of the module.
Definition Module.cc:214
void setPropertyFlags(unsigned int propertyFlags)
Sets the flags for the module properties.
Definition Module.cc:208
Module()
Constructor.
Definition Module.cc:30
@ c_ParallelProcessingCertified
This module can be run in parallel processing mode safely (All I/O must be done through the data stor...
Definition Module.h:80
Class to identify a wire inside the CDC.
Definition WireID.h:34
void addParam(const std::string &name, T &paramVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
Definition Module.h:559
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition Module.h:649
Abstract base class for different kinds of events.
STL namespace.