Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DBCALShower_factory_IU.cc
Go to the documentation of this file.
1 /*
2  * DBCALShower_factory_IU.cc
3  * (formerly DBCALShower_factory.cc)
4  *
5  * Created by Matthew Shepherd on 3/24/11.
6  *
7  */
8 
10 #include "DBCALCluster.h"
11 
12 #include "DANA/DApplication.h"
13 
14 #include "units.h"
15 
16 #include "TH2F.h"
17 #include "TROOT.h"
18 #include "TDirectory.h"
19 
20 #include <thread>
21 
22 #include <DMatrix.h>
23 #include <DMatrixDSym.h>
24 
25 
27  // defaults set to minimize probles if new values are not loaded
29  const_term = 1; ///< default to make no change to energy
30  first_term_scale_factor = 0; ///< default to make no change to energy
31  first_exp_param0 = 1; ///< default to make no change to energy
32  first_exp_param1 = 1; ///< default to make no change to energy
33  second_term_scale_factor = 0;///< default to make no change to energy
34  second_exp_const_term = 1; ///< default to make no change to energy
35  second_exp_scale_factor = 1; ///< default to make no change to energy
36  second_exp_param0 = 1; ///< default to make no change to energy
37  second_exp_param1 = 1; ///< default to make no change to energy
38  VERBOSE = 0; ///< >0 once off info ; >2 event by event ; >3 everything
39  COVARIANCEFILENAME = ""; ///< Setting the filename will take precidence over the CCDB. Files must end in ij.txt, where i and j are integers corresponding to the element of the matrix
40 
41  if (gPARMS){
42  gPARMS->SetDefaultParameter("BCAL:LOAD_NONLIN_CCDB", LOAD_CCDB_CONSTANTS);
43  /// use to set energy corrections on command line
44  gPARMS->SetDefaultParameter("BCAL:const_term", const_term);
45  gPARMS->SetDefaultParameter("BCAL:first_term_scale_factor", first_term_scale_factor);
46  gPARMS->SetDefaultParameter("BCAL:first_exp_param0", first_exp_param0);
47  gPARMS->SetDefaultParameter("BCAL:first_exp_param1", first_exp_param1);
48  gPARMS->SetDefaultParameter("BCAL:second_term_scale_factor", second_term_scale_factor);
49  gPARMS->SetDefaultParameter("BCAL:second_exp_const_term", second_exp_const_term);
50  gPARMS->SetDefaultParameter("BCAL:second_exp_scale_factor", second_exp_scale_factor);
51  gPARMS->SetDefaultParameter("BCAL:second_exp_param0", second_exp_param0);
52  gPARMS->SetDefaultParameter("BCAL:second_exp_param1", second_exp_param1);
53  gPARMS->SetDefaultParameter("DBCALShower:VERBOSE", VERBOSE, "Verbosity level for DBCALShower objects and factories");
54  gPARMS->SetDefaultParameter("DBCALShower:COVARIANCEFILENAME", COVARIANCEFILENAME, "File name for covariance files");
55  }
56 }
57 
58 jerror_t DBCALShower_factory_IU::brun(JEventLoop *loop, int32_t runnumber) {
59  DApplication* app = dynamic_cast<DApplication*>(loop->GetJApplication());
60  DGeometry* geom = app->GetDGeometry(runnumber);
61  geom->GetTargetZ(m_zTarget);
62 
63  //by default, energy correction parameters are obtained through ccdb
64  if(LOAD_CCDB_CONSTANTS > 0.5){
65  map<string, double> shower_calib2;
66 
67 
68  if (loop->GetCalib("BCAL/shower_calib2", shower_calib2)){
69  jerr << " Error loading BCAL nonlinear correction parameters from CCDB\n";
70  } else {
71  const_term = shower_calib2["const_term"];
72  first_term_scale_factor = shower_calib2["first_term_scale_factor"];
73  first_exp_param0 = shower_calib2["first_exp_param0"];
74  first_exp_param1 = shower_calib2["first_exp_param1"];
75  second_term_scale_factor = shower_calib2["second_term_scale_factor"];
76  second_exp_const_term = shower_calib2["second_exp_const_term"];
77  second_exp_scale_factor = shower_calib2["second_exp_scale_factor"];
78  second_exp_param0 = shower_calib2["second_exp_param0"];
79  second_exp_param1 = shower_calib2["second_exp_param1"];
80  }
81 
82  }
83 
84  if(VERBOSE>0) {
85  printf("%20s = %f\n","const_term",const_term);
86  printf("%20s = %f\n","first_term_scale_factor",first_term_scale_factor);
87  printf("%20s = %f\n","first_exp_param0",first_exp_param0);
88  printf("%20s = %f\n","first_exp_param1",first_exp_param1);
89  printf("%20s = %f\n","second_term_scale_factor",second_term_scale_factor);
90  printf("%20s = %f\n","second_exp_const_term",second_exp_const_term);
91  printf("%20s = %f\n","second_exp_scale_factor",second_exp_scale_factor);
92  printf("%20s = %f\n","second_exp_param0",second_exp_param0);
93  printf("%20s = %f\n","second_exp_param1",second_exp_param1);
94  }
95 
96  jerror_t result = LoadCovarianceLookupTables(loop);
97  if (result!=NOERROR) return result;
98 
99  // load BCAL geometry
100  vector<const DBCALGeometry *> BCALGeomVec;
101  loop->Get(BCALGeomVec);
102  if(BCALGeomVec.size() == 0)
103  throw JException("Could not load DBCALGeometry object!");
104  dBCALGeom = BCALGeomVec[0];
105 
106  return NOERROR;
107 }
108 
109 
111  // delete lookup tables to prevent memory leak
112  for (int i=0; i<5; i++) {
113  for (int j=0; j<=i; j++) {
114  delete CovarianceLookupTable[i][j];
115  CovarianceLookupTable[i][j] = nullptr;
116  }
117  }
118  return NOERROR;
119 }
120 
121 
122 jerror_t
123 DBCALShower_factory_IU::evnt( JEventLoop *loop, uint64_t eventnumber ){
124 
125  vector< const DBCALCluster* > clusters;
126  loop->Get( clusters );
127 
128  // loop through and fill the shower structure from the cluster
129  // right now just a simple 1 to 1 correspondence with
130  // an overall energy correction
131 
132  for( vector< const DBCALCluster* >::const_iterator clItr = clusters.begin();
133  clItr != clusters.end();
134  ++clItr ){
135 
136  if( isnan((**clItr).t()) == 1 || isnan((**clItr).theta()) == 1 || isnan((**clItr).phi()) == 1 || isnan((**clItr).rho()) == 1 ) continue;
137 
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());
145 
146  DBCALShower* shower = new DBCALShower();
147 
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;
155  shower->z = rho * cosTh + m_zTarget;
156  shower->Q = (**clItr).Q();
157 
158  //DBCALCluster::t() returns the time at the inner radius
159  //so we need to make an adjustment so that the shower t is the time at
160  //the shower location (x,y,z)
161  double t = (**clItr).t();
162  double inner_rad = dBCALGeom->GetBCAL_inner_rad();
163  double dist_in_BCAL = rho - inner_rad/sinTh;
164  t = t + dist_in_BCAL/(30*k_cm/k_nsec);
165  shower->t = t;
166 
167  // shower widths for further selection in REST
168  shower->sigLong = (**clItr).sigRho();
169  shower->sigTrans = (**clItr).sigPhi();
170  shower->sigTheta = (**clItr).sigTheta();
171 // shower->sigTime = (**clItr).sigT();
172  shower->rmsTime = (**clItr).rmsTime();
173 
174  shower->N_cell = (**clItr).nCells();
175 
176  // Non-linear energy corrections can be found at https://logbooks.jlab.org/entry/3469359
177 
179 
180  // Get covariance matrix and uncertainties
181  FillCovarianceMatrix(shower);
182  if (VERBOSE>2) {
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());
191  }
192 
193  shower->AddAssociatedObject(*clItr);
194 
195  _data.push_back( shower );
196  }
197 
198  return NOERROR;
199 }
200 
201 
202 jerror_t
204  /// This function takes a BCALShower object and using the internal variables
205  /// overwrites any existing covaraince matrix using lookup tables.
206 
207  // Get edges of lookup table histograms (assume that all histograms have the same limits.)
208  TAxis *xaxis = CovarianceLookupTable[0][0]->GetXaxis();
209  TAxis *yaxis = CovarianceLookupTable[0][0]->GetYaxis();
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());
214 
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;
220 
221  // Adjust values: in order to use Interpolate() must be within histogram range
222  if (Elookup<minElookup) Elookup=minElookup;
223  if (Elookup>maxElookup) Elookup=maxElookup-0.0001; // move below edge, on edge doesn't work.
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);
228 
229  DMatrixDSym ErphiztCovariance(5);
230  for (int i=0; i<5; i++) {
231  for (int j=0; j<=i; j++) {
232  float val = CovarianceLookupTable[i][j]->Interpolate(Elookup, thlookup);
233  if (i==0 && j==0) val *= shower_E; // E variance is divided by energy in CCDB
234  ErphiztCovariance(i,j) = ErphiztCovariance(j,i) = val;
235  }
236  }
237 
238  float shower_phi = atan2(shower->y,shower->x);
239  float cosPhi = cos(shower_phi);
240  float sinPhi = sin(shower_phi);
241  DMatrix rotationmatrix(5,5);
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;
249 
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++)
254  shower->ExyztCovariance(i, j) = D(i, j);
255  }
256  if (VERBOSE>2) {printf("(E,x,y,z,t) "); shower->ExyztCovariance.Print(); }
257 
258  return NOERROR;
259 }
260 
261 
262 jerror_t
264  // Note that there's no error checking that the lookup tables have been loaded correctly!!
265  std::thread::id this_id = std::this_thread::get_id();
266  stringstream idstring;
267  idstring << this_id;
268  if (VERBOSE>0) printf("DBCALShower_factory_IU::LoadCovarianceLookupTables(): Thread %s\n",idstring.str().c_str());
269 
270  bool USECCDB=0;
271  // if filename specified try to use filename else get info from CCDB
272  if (COVARIANCEFILENAME == "") USECCDB=1;
273 
274  map<string,string> covariance_data;
275  if (USECCDB) {
276  // load information for covariance matrix
277  if (eventLoop->GetJCalibration()->GetCalib("/BCAL/shower_covariance", covariance_data)) {
278  jerr << "Error loading /BCAL/shower_covariance !" << endl;
279  return ERROR_OPENING_EVENT_SOURCE;
280  }
281  if (covariance_data.size() == 15) { // there are 15 elements in the covariance matrix
282  // for example, print it all out
283  if (VERBOSE>0) {
284  for(auto element : covariance_data) {
285  cout << "\nTEST: " << element.first << " = " << element.second << endl;
286  }
287  }
288  } else {
289  jerr << "Wrong number of elements /BCAL/shower_covariance !" << endl;
290  return ERROR_OPENING_EVENT_SOURCE;
291  }
292  }
293 
294  for (int i=0; i<5; i++) {
295  for (int j=0; j<=i; j++) {
296 
297  japp->RootWriteLock();
298  // change directory to memory so that histograms are not saved to file
299  TDirectory *savedir = gDirectory;
300 
301  // Read in string
302  ifstream ifs;
303  string line;
304  stringstream ss;
305  if (USECCDB) {
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()]);
310  } else {
311  char filename[255];
312  sprintf(filename,"%s%i%i.txt",COVARIANCEFILENAME.c_str(),i,j);
313  if (VERBOSE>0) cout << filename << std::endl;
314  ifs.open(filename);
315  if (! ifs.is_open()) {
316  jerr << " Error: Cannot open file! " << filename << std::endl;
317  return ERROR_OPENING_EVENT_SOURCE;
318  }
319  getline(ifs, line, '\n');
320  ss.str(line);
321  if (VERBOSE>1) cout << "Using file " <<line<<endl;
322  }
323 
324  // Parse string
325  int nxbins, nybins;
326  ss>>nxbins;
327  ss>>nybins;
328  if (VERBOSE>1) printf("bins (%i,%i)\n",nxbins,nybins);
329  Float_t xbins[nxbins+1];
330  Float_t ybins[nybins+1];
331  for (int count=0; count<=nxbins; count++) {
332  ss>>xbins[count];
333  if (VERBOSE>1) printf("(%i,%f) ",count,xbins[count]);
334  }
335  if (VERBOSE>1) printf("\n");
336  for (int count=0; count<=nybins; count++) {
337  ss>>ybins[count];
338  if (VERBOSE>1) printf("(%i,%f) ",count,ybins[count]);
339  }
340  if (VERBOSE>1) printf("\n");
341  int xbin=1;
342  double cont;
343  int ybin=1;
344  // create histogram
345  char histname[255];
346  sprintf(histname,"covariance_%i%i_thread%s",i,j,idstring.str().c_str());
347  CovarianceLookupTable[i][j] = new TH2F(histname,"Covariance histogram",nxbins,xbins,nybins,ybins);
348  CovarianceLookupTable[i][j]->SetDirectory(nullptr);
349  // fill histogram
350  while(ss>>cont){
351  if (VERBOSE>1) printf("(%2i,%2i) (%2i,%2i) %e ",i,j,xbin,ybin,cont);
352  CovarianceLookupTable[i][j]->SetBinContent(xbin,ybin,cont);
353  ybin++;
354  if (ybin>nybins) { xbin++; ybin=1; }
355  }
356  if (VERBOSE>1) printf("\n");
357  // Close file
358  ifs.close();
359 
360  savedir->cd();
361  japp->RootUnLock();
362 
363  }
364  }
365  return NOERROR;
366 }
367 
const DBCALGeometry * dBCALGeom
TMatrixD DMatrix
Definition: DMatrix.h:14
jerror_t FillCovarianceMatrix(DBCALShower *shower)
TMatrixFSym ExyztCovariance
Definition: DBCALShower.h:34
sprintf(text,"Post KinFit Cut")
TString filename
TMatrixDSym DMatrixDSym
Definition: DMatrixDSym.h:13
JApplication * japp
const double k_cm
Definition: units.h:14
DGeometry * GetDGeometry(unsigned int run_number)
jerror_t LoadCovarianceLookupTables(JEventLoop *eventLoop)
float GetBCAL_inner_rad() const
jerror_t brun(JEventLoop *loop, int32_t runnumber)
double sqrt(double)
double sin(double)
const double k_nsec
Definition: units.h:67
TDirectory * savedir
printf("string=%s", string)
bool GetTargetZ(double &z_target) const
z-location of center of target
Definition: DGeometry.cc:1933
jerror_t evnt(JEventLoop *loop, uint64_t eventnumber)