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 
 189 void CalibProcessor::processRunHeader( LCRunHeader* run)
 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 }