Skip to content

Commit 7e2a5e1

Browse files
aaustaaust
authored
Merge pull request #787 from JeffersonLab/ijaegle-nonlinear
Add CPP non-linear correction
2 parents 8af6d56 + 95e4f86 commit 7e2a5e1

File tree

2 files changed

+82
-7
lines changed

2 files changed

+82
-7
lines changed

‎src/libraries/FCAL/DFCALShower_factory.cc

Lines changed: 79 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -38,6 +38,9 @@ DFCALShower_factory::DFCALShower_factory()
3838

3939
USE_RING_E_CORRECTION_V2=false;
4040
gPARMS->SetDefaultParameter("FCAL:USE_RING_E_CORRECTION_V2",USE_RING_E_CORRECTION_V2);
41+
42+
USE_CPP_E_CORRECTION=false;
43+
gPARMS->SetDefaultParameter("FCAL:USE_CPP_CORRECTION",USE_CPP_E_CORRECTION);
4144

4245
SHOWER_ENERGY_THRESHOLD = 50*k_MeV;
4346
gPARMS->SetDefaultParameter("FCAL:SHOWER_ENERGY_THRESHOLD", SHOWER_ENERGY_THRESHOLD);
@@ -79,15 +82,22 @@ DFCALShower_factory::DFCALShower_factory()
7982
SHOWER_POSITION_LOG = true;
8083
USE_RING_E_CORRECTION_V1 = true;
8184
USE_RING_E_CORRECTION_V2 = false;
85+
USE_CPP_E_CORRECTION = false;
8286
} else if (USE_NONLINEAR_CORRECTION_TYPE == 2) {
8387
expfit_param1 = 2;
8488
expfit_param1 = 0;
8589
expfit_param1 = 0;
8690
SHOWER_POSITION_LOG = true;
8791
USE_RING_E_CORRECTION_V1 = false;
8892
USE_RING_E_CORRECTION_V2 = true;
93+
USE_CPP_E_CORRECTION = false;
94+
} else if (USE_NONLINEAR_CORRECTION_TYPE == 3) {
95+
SHOWER_POSITION_LOG = true;
96+
USE_RING_E_CORRECTION_V1 = false;
97+
USE_RING_E_CORRECTION_V2 = false;
98+
USE_CPP_E_CORRECTION = true;
8999
}
90-
100+
91101
gPARMS->SetDefaultParameter("FCAL:P0", timeConst0);
92102
gPARMS->SetDefaultParameter("FCAL:P1", timeConst1);
93103
gPARMS->SetDefaultParameter("FCAL:P2", timeConst2);
@@ -197,6 +207,7 @@ jerror_t DFCALShower_factory::brun(JEventLoop *loop, int32_t runnumber)
197207
jout << Form(" %s", str_coef[0].Data()) << nonlinear_correction_type[0];
198208
jout << endl;
199209
}
210+
200211
if (nonlinear_correction_type[0] == 0) {
201212
LOAD_NONLIN_CCDB = true;
202213
} else if (nonlinear_correction_type[0] == 1) {
@@ -207,6 +218,7 @@ jerror_t DFCALShower_factory::brun(JEventLoop *loop, int32_t runnumber)
207218
SHOWER_POSITION_LOG = true;
208219
USE_RING_E_CORRECTION_V1 = true;
209220
USE_RING_E_CORRECTION_V2 = false;
221+
USE_CPP_E_CORRECTION = false;
210222
} else if (nonlinear_correction_type[0] == 2) {
211223
LOAD_NONLIN_CCDB = true;
212224
expfit_param1 = 2;
@@ -215,12 +227,21 @@ jerror_t DFCALShower_factory::brun(JEventLoop *loop, int32_t runnumber)
215227
SHOWER_POSITION_LOG = true;
216228
USE_RING_E_CORRECTION_V1 = false;
217229
USE_RING_E_CORRECTION_V2 = true;
230+
USE_CPP_E_CORRECTION = false;
231+
} else if (nonlinear_correction_type[0] == 3) {
232+
LOAD_NONLIN_CCDB = true;
233+
SHOWER_POSITION_LOG = true;
234+
USE_RING_E_CORRECTION_V1 = false;
235+
USE_RING_E_CORRECTION_V2 = false;
236+
USE_CPP_E_CORRECTION = true;
218237
}
219238
}
220-
239+
221240
// but allow these to be overridden by command line parameters
222241
energy_dependence_correction_vs_ring.clear();
223242
nonlinear_correction.clear();
243+
nonlinear_correction_cpp.clear();
244+
block_to_square.clear();
224245
if(LOAD_NONLIN_CCDB) {
225246
map<string, double> shower_calib_piecewise;
226247
loop->GetCalib("FCAL/shower_calib_piecewise", shower_calib_piecewise);
@@ -232,7 +253,9 @@ jerror_t DFCALShower_factory::brun(JEventLoop *loop, int32_t runnumber)
232253
expfit_param3 = shower_calib_piecewise["expfit_param3"];
233254
m_beamSpotX = 0;
234255
m_beamSpotY = 0;
235-
256+
//expfit_param1 = 1.10358;
257+
//expfit_param2 = 0.31385;
258+
//expfit_param3 = -2.02585;
236259
if(debug_level>0) {
237260
jout << "cutoff_energy = " << cutoff_energy << endl;
238261
jout << "linfit_slope = " << linfit_slope << endl;
@@ -257,7 +280,7 @@ jerror_t DFCALShower_factory::brun(JEventLoop *loop, int32_t runnumber)
257280
}
258281
}
259282
loop->GetCalib("FCAL/nonlinear_correction", nonlinear_correction);
260-
if (nonlinear_correction.size() > 0 && nonlinear_correction[0][0] != 0) {
283+
if (nonlinear_correction.size() > 0) {
261284
m_beamSpotX = beam_spot.at("x");
262285
m_beamSpotY = beam_spot.at("y");
263286
if (debug_level > 0) {
@@ -271,8 +294,29 @@ jerror_t DFCALShower_factory::brun(JEventLoop *loop, int32_t runnumber)
271294
}
272295
}
273296
}
297+
298+
loop->GetCalib("FCAL/block_to_square", block_to_square);
299+
if (block_to_square.size() > 0) {
300+
if (debug_level > 0) {
301+
for (int i = 0; i < (int) block_to_square.size(); i ++) {
302+
jout << block_to_square[i];
303+
}
304+
jout << endl;
305+
}
306+
}
307+
308+
loop->GetCalib("FCAL/nonlinear_correction_cpp", nonlinear_correction_cpp);
309+
if (nonlinear_correction_cpp.size() > 0) {
310+
m_beamSpotX = beam_spot.at("x");
311+
m_beamSpotY = beam_spot.at("y");
312+
if (debug_level > 0) {
313+
for (int i = 0; i < (int) nonlinear_correction_cpp.size(); i ++) {
314+
jout << nonlinear_correction_cpp[i];
315+
}
316+
jout << endl;
317+
}
318+
}
274319
}
275-
276320
if (LOAD_TIMING_CCDB) {
277321
// Get timing correction polynomial, J. Mirabelli 10/31/17
278322
map<string,double> timing_correction;
@@ -515,7 +559,11 @@ jerror_t DFCALShower_factory::evnt(JEventLoop *eventLoop, uint64_t eventnumber)
515559
//int MAXITER = 1000;
516560

517561
DVector3 posInCal = cluster->getCentroid();
518-
562+
int block = cluster->getChannelEmax();
563+
int square_nb = -1;
564+
if (USE_CPP_E_CORRECTION)
565+
square_nb = block_to_square[block];
566+
519567
float x0 = posInCal.Px();
520568
float y0 = posInCal.Py();
521569
double Eclust = cluster->getEnergy();
@@ -533,6 +581,8 @@ jerror_t DFCALShower_factory::evnt(JEventLoop *eventLoop, uint64_t eventnumber)
533581
double Egamma = Eclust;
534582
Ecorrected = 0;
535583

584+
585+
536586
// block properties
537587
double radiation_length=FCAL_RADIATION_LENGTH;
538588
double shower_offset=FCAL_SHOWER_OFFSET;
@@ -632,7 +682,29 @@ jerror_t DFCALShower_factory::evnt(JEventLoop *eventLoop, uint64_t eventnumber)
632682
// if all C=D=E=0 by mistake then Egamma = - Eclust
633683
Egamma = Eclust / (C - exp(-D * Eclust + E)); // Non-linear part
634684
}
635-
} // End Correction method I
685+
//cout <<"Eclust " << Eclust << " Egamma " << Egamma << " A " << A << " B " << B << " C " << C << " D " << D << " E " << E << endl;
686+
} // End Correction method I
687+
688+
if (USE_CPP_E_CORRECTION && !USE_RING_E_CORRECTION_V2 && !USE_RING_E_CORRECTION_V1) {
689+
double scalef = nonlinear_correction_cpp[0];
690+
if (square_nb >= 0 && square_nb <= 13) {
691+
double Eshift = 0;
692+
if (square_nb == 10) {
693+
Eshift = atan(nonlinear_correction_cpp[4 + square_nb * 3] * Egamma + nonlinear_correction_cpp[5 + square_nb * 3]);
694+
Eshift *= nonlinear_correction_cpp[3 + square_nb * 3] * Eshift;
695+
} else {
696+
Eshift = nonlinear_correction_cpp[3 + square_nb * 3] * atan(nonlinear_correction_cpp[4 + square_nb * 3] * Egamma + nonlinear_correction_cpp[5 + square_nb * 3]);
697+
}
698+
Eshift = scalef * Eshift;
699+
if (Eshift > 0.) {
700+
Egamma *= (1. + nonlinear_correction_cpp[1] * 1.e-2 * Egamma + nonlinear_correction_cpp[2] * 1.e-2 * Egamma * Egamma) / Eshift;
701+
} else {
702+
if (VERBOSE > 3) jerr << "CPP nonlinear correction has a wrong Eshift" << endl;
703+
}
704+
} else {
705+
if (VERBOSE > 3) jerr << "CPP nonlinear correction has no square_nb" << endl;
706+
}
707+
}
636708
}
637709
//End energy dependence correction
638710

‎src/libraries/FCAL/DFCALShower_factory.h

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -68,6 +68,7 @@ class DFCALShower_factory:public JFactory<DFCALShower>{
6868
bool LOAD_CCDB_CONSTANTS;
6969
bool USE_RING_E_CORRECTION_V1;
7070
bool USE_RING_E_CORRECTION_V2;
71+
bool USE_CPP_E_CORRECTION;
7172
bool SHOWER_POSITION_LOG;
7273
double SHOWER_ENERGY_THRESHOLD;
7374
double cutoff_energy;
@@ -81,6 +82,8 @@ class DFCALShower_factory:public JFactory<DFCALShower>{
8182
vector < vector <double> > energy_dependence_correction_vs_ring;
8283
//V2 based on E/P
8384
vector < vector <double> > nonlinear_correction;
85+
vector <double> nonlinear_correction_cpp;
86+
vector <int> block_to_square;
8487
//Choose directly in DB which non-linear correction version to use
8588
vector <int> nonlinear_correction_type;
8689

0 commit comments

Comments
 (0)