## Attachment 'CalibProcessor.cc'

```   1 #include "CalibProcessor.h"
2
3 using namespace lcio ;
4 using namespace marlin ;
5
6 /* **** ************************************* **** */
7 /* ****  Additional Structures and Functions  **** */
8 /* ****  Needed to explain a Gauss fit to gsl **** */
9 /* **** ************************************* **** */
10
11
12 struct data {
13   int n;
14   double* x;
15   double* y;
16 };
17
18
19
20 int GausFunc(const gsl_vector* par, void* d, gsl_vector* f) {
21
22 /* ------------------------------------------------ */
23 /* Function to fit:                                 */
24 /* g = norm \ (sqrt(2*pi)*sigma )                   */
25 /*   * exp [- (pow((x-mean),2)) / (2*pow(sigma,2))] */
26 /*                                                  */
27 /* Function to minimise:                            */
28 /* f = norm \ (sqrt(2*pi)*sigma )                   */
29 /*   * exp [- (pow((x-mean),2)) / (2*pow(sigma,2))] */
30 /*   - g                                            */
31 /* ------------------------------------------------ */
32
33   double mean  = gsl_vector_get(par,0);
34   double sigma = gsl_vector_get(par,1);
35   double norm  = gsl_vector_get(par,2);
36
37   int n        = ((struct data*)d)->n;
38   double* t    = ((struct data*)d)->x;
39   double* g    = ((struct data*)d)->y;
40   double fit   = 0.0;
41
42   for (int i(0); i < n; i++) {
43     fit =  norm / (sqrt(2*M_PI*pow(sigma,2) ) ) * exp (- (pow((t[i]-mean),2))/(2*pow(sigma,2)) )
44 	   - g[i];
45     gsl_vector_set(f,i,fit);
46   }
47
48   return GSL_SUCCESS;
49 }
50
51
52 int dGausFunc(const gsl_vector* par, void* d,  gsl_matrix* J) {
53
54   double mean  = gsl_vector_get(par,0);
55   double sigma = gsl_vector_get(par,1);
56   double norm = gsl_vector_get(par,2);
57
58   int n       = ((struct data*)d)->n;
59   double* t    = ((struct data*)d)->x;
60
61   // calculate Jacobi's matrix J[i][0] = dGaus/dMean, J[i][1] = dGaus/dSigma, J[i][2] = dGaus/dNormalisation
62   for (int i(0); i < n; i++) {
63     gsl_matrix_set(J,i,0,( ( norm*(t[i]-mean)/sqrt(2*M_PI*pow(sigma,6))
64 			     * exp (- (pow((t[i]-mean),2))/(2*pow(sigma,2))) ) ) );
65     gsl_matrix_set(J,i,1,( (norm/sqrt(2*M_PI*pow(sigma,2))
66 			    * exp(-pow(t[i]-mean,2)/(2*pow(sigma,2)))
67 			    * ((pow((t[i]-mean),2)/(pow(sigma,3))) - 1/sigma) ) ) );
68     gsl_matrix_set(J,i,2,( (1/sqrt(2*M_PI*pow(sigma,2))
69 			    * exp(-pow(t[i]-mean,2)/(2*pow(sigma,2))) ) ) );
70   }
71
72   return GSL_SUCCESS;
73 }
74
75 int fdfGausFunc(const gsl_vector* par, void* d, gsl_vector* f, gsl_matrix* J) {
76
77   GausFunc(par, d, f);
78   dGausFunc(par, d, J);
79
80   return GSL_SUCCESS;
81
82 }
83
84
85 /* **** ****************************** **** */
86 /* **** Start of Calibration Processor **** */
87 /* **** ****************************** **** */
88
89
90 CalibProcessor aCalibProcessor ;
91
92
93 CalibProcessor::CalibProcessor() : Processor("CalibProcessor")
94 {
95   // processor description
96   _description = "CalibProcessor calculates new sampling fraction coefficients c1, c2, c3 needed for energy recalibration of the ECAL and HCAL." ;
97
98   registerProcessorParameter( "colNameHCALHits" ,
99 			      "Name of the HCAL hit collection" ,
100 			      _colNameHCALHits ,
101 			      std::string("HCAL") ) ;
102
103   registerProcessorParameter( "colNameECALHits" ,
104 			      "Name of the ECAL hit collection" ,
105 			      _colNameECALHits ,
106 			      std::string("ECAL") ) ;
107
108
109   registerProcessorParameter("CMS_Energy" ,
110 			     "CMS energy [GeV]" ,
111 			     _E_CMS ,
112 			     float(500) );
113
114   registerProcessorParameter("N_Bins" ,
115 			     "number of bins" ,
116 			     _Nbins ,
117 			     int(100) );
118
119
120   std::vector<float> defaultValuesForCalibrECAL;
121   defaultValuesForCalibrECAL.push_back(1.0);
122   defaultValuesForCalibrECAL.push_back(2.0);
123
124   registerProcessorParameter("CalibrECAL" ,
125 			     "CalibCoeffEcal" ,
126 			     _CalibrECAL ,
127 			     defaultValuesForCalibrECAL);
128
129
130   std::vector<float> defaultValuesForCalibrHCAL;
131   defaultValuesForCalibrHCAL.push_back(3.0);
132
133   registerProcessorParameter("CalibrHCAL" ,
134 			     "CalibCoeffHcal" ,
135 			     _CalibrHCAL ,
136 			     defaultValuesForCalibrHCAL);
137
138 }
139
140
141
142 void CalibProcessor::init()
143 {
144
145 	// usually a good idea to
146 	printParameters();
147
148 	_nRun = -1 ;
149 	_nEvt = 0 ;
150
151
152 	// Number of bins in Historgrams are set here
153 	int Nbin_EHistogram2d_Ecal = _Nbins;
154 	int Nbin_EHistogram2d_Hcal = _Nbins;
155 	int Nbin_PartHistogram     = _Nbins / 4;
156
157 	// x-axis: Energy in ECAL: y-axis: Energy in HCAL
158 	_EHistogram2d_Ecal_Hcal = gsl_histogram2d_alloc ( Nbin_EHistogram2d_Ecal, Nbin_EHistogram2d_Hcal );
159 	// x-axis: Energy in HCAL-ECAL; y-axis: Energy in HCAL+ECAL
160 	_EHistogram2d_Rot = gsl_histogram2d_alloc ( Nbin_EHistogram2d_Ecal, Nbin_EHistogram2d_Hcal );
161
162
163 	// Histograms for the bins of the cloud
164 	_PartHistogram_1 = gsl_histogram_alloc ( Nbin_PartHistogram );
165 	_PartHistogram_2 = gsl_histogram_alloc ( Nbin_PartHistogram );
166 	_PartHistogram_3 = gsl_histogram_alloc ( Nbin_PartHistogram );
167 	_PartHistogram_4 = gsl_histogram_alloc ( Nbin_PartHistogram );
168 	_PartHistogram_5 = gsl_histogram_alloc ( Nbin_PartHistogram );
169 	_PartHistogram_6 = gsl_histogram_alloc ( Nbin_PartHistogram );
170 	_PartHistogram_7 = gsl_histogram_alloc ( Nbin_PartHistogram );
171 	_PartHistogram_8 = gsl_histogram_alloc ( Nbin_PartHistogram );
172
173 	_histoVector.push_back(_PartHistogram_1);
174 	_histoVector.push_back(_PartHistogram_2);
175 	_histoVector.push_back(_PartHistogram_3);
176 	_histoVector.push_back(_PartHistogram_4);
177 	_histoVector.push_back(_PartHistogram_5);
178 	_histoVector.push_back(_PartHistogram_6);
179 	_histoVector.push_back(_PartHistogram_7);
180 	_histoVector.push_back(_PartHistogram_8);
181
182 	for ( int sliceCounter=0; sliceCounter<8; sliceCounter++) {
183 	  nEntries[sliceCounter] = 0;
184 	}
185 }
186
187
188
190 {
191   _nRun++ ;
192 }
193
194
195
196 void CalibProcessor::processEvent( LCEvent * evt )
197 {
198   static bool firstEvent = true ;
199
200   if ( firstEvent==true ) {
201     std::cout << "CalibProcessor called for first event" << std::endl;
202
203   }
204
205
206 #ifdef MARLIN_USE_AIDA
207
208   try {
209
210     double entry_ecal=0;
211     double entry_hcal=0;
212
213     // create instances of AIDA CLOUDS
214     static AIDA::ICloud1D* Slice_1;
215     static AIDA::ICloud1D* Slice_2;
216     static AIDA::ICloud1D* Slice_3;
217     static AIDA::ICloud1D* Slice_4;
218     static AIDA::ICloud1D* Slice_5;
219     static AIDA::ICloud1D* Slice_6;
220     static AIDA::ICloud1D* Slice_7;
221     static AIDA::ICloud1D* Slice_8;
222     static AIDA::ICloud2D* Cloud;
223     static AIDA::ICloud2D* CloudRot;
224
225     if ( isFirstEvent() ) {
226
227       /* **** ************************ **** */
228       /* **** BEGIN HISTOGRAM SETTINGS **** */
229       /* **** ************************ **** */
230
231       /* -------------- */
232       /* GSL HISTOGRAMS */
233       /* -------------- */
234
235       double range_min_EHistogram2d_Ecal_Hcal = 0;
236       double range_max_EHistogram2d_Ecal_Hcal = _E_CMS + 100;
237       gsl_histogram2d_set_ranges_uniform ( _EHistogram2d_Ecal_Hcal,				                 // histogram
238 					   range_min_EHistogram2d_Ecal_Hcal, range_max_EHistogram2d_Ecal_Hcal,	 // x-axis ranges
239 					   range_min_EHistogram2d_Ecal_Hcal, range_max_EHistogram2d_Ecal_Hcal ); // y-axis ranges
240
241       double range_min_x_EHistogram2d_Rot = -(_E_CMS+100) / 2;
242       double range_max_x_EHistogram2d_Rot = +(_E_CMS+100) / 2;
243       double range_min_y_EHistogram2d_Rot = 0;
244       double range_max_y_EHistogram2d_Rot = _E_CMS+100;
245       gsl_histogram2d_set_ranges_uniform ( _EHistogram2d_Rot,			                                 // histogram
246 					   range_min_x_EHistogram2d_Rot, range_max_x_EHistogram2d_Rot,		 // x-axis ranges
247 					   range_min_y_EHistogram2d_Rot, range_max_y_EHistogram2d_Rot );	 // y-axis ranges
248
249       double range_min_PartHistogram = 0;
250       double range_max_PartHistogram = _E_CMS+100;
251       gsl_histogram_set_ranges_uniform ( _PartHistogram_1, range_min_PartHistogram, range_max_PartHistogram );
252       gsl_histogram_set_ranges_uniform ( _PartHistogram_2, range_min_PartHistogram, range_max_PartHistogram );
253       gsl_histogram_set_ranges_uniform ( _PartHistogram_3, range_min_PartHistogram, range_max_PartHistogram );
254       gsl_histogram_set_ranges_uniform ( _PartHistogram_4, range_min_PartHistogram, range_max_PartHistogram );
255       gsl_histogram_set_ranges_uniform ( _PartHistogram_5, range_min_PartHistogram, range_max_PartHistogram );
256       gsl_histogram_set_ranges_uniform ( _PartHistogram_6, range_min_PartHistogram, range_max_PartHistogram );
257       gsl_histogram_set_ranges_uniform ( _PartHistogram_7, range_min_PartHistogram, range_max_PartHistogram );
258       gsl_histogram_set_ranges_uniform ( _PartHistogram_8, range_min_PartHistogram, range_max_PartHistogram );
259
260       // create AIDA clouds
261       Slice_1 = AIDAProcessor::histogramFactory(this)->createCloud1D( "Slice_1", "Hits in Slice_1", -1 );
262       Slice_2 = AIDAProcessor::histogramFactory(this)->createCloud1D( "Slice_2", "Hits in Slice_2", -1 );
263       Slice_3 = AIDAProcessor::histogramFactory(this)->createCloud1D( "Slice_3", "Hits in Slice_3", -1 );
264       Slice_4 = AIDAProcessor::histogramFactory(this)->createCloud1D( "Slice_4", "Hits in Slice_4", -1 );
265       Slice_5 = AIDAProcessor::histogramFactory(this)->createCloud1D( "Slice_5", "Hits in Slice_5", -1 );
266       Slice_6 = AIDAProcessor::histogramFactory(this)->createCloud1D( "Slice_6", "Hits in Slice_6", -1 );
267       Slice_7 = AIDAProcessor::histogramFactory(this)->createCloud1D( "Slice_7", "Hits in Slice_7", -1 );
268       Slice_8 = AIDAProcessor::histogramFactory(this)->createCloud1D( "Slice_8", "Hits in Slice_8", -1 );
269
270       Cloud = AIDAProcessor::histogramFactory(this)->createCloud2D( "Cloud", "x: Eecal, y: Ehcal", -1 );
271       CloudRot = AIDAProcessor::histogramFactory(this)->createCloud2D( "CloudRot", "x: Ehcal-Eecal, y: Ehcal+Eecal", -1 );
272
273     }
274     /* **** ************************* **** */
275     /* **** END OF HISTOGRAM SETTINGS **** */
276     /* **** ************************* **** */
277
278
279     std::vector< std::string >::const_iterator iter;
280     const std::vector< std::string >* ColNames = evt->getCollectionNames();
281
282
283
284     /**** ************************* **** */
285     /**** loop over all Collections **** */
286     /**** ************************* **** */
287     for ( iter = ColNames->begin() ; iter != ColNames->end() ; iter++ ) {
288
289       LCCollection* col = evt->getCollection( *iter ) ;
290
291
292
293       /* *** ********************************* *** */
294       /* *** Process Hadronic Calorimeter Hits *** */
295       /* *** ********************************* *** */
296
297       if ( (col->getTypeName() == LCIO::CALORIMETERHIT) && (*iter == _colNameHCALHits) ) {
298
299 	int nHcalElements = col->getNumberOfElements();
300
301 	double sum_hcal_energy = 0;
302
303 	//looping over all hits
304 	for (int hitCounter = 0 ; hitCounter < nHcalElements ; ++hitCounter ) {
305
306 	  CalorimeterHit* CaloHit =  static_cast<CalorimeterHit*>(col->getElementAt(hitCounter));// CLUSTER
307 	  double energy = CaloHit->getEnergy();
308 	  sum_hcal_energy += energy; //*0.7 to fake miscalibration
309
310 	}
311 	entry_hcal=sum_hcal_energy;
312       }// end of processing hits in hcal
313
314
315
316       /* *** **************************************** *** */
317       /* *** Process Electromagnetic Calorimeter Hits *** */
318       /* *** **************************************** *** */
319
320       if ( (col->getTypeName() == LCIO::CALORIMETERHIT) && (*iter == _colNameECALHits) ) {
321
322 	int nEcalElements = col->getNumberOfElements();
323
324 	double sum_ecal_energy=0;
325
326 	//looping over all hits
327 	for (int hitCounter = 0; hitCounter < nEcalElements; ++hitCounter ) {
328
329 	  CalorimeterHit* CaloHit = static_cast<CalorimeterHit*>(col->getElementAt(hitCounter));
330 	  double energy = CaloHit->getEnergy();
331 	  sum_ecal_energy += energy;
332 	}
333 	entry_ecal = sum_ecal_energy;
334       }// end of processing hits in ECAL
335     }// end of for loop over collections
336
337     // fill 2d-histogram with hcal and ecal data
338     gsl_histogram2d_increment ( _EHistogram2d_Ecal_Hcal, entry_ecal, entry_hcal );
339
340     // coordinate transformation to 'rotate' the entries in x-y plane:
341     double x_rot = entry_hcal - entry_ecal;
342     double y_rot = entry_hcal + entry_ecal;
343
344     // fill 2d-histogram using new coordinates:
345     gsl_histogram2d_increment ( _EHistogram2d_Rot, x_rot, y_rot );
346     // fill 2d-clouds:
347     Cloud->fill(entry_ecal, entry_hcal);
348     CloudRot->fill(x_rot, y_rot);
349
350
351     // fill the entries of the data sample into 8 binnings :
352     // FIXME: NO HARDCODED NUMBERS! -> automated binning
353     // change in number of bins will cause further changes in the code :-)
354     if ( x_rot < -100 ) {
355       Slice_1->fill(y_rot);
356       gsl_histogram_increment ( _PartHistogram_1, y_rot );
357       nEntries[0]++;
358     }
359     else if ( (x_rot >= -100 ) && (x_rot < -50) ) {
360       Slice_2->fill(y_rot);
361       gsl_histogram_increment ( _PartHistogram_2, y_rot );
362       nEntries[1] ++;
363     }
364     else if ( (x_rot >= -50 ) && (x_rot < 0) ) {
365       Slice_3->fill(y_rot);
366       gsl_histogram_increment ( _PartHistogram_3, y_rot );
367       nEntries[2]++;
368     }
369     else if ( (x_rot >= 0 ) && (x_rot < 25) ) {
370       Slice_4->fill(y_rot);
371       gsl_histogram_increment ( _PartHistogram_4, y_rot );
372       nEntries[3]++;
373     }
374     else if ( (x_rot >= 25 ) && (x_rot < 50) ) {
375       Slice_5->fill(y_rot);
376       gsl_histogram_increment ( _PartHistogram_5, y_rot );
377       nEntries[4]++;
378     }
379     else if ( (x_rot >= 50 ) && (x_rot < 100) ) {
380       Slice_6->fill(y_rot);
381       gsl_histogram_increment ( _PartHistogram_6, y_rot );
382       nEntries[5]++;
383     }
384     else if ( (x_rot >= 100 ) && (x_rot < 150) ) {
385       Slice_7->fill(y_rot);
386       gsl_histogram_increment ( _PartHistogram_7, y_rot );
387       nEntries[6]++;
388     }
389     else if ( x_rot > 150 ) {
390       Slice_8->fill(y_rot);
391       gsl_histogram_increment ( _PartHistogram_8, y_rot );
392       nEntries[7]++;
393     }
394
395   }	// end of try
396
397
398   catch(DataNotAvailableException &e){
399     std::cout << "!! CalibProcessor: collection not available !!" << std::endl ;
400   };
401
402 #endif
403
404
405   _nEvt ++;
406   firstEvent = false ;
407
408 }	// end of process event block
409
410
411
412 void CalibProcessor::end()
413 {
414
415   /* *** *********************************************** *** */
416   /* *** Perform Gauss fit to find mean energy per Slice *** */
417   /* *** *********************************************** *** */
418
419   // local variables
420   int npar = 3; // 8 times 3 parameters to fit:
421   double mean[8]; double sigma[8]; double norm[8];
422
423
424   // converging criteria
425   const double abs_error = 1e-4;
426   const double rel_error = 1e-4;
427
428   //find starting value for fit: mean = bin with maximum entry / normalisation = number of entries
429   double max_bin_PartHistogram_lower_range_limit[8];
430   double max_bin_PartHistogram_upper_range_limit[8];
431   int    max_bin_PartHistogram[8];
432
433   for ( int sliceCounter = 0; sliceCounter<8; sliceCounter++ ) {
434
435     max_bin_PartHistogram[sliceCounter] = gsl_histogram_max_bin (_histoVector.at(sliceCounter));
436
437     gsl_histogram_get_range (_histoVector.at(sliceCounter), max_bin_PartHistogram[sliceCounter], &max_bin_PartHistogram_lower_range_limit[sliceCounter], &max_bin_PartHistogram_upper_range_limit[sliceCounter]);
438
439     mean[sliceCounter] = (max_bin_PartHistogram_lower_range_limit[sliceCounter]
440 			  + (max_bin_PartHistogram_upper_range_limit[sliceCounter]))/2;
441
442
443     int nBins = gsl_histogram_bins (_histoVector.at(sliceCounter));
444     double xValue[nBins];
445     double yValue[nBins];
446
447     for(int binCounter = 0; binCounter < nBins; binCounter++) {
448
449       double lowerEdge = 0.0;
450       double upperEdge = 0.0;
451       // int status = gsl_histogram_get_range(_histoVector.at(sliceCounter),binCounter,&lowerEdge,&upperEdge);
452
453       xValue[binCounter]= ( upperEdge - lowerEdge ) / 2.0;
454       yValue[binCounter]=gsl_histogram_get(_histoVector.at(sliceCounter),binCounter);
455
456     }
457
458     gsl_multifit_function_fdf fitfunc;
459
460     const gsl_multifit_fdfsolver_type* T = gsl_multifit_fdfsolver_lmsder;
461     gsl_multifit_fdfsolver* Solver = gsl_multifit_fdfsolver_alloc(T,nBins,npar);
462     gsl_matrix* covar = gsl_matrix_alloc(npar,npar);   // covariance matrix
463
464     int status   = 0;
465     int iter     = 0;
466     int max_iter = 20;
467
468     //initialize fit starting values
469     double par_init[npar];
470     par_init[0] = mean[sliceCounter];
471     par_init[1] = 20;
472     par_init[2] = nEntries[sliceCounter];
473
474     data DataSet;
475     DataSet.n = nBins;
476     DataSet.x = &xValue[0];
477     DataSet.y = &yValue[0];
478
479     fitfunc.f      = &GausFunc;
480     fitfunc.df     = &dGausFunc;
481     fitfunc.fdf    = &fdfGausFunc;
482     fitfunc.n      = nBins;
483     fitfunc.p      = npar;
484     fitfunc.params = &DataSet;
485
486     gsl_vector_view pinit = gsl_vector_view_array(par_init,npar);
487     gsl_multifit_fdfsolver_set(Solver,&fitfunc,&pinit.vector);
488
489
490     // perform fit
491     do {
492       iter++;
493       std::cout << "Gauss Fit Iteration started ... ... "<<std::endl;
494       status = gsl_multifit_fdfsolver_iterate(Solver);
495
496       // DEBUG: Fit status
497       std::cout << "iteration: " << iter
498 		<<" slice: " << sliceCounter
499 		<<" mean: " << gsl_vector_get(Solver->x,0)
500 		<<" sigma: " << gsl_vector_get(Solver->x,1)
501 		<<" normalisation: " << gsl_vector_get(Solver->x,2)<<std::endl;
502       std::cout << "--- DONE ---" << std::endl;
503
504       if (status) break;
505       status = gsl_multifit_test_delta (Solver->dx, Solver->x,abs_error,rel_error);
506
507     } while ( status==GSL_CONTINUE && iter < max_iter);
508
509     gsl_multifit_covar (Solver->J, rel_error, covar);
510
511
512     mean[sliceCounter] = gsl_vector_get(Solver->x,0); // mean
513     sigma[sliceCounter] = gsl_vector_get(Solver->x,1); // sigm
514     norm[sliceCounter] = gsl_vector_get(Solver->x,2); // normalisation
515
516
517     gsl_multifit_fdfsolver_free(Solver);
518     gsl_matrix_free(covar);
519   }
520
521   //bin centers of (H-E) bins of the "rotated cloud" according to cuts set in the various if statments above
522   // FIXME: NO HARD CODED NUMBERS!
523   double x[8] = {-200, -75, -25, 12.5, 37.5, 75, 125, 225}; // bin center [GeV]
524   double xErr[8] = {100, 25, 25, 12.5, 12.5, 25, 25, 75}; // bin width [GeV]
525
526   // fit y(x) = b + a*x for the dataset (x, y)
527   int xstride = 1;
528   int ystride = 1;
529   int nSlices = 8;
530   double a, b, cov00, cov01, cov11, sumsq;
531   gsl_fit_linear (x, xstride, mean, ystride, nSlices, &b, &a, &cov00, &cov01, &cov11, &sumsq);
532
533
534   // calculate a_prime, b_prime in (x=Eecal,y=Ehacal) coordinates
535   // the function describing the edge of the contour plot is f(Eecal)= b_prime + a_prime * Eecal
536   double a_prime, b_prime;
537   a_prime = - (1.0+a) / (1.0-a); // slope
538   b_prime = b / (1.0-a);			// y-offset
539
540
541   //double sampling_fraction_coeff[3] = {33.0235, 93.5682, 21.19626};
542   double sampling_fraction_coeff_calib[3];
543   double f = _E_CMS/b_prime;
544
545   // a_prime is considered to be always negative. use (-1)*a_prime to get abs(a_prime)
546   sampling_fraction_coeff_calib[0] = f*(-1)*a_prime*_CalibrECAL.at(0);
547   sampling_fraction_coeff_calib[1] = f*(-1)*a_prime*_CalibrECAL.at(1);
548   sampling_fraction_coeff_calib[2] = f*_CalibrHCAL.at(0);
549
550
551
552   /* *********************** */
553   /*     PROCESSOR OUTPUT    */
554   /* *********************** */
555
556   /* --------------- */
557   /* PRINT TO SCREEN */
558   /* --------------- */
559   std::cout << std::endl;
560   std::cout << std::endl;
561   std::cout << "******************************" << std::endl;
562   std::cout << "*********  RESULTS   *********" << std::endl;
563   std::cout << "******************************" << std::endl;
564   std::cout << "**  Parameters determined by gsl_fit_linear:" << std::endl;
565   std::cout << "**  a = "<< a << std::endl;
566   std::cout << "**  b = "<< b << std::endl;
567   std::cout << "**  Parameters describing edge. f(x)=a*x+b:" << std::endl;
568   std::cout << "**  a = "<< a_prime << std::endl;
569   std::cout << "**  b = "<< b_prime << std::endl;
570
571
572   std::cout << "**  Old sampling fraction coefficients:" << std::endl;
573   for (int i=0 ; i<2 ; i++) {
574     std::cout<< "**  c" << i+1 << " = " << _CalibrECAL.at(i) << std::endl;
575   }
576   std::cout << "**  c" << 3 << " = " << _CalibrHCAL.at(0) << std::endl;
577
578
579   std::cout << "**  New sampling fraction coefficients:" << std::endl;
580   for (int i=0; i<3; i++) {
581     std::cout<< "**  c" << i+1 << "_calib = " << sampling_fraction_coeff_calib[i] << std::endl;
582   }
583   std::cout << "**  " <<  std::endl;
584   std::cout << "**  f = "<< f << std::endl;
585   std::cout << "**  (-1)*f*a_prime = "<< (-1)*f*a_prime << std::endl;
586   std::cout << "******************************" << std::endl;
587
588
589   /* --------------- */
590   /* PRINT TO FILE   */
591   /* --------------- */
592
593   //creating textfile 'Calibration-coefficients.txt' to write output into
594   std::ofstream myfile;
595   myfile.open ("CalibProcessor_results.txt");
596
597   myfile << "Energy per Slice [GeV] \n";
598   for (int i=0; i<8; i++) {
599     myfile << "slice" << i+1
600 	   << " meanX " << x[i] << " errX "<< xErr[i]
601 	   << " meanY " << mean[i] << " sigma " << sigma[i] << "\n";
602   }
603
604   myfile << "Parameters describing edge. f(x)=a*x+b: \n";
605   myfile << "a " << a << " b " << b << " \n";
606
607   myfile << "New calibration coefficients: \n";
608   for (int i=0; i<3; i++) {
609     myfile << "c" << i+1 << "_calib " << sampling_fraction_coeff_calib[i] << "\n";
610   }
611   myfile.close();
612
613   /* **************** */
614   /* FREEE DISK SPACE */
615   /* **************** */
616   gsl_histogram2d_free ( _EHistogram2d_Ecal_Hcal);
617   gsl_histogram2d_free ( _EHistogram2d_Rot);
618   gsl_histogram_free (_PartHistogram_1 );
619   gsl_histogram_free (_PartHistogram_2 );
620   gsl_histogram_free (_PartHistogram_3 );
621   gsl_histogram_free (_PartHistogram_4 );
622   gsl_histogram_free (_PartHistogram_5 );
623   gsl_histogram_free (_PartHistogram_6 );
624   gsl_histogram_free (_PartHistogram_7 );
625   gsl_histogram_free (_PartHistogram_8 );
626
627
628   // Improve the code. Thank you.
629   // Questions: viktor@access.uzh.ch
630   // Sept. 2007
631 }
```

## Attached Files

To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.
• [get | view] (2008-05-13 08:17:42, 20.9 KB) [[attachment:CalibProcessor.cc]]
• [get | view] (2008-04-29 11:06:29, 26.1 KB) [[attachment:CalibProcessor_data_histo_contour.eps]]
• [get | view] (2008-04-28 13:41:33, 441.3 KB) [[attachment:documentation]]
• [get | view] (2008-04-28 13:39:36, 589.8 KB) [[attachment:tarball]]
All files | Selected Files: delete move to page copy to page