 |
Belle II Software
release-05-01-25
|
15 #include "trg/ecl/modules/trgeclTimingCal/TrgEclTimingCalibration.h"
17 #include "TDecompLU.h"
33 "0 : Raw data 1 : Simulation data", 0);
36 "0 : Beam data 1 : Cosmic ray data", 0);
39 "0 : Not using 3D bhabha veto bit 1 : Using 3D bhabha veto bit ", 1);
42 "Reference TCID. Default setting : 184", 184);
45 "A TC energy calibration factor (GeV/ADC), default 0.00525", 0.00525);
48 "Number of TC cut : Events will be used < ntc", 999);
51 "High energy cut", 9999.0);
54 "Low energy cut", 0.0);
57 "Flag including FWD end cap", 1);
60 "Flag including Barrel", 1);
63 "Flag including BWD end cap", 1);
72 "iteration chisq cut", 1000.0);
98 TCHitUnpacker.isRequired();
100 TRGECLEvtStore.isRequired();
136 if (nEvtStore <= 0)
return;
137 if (TRGECLEvtStore[0] -> getCheckSum() != 0)
return;
138 if (TRGECLEvtStore[0] -> get3DBhabhaV() != 1)
return;
144 if (nEvtStore <= 0)
return;
145 if (
TRGECLTrg[0] -> get3DBhabha() != 1)
return;
165 if (ntc < 2 || ntc >
cut_ntc)
return 0;
166 for (
int itc = 0; itc < ntc; itc++) {
167 TCId.push_back(TCHit[itc]->getTCId());
168 TCEnergy.push_back(TCHit[itc]->getEnergyDep());
169 TCTiming.push_back((
double) TCHit[itc]->getTimeAve());
174 if (ntc < 2 || ntc >
cut_ntc)
return 0;
175 for (
int itc = 0; itc < ntc; itc++) {
176 if (TCHit[itc] -> getChecksum() != 0)
return 0;
177 if (!(TCHit[itc] -> getHitWin() == 3 || TCHit[itc] -> getHitWin() == 4))
continue;
178 TCId.push_back(TCHit[itc]->getTCId());
180 TCTiming.push_back((
double) TCHit[itc]->getTCTime());
189 int flag_event_used = 0;
190 int ntc =
TCId.size();
192 for (
int itc1 = 0; itc1 < ntc; itc1++) {
195 int m_TCID1 =
TCId[itc1] - 1;
198 for (
int itc2 = itc1 + 1; itc2 < ntc; itc2++) {
199 int m_TCID2 =
TCId[itc2] - 1;
205 double sigma2 = 5.0 * 5.0;
210 _Matrix[m_TCID1][m_TCID2] += (-1.0 / sigma2);
211 _Matrix[m_TCID2][m_TCID1] += (-1.0 / sigma2);
213 _Matrix[m_TCID1][m_TCID1] += (1.0 / sigma2);
214 _Matrix[m_TCID2][m_TCID2] += (1.0 / sigma2);
217 Chisq1[m_TCID1][m_TCID2] += (1.0 / sigma2);
218 Chisq1[m_TCID2][m_TCID1] += (1.0 / sigma2);
222 double Vector_component = 0;
233 Vector_component /= sigma2;
234 _Vector[m_TCID1] += Vector_component;
235 _Vector[m_TCID2] += -Vector_component;
237 Chisq2[m_TCID1][m_TCID2] += -2.0 * Vector_component;
238 Chisq2[m_TCID2][m_TCID1] += 2.0 * Vector_component;
240 Chisq3[m_TCID1][m_TCID2] += Vector_component * Vector_component * sigma2;
241 Chisq3[m_TCID2][m_TCID1] += Vector_component * Vector_component * sigma2;
250 return flag_event_used;
259 const int nTC_tot = 576;
260 TMatrixDSym mat_A(575);
261 mat_A.ResizeTo(575, 575, 0.0);
266 for (
int i = 0; i < 576; i++) {
267 if (i ==
TC_ref - 1)
continue;
269 if (i >
TC_ref - 1) indexi = i - 1;
271 for (
int j = 0; j < 576; j++) {
272 if (j ==
TC_ref - 1)
continue;
274 if (j >
TC_ref - 1) indexj = j - 1;
275 mat_A[indexi][indexj] =
_Matrix[i][j];
291 for (
int i = 0; i < 576; i++) {
300 TVectorD x = lu.Solve(vec_b, b_solve);
304 for (
int i = 0; i < nTC_tot; i++) {
308 }
else if (i ==
TC_ref - 1) {
314 for (
int i = 0; i < nTC_tot; i++) {
324 for (
int iTCID = 0; iTCID < 576; iTCID++) {
341 for (
int i = 0; i < 576; i++) {
347 _Matrix.ResizeTo(576, 576, 0.0);
349 for (
int i = 0; i < 576; i++) {
351 for (
int j = 0; j < 576; j++) {
360 for (
int i = 0; i < 576; i++) {
361 Chisq1[i].assign(576, 0.0);
362 Chisq2[i].assign(576, 0.0);
363 Chisq3[i].assign(576, 0.0);
368 for (
int i = 0; i < 576; i++) {
373 for (
int i = 0; i < 80; i++) {
378 for (
int i = 80; i < 512; i++) {
383 for (
int i = 512; i < 576; i++) {
395 TFile* tf =
new TFile(fname,
"RECREATE");
396 TTree* tr =
new TTree(
"tree",
"TC time calibration result");
397 vector<double> time_offset;
398 vector<double> time_offset_err;
400 vector<double> chisq;
402 time_offset.resize(576);
403 time_offset_err.resize(576);
407 double timeoffsetave = 0.0;
408 tr -> Branch(
"TCID", &TCID);
409 tr -> Branch(
"TimeOffset", &time_offset);
410 tr -> Branch(
"TimeOffAve", &timeoffsetave);
411 tr -> Branch(
"TimeOffsetErr", &time_offset_err);
412 tr -> Branch(
"Chisq", &chisq);
413 tr -> Branch(
"Nevent", &Nevent);
415 tr -> Branch(
"Matrix", &
_Matrix);
416 tr -> Branch(
"Vector", &
_Vector);
417 tr -> Branch(
"ChisqComponent1", &
Chisq1);
418 tr -> Branch(
"ChisqComponent2", &
Chisq2);
419 tr -> Branch(
"ChisqComponent3", &
Chisq3);
422 for (
int i = 0; i < 576; i++) {
430 timeoffsetave /= 576.0;
451 for (
int i = 0; i < 576; i++) {
453 std::vector<double> offsetdiff;
454 offsetdiff.assign(576, 0.0);
455 for (
int j = 0; j < 576; j++) {
456 if (i == j || (
b_Inc_TC[j] ==
false && j !=
TC_ref - 1))
continue;
462 for (
int i = 0; i < 576; i++) {
472 TFile* tf =
new TFile(str_fname_tmp,
"READ");
473 auto tr = (TTree*) tf -> Get(
"tree");
474 vector<double>* TimeOffset_tmp = 0;
475 tr -> SetBranchAddress(
"TimeOffset", &TimeOffset_tmp);
std::vector< std::vector< double > > Chisq2
chisq component2
int fInclude_BR
Flag include barrel 0 : exclude BR 1 : include BR.
int fCalType
Calibration type 0 : beam 1 : cosmic.
virtual void beginRun() override
Begin Run function.
std::vector< double > chisq_result
Chisq values based on time offset obtained from "Solve_Matrix" function.
virtual void initialize() override
initialize function
void setDescription(const std::string &description)
Sets the description of the module.
double cut_chisq
TC chisq cut (iteration mode)
int Solve_Matrix()
Solve matrix.
int cut_ntc
The number of TC cut.
void Save_Result()
Save time offset and chisq as a root file.
void InitParams()
Parameters initialization.
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
@ c_ParallelProcessingCertified
This module can be run in parallel processing mode safely (All I/O must be done through the data stor...
std::vector< bool > b_Inc_TC
Flag TC included.
std::vector< int > Nevent_TC
The number of uesed TC hit (TC by TC)
int fInclude_FWD
Flag include forward endcap 0 : exclude FWD 1 : include FWD.
std::vector< double > TimeOffset
Time offset Input time offset (iteration mode)
std::vector< double > tcal_result
Time offset calibration result.
std::string str_default_name
Default output root file name.
double cut_low_energy
TC energy cut lower limit.
int fInclude_BWD
Flag include backward endcap 0 : exclude BWD 1 : include BWD.
int TC_ref
Reference TC whose time offset set to be 0.
double cut_high_energy
TC energy cut upper limit.
virtual void event() override
Event function.
int f3Dbhabha_veto
Flag 3D bhabha veto 0 : not use 1 : use.
std::vector< std::vector< double > > Chisq1
chisq component1
TVectorD _Vector
Vector component.
std::vector< double > TCTiming
TC time.
void setPropertyFlags(unsigned int propertyFlags)
Sets the flags for the module properties.
std::vector< std::vector< double > > Chisq3
chisq component3
void TRGECLTimingCalClear()
Clear vector every event.
int fSimulation
Flag simulation 0 : real data 1 : simulation data.
void GetTimeOffset()
Get previous time offsets (iteration mode)
Abstract base class for different kinds of events.
std::vector< TVector3 > TCPosition
TC position from TRGECLMap (cosmic calibration)
double TCEnergyCalibrationConstant
ADC to GeV energy conversion factor.
TVector3 getTCPosition(int)
TC position (cm)
int FlagMatrixSolved
Flag matrix solved (0 : not solved 1 : solved)
virtual void endRun() override
End Run function.
int fIteration
Flag iteration 0 : first calbiration 1 : iteration.
TMatrixDSym _Matrix
Matrix Marix component.
TRGECLTimingCalModule()
Constructor.
std::string str_timeoffset_fname
Input time offset file name (iteration mode)
virtual void terminate() override
Terminate function.
void addParam(const std::string &name, T ¶mVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
void Set_TCposition()
Set TC position from TRGECLMap.
int nevt_selected
The number of selected events.
Accessor to arrays stored in the data store.
int FillMatrix()
Fill matrix and vector components.
std::vector< int > TCId
Trigger Cell TCID.
std::string str_ofilename
Output root file name.
void Calculate_Chisq()
Calculate chisq based on time offset obtained from "Solve_Matrix" function.
std::vector< bool > b_Exclude_TC
Flag TC excluded.
int getEntries() const
Get the number of objects in the array.
int nevt_read
The number of read events.
int ReadData()
Fill data from.
std::vector< double > tcal_result_err
Time offset error.
std::vector< double > TCEnergy
TC energy.
virtual ~TRGECLTimingCalModule()
Destructor.