18 #include "TDirectory.h"
44 gPARMS->SetDefaultParameter(
"BCAL:const_term",
const_term);
53 gPARMS->SetDefaultParameter(
"DBCALShower:VERBOSE",
VERBOSE,
"Verbosity level for DBCALShower objects and factories");
54 gPARMS->SetDefaultParameter(
"DBCALShower:COVARIANCEFILENAME",
COVARIANCEFILENAME,
"File name for covariance files");
65 map<string, double> shower_calib2;
68 if (loop->GetCalib(
"BCAL/shower_calib2", shower_calib2)){
69 jerr <<
" Error loading BCAL nonlinear correction parameters from CCDB\n";
97 if (result!=NOERROR)
return result;
100 vector<const DBCALGeometry *> BCALGeomVec;
101 loop->Get(BCALGeomVec);
102 if(BCALGeomVec.size() == 0)
103 throw JException(
"Could not load DBCALGeometry object!");
112 for (
int i=0; i<5; i++) {
113 for (
int j=0; j<=i; j++) {
125 vector< const DBCALCluster* > clusters;
126 loop->Get( clusters );
132 for( vector< const DBCALCluster* >::const_iterator clItr = clusters.begin();
133 clItr != clusters.end();
136 if( isnan((**clItr).t()) == 1 || isnan((**clItr).theta()) == 1 || isnan((**clItr).phi()) == 1 || isnan((**clItr).rho()) == 1 )
continue;
138 float cosTh = cos( (**clItr).theta() );
139 float sinTh =
sin( (**clItr).theta() );
140 float cosPhi = cos( (**clItr).phi() );
141 float sinPhi =
sin( (**clItr).phi() );
142 float rho = (**clItr).rho();
143 if (
VERBOSE>2)
printf(
"%4lu cluster: E=%10.6f th=%10.6f phi=%10.6f rho=%10.6f t=%10.6f\n",eventnumber,
144 (**clItr).E(),(**clItr).theta()/3.14159265*180,(**clItr).phi()/3.14159265*180,(**clItr).rho(),(**clItr).t());
148 shower->E_raw = (**clItr).E();
149 shower->E_preshower = (**clItr).E_preshower();
150 shower->E_L2 = (**clItr).E_L2();
151 shower->E_L3 = (**clItr).E_L3();
152 shower->E_L4 = (**clItr).E_L4();
153 shower->x = rho * sinTh * cosPhi;
154 shower->y = rho * sinTh * sinPhi;
156 shower->Q = (**clItr).Q();
161 double t = (**clItr).t();
163 double dist_in_BCAL = rho - inner_rad/sinTh;
168 shower->sigLong = (**clItr).sigRho();
169 shower->sigTrans = (**clItr).sigPhi();
170 shower->sigTheta = (**clItr).sigTheta();
172 shower->rmsTime = (**clItr).rmsTime();
174 shower->N_cell = (**clItr).nCells();
183 printf(
"shower: E=%10.6f x=%10.6f y=%10.6f z=%10.6f t=%10.6f\n",
184 shower->E,shower->x,shower->y,shower->z,shower->t);
185 printf(
"shower: dE=%10.6f dx=%10.6f dy=%10.6f dz=%10.6f dt=%10.6f\n",
186 shower->EErr(),shower->xErr(),shower->yErr(),shower->zErr(),shower->tErr());
187 printf(
"shower: Ex=%10.6f Ey=%10.6f Ez=%10.6f Et=%10.6f xy=%10.6f\n",
188 shower->EXcorr(),shower->EYcorr(),shower->EZcorr(),shower->ETcorr(),shower->XYcorr());
189 printf(
"shower: xz=%10.6f xt=%10.6f yz=%10.6f yt=%10.6f zt=%10.6f\n\n",
190 shower->XZcorr(),shower->XTcorr(),shower->YZcorr(),shower->YTcorr(),shower->ZTcorr());
193 shower->AddAssociatedObject(*clItr);
195 _data.push_back( shower );
210 float minElookup = xaxis->GetBinLowEdge(1);
211 float maxElookup = xaxis->GetBinUpEdge(xaxis->GetNbins());
212 float minthlookup = yaxis->GetBinLowEdge(1);
213 float maxthlookup = yaxis->GetBinUpEdge(yaxis->GetNbins());
215 float shower_E = shower->
E;
216 float shower_r =
sqrt(shower->
x*shower->
x + shower->
y*shower->
y);
217 float shower_theta = atan2(shower_r,shower->
z-
m_zTarget);
218 float thlookup = shower_theta/3.14159265*180;
219 float Elookup = shower_E;
222 if (Elookup<minElookup) Elookup=minElookup;
223 if (Elookup>maxElookup) Elookup=maxElookup-0.0001;
224 if (thlookup<minthlookup) thlookup=minthlookup;
225 if (thlookup>maxthlookup) thlookup=maxthlookup-0.0001;
226 if (
VERBOSE>3)
printf(
"lookup (E,theta)=(%f,%F) limits (%f,%f) (%f,%f)\n",
227 Elookup,thlookup,minElookup,maxElookup,minthlookup,maxthlookup);
230 for (
int i=0; i<5; i++) {
231 for (
int j=0; j<=i; j++) {
233 if (i==0 && j==0) val *= shower_E;
234 ErphiztCovariance(i,j) = ErphiztCovariance(j,i) = val;
238 float shower_phi = atan2(shower->
y,shower->
x);
239 float cosPhi = cos(shower_phi);
240 float sinPhi =
sin(shower_phi);
242 rotationmatrix(0,0) = 1;
243 rotationmatrix(3,3) = 1;
244 rotationmatrix(4,4) = 1;
245 rotationmatrix(1,1) = cosPhi;
246 rotationmatrix(1,2) = -sinPhi;
247 rotationmatrix(2,1) = sinPhi;
248 rotationmatrix(2,2) = cosPhi;
250 if (
VERBOSE>3) {
printf(
"(E,r,phi,z,t) "); ErphiztCovariance.Print(); }
251 DMatrixDSym &D = ErphiztCovariance.Similarity(rotationmatrix);
252 for (
int i=0; i<5; i++) {
253 for (
int j=0; j<5; j++)
265 std::thread::id this_id = std::this_thread::get_id();
266 stringstream idstring;
268 if (
VERBOSE>0)
printf(
"DBCALShower_factory_IU::LoadCovarianceLookupTables(): Thread %s\n",idstring.str().c_str());
274 map<string,string> covariance_data;
277 if (eventLoop->GetJCalibration()->GetCalib(
"/BCAL/shower_covariance", covariance_data)) {
278 jerr <<
"Error loading /BCAL/shower_covariance !" << endl;
279 return ERROR_OPENING_EVENT_SOURCE;
281 if (covariance_data.size() == 15) {
284 for(
auto element : covariance_data) {
285 cout <<
"\nTEST: " << element.first <<
" = " << element.second << endl;
289 jerr <<
"Wrong number of elements /BCAL/shower_covariance !" << endl;
290 return ERROR_OPENING_EVENT_SOURCE;
294 for (
int i=0; i<5; i++) {
295 for (
int j=0; j<=i; j++) {
297 japp->RootWriteLock();
299 TDirectory *
savedir = gDirectory;
306 stringstream matrixname;
307 matrixname <<
"covmatrix_" << i << j;
308 if (
VERBOSE>1) cout <<
"Using CCDB \"" << matrixname.str() <<
"\" " << covariance_data[matrixname.str()] << endl;
309 ss.str(covariance_data[matrixname.str()]);
313 if (
VERBOSE>0) cout << filename << std::endl;
315 if (! ifs.is_open()) {
316 jerr <<
" Error: Cannot open file! " << filename << std::endl;
317 return ERROR_OPENING_EVENT_SOURCE;
319 getline(ifs, line,
'\n');
321 if (
VERBOSE>1) cout <<
"Using file " <<line<<endl;
329 Float_t xbins[nxbins+1];
330 Float_t ybins[nybins+1];
331 for (
int count=0; count<=nxbins; count++) {
336 for (
int count=0; count<=nybins; count++) {
346 sprintf(histname,
"covariance_%i%i_thread%s",i,j,idstring.str().c_str());
348 CovarianceLookupTable[i][j]->SetDirectory(
nullptr);
351 if (
VERBOSE>1)
printf(
"(%2i,%2i) (%2i,%2i) %e ",i,j,xbin,ybin,cont);
352 CovarianceLookupTable[i][j]->SetBinContent(xbin,ybin,cont);
354 if (ybin>nybins) { xbin++; ybin=1; }
const DBCALGeometry * dBCALGeom
jerror_t FillCovarianceMatrix(DBCALShower *shower)
TMatrixFSym ExyztCovariance
sprintf(text,"Post KinFit Cut")
double second_term_scale_factor
double first_term_scale_factor
double LOAD_CCDB_CONSTANTS
double second_exp_scale_factor
DGeometry * GetDGeometry(unsigned int run_number)
jerror_t LoadCovarianceLookupTables(JEventLoop *eventLoop)
float GetBCAL_inner_rad() const
TH2F * CovarianceLookupTable[5][5]
jerror_t brun(JEventLoop *loop, int32_t runnumber)
string COVARIANCEFILENAME
printf("string=%s", string)
double second_exp_const_term
bool GetTargetZ(double &z_target) const
z-location of center of target
jerror_t evnt(JEventLoop *loop, uint64_t eventnumber)