EcalTPGParamBuilder.cc

Go to the documentation of this file.

#include "EcalTPGParamBuilder.h"

#include "CalibCalorimetry/EcalTPGTools/plugins/EcalTPGDBApp.h"

#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
#include "Geometry/Records/interface/CaloGeometryRecord.h"
#include "Geometry/EcalMapping/interface/EcalElectronicsMapping.h"
#include "Geometry/EcalMapping/interface/EcalMappingRcd.h"

#include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
#include "DataFormats/EcalDetId/interface/EBDetId.h"
#include "DataFormats/EcalDetId/interface/EEDetId.h"

#include "CondFormats/DataRecord/interface/EcalIntercalibConstantsRcd.h"
#include "CondFormats/EcalObjects/interface/EcalADCToGeVConstant.h"
#include "CondFormats/DataRecord/interface/EcalADCToGeVConstantRcd.h"
#include "CondFormats/EcalObjects/interface/EcalMGPAGainRatio.h"
#include "CondFormats/DataRecord/interface/EcalGainRatiosRcd.h"
#include "CondFormats/DataRecord/interface/EcalPedestalsRcd.h"
#include "CondFormats/EcalObjects/interface/EcalTPGPedestals.h"
#include "CondFormats/DataRecord/interface/EcalTPGPedestalsRcd.h"

//modif-alex-27-july-2015
#include "CalibCalorimetry/EcalLaserCorrection/interface/EcalLaserDbService.h"
#include "CalibCalorimetry/EcalLaserCorrection/interface/EcalLaserDbRecord.h"

#include "SimCalorimetry/EcalSimAlgos/interface/EcalSimParameterMap.h"
#include "SimCalorimetry/EcalSimAlgos/interface/EBShape.h"
#include "SimCalorimetry/EcalSimAlgos/interface/EEShape.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include <iostream>
#include <string>
#include <sstream>
#include <vector>
#include <ctime>

#include <TF1.h>
#include <TH2F.h>
#include <TFile.h>
#include <TNtuple.h>
#include <iomanip>
#include <fstream>


using namespace std;

double oneOverEtResolEt(double *x, double *par) { 
  double Et = x[0] ;
  if (Et<1e-6) return 1./par[1] ; // to avoid division by 0.
  double resolEt_overEt = sqrt( (par[0]/sqrt(Et))*(par[0]/sqrt(Et)) + (par[1]/Et)*(par[1]/Et) + par[2]*par[2] ) ;
  return 1./(Et*resolEt_overEt) ;
}

EcalTPGParamBuilder::EcalTPGParamBuilder(edm::ParameterSet const& pSet)
  : xtal_LSB_EB_(0), xtal_LSB_EE_(0), nSample_(5), complement2_(7), useDBShape_(true)
{
  ped_conf_id_=0;
  lin_conf_id_=0;
  lut_conf_id_=0;
  wei_conf_id_=0;
  fgr_conf_id_=0;
  sli_conf_id_=0;
  spi_conf_id_=0; //modif-alex 21/01/11
  del_conf_id_=0; //modif-alex 21/01/11
  bxt_conf_id_=0;
  btt_conf_id_=0;
  bst_conf_id_=0;
  tag_="";
  version_=0;

  m_write_ped=1;
  m_write_lin=1;
  m_write_lut=1;
  m_write_wei=1;
  m_write_fgr=1;
  m_write_sli=1;
  m_write_spi=1; //modif-alex 21/01/11
  m_write_del=1; //modif-alex 21/01/11
  m_write_bxt=1;
  m_write_btt=1;
  m_write_bst=1;

  writeToDB_  = pSet.getParameter<bool>("writeToDB") ;
  DBEE_ = pSet.getParameter<bool>("allowDBEE") ;
  string DBsid    = pSet.getParameter<std::string>("DBsid") ;
  string DBuser   = pSet.getParameter<std::string>("DBuser") ;
  string DBpass   = pSet.getParameter<std::string>("DBpass") ;
  //uint32_t DBport = pSet.getParameter<unsigned int>("DBport") ;
  
  tag_   = pSet.getParameter<std::string>("TPGtag") ;
  version_ = pSet.getParameter<unsigned int>("TPGversion") ;

  m_write_ped = pSet.getParameter<unsigned int>("TPGWritePed") ;
  m_write_lin = pSet.getParameter<unsigned int>("TPGWriteLin") ;
  m_write_lut = pSet.getParameter<unsigned int>("TPGWriteLut") ;
  m_write_wei = pSet.getParameter<unsigned int>("TPGWriteWei") ;
  m_write_fgr = pSet.getParameter<unsigned int>("TPGWriteFgr") ;
  m_write_sli = pSet.getParameter<unsigned int>("TPGWriteSli") ;
  m_write_spi = pSet.getParameter<unsigned int>("TPGWriteSpi") ; //modif-alex 21/01/11
  m_write_del = pSet.getParameter<unsigned int>("TPGWriteDel") ; //modif-alex 21/01/11
  m_write_bxt = pSet.getParameter<unsigned int>("TPGWriteBxt") ;
  m_write_btt = pSet.getParameter<unsigned int>("TPGWriteBtt") ;
  m_write_bst = pSet.getParameter<unsigned int>("TPGWriteBst") ;

  btt_conf_id_=m_write_btt;
  bxt_conf_id_=m_write_bxt;
  bst_conf_id_=m_write_bst;

  if(m_write_ped != 0 && m_write_ped != 1 ) ped_conf_id_=m_write_ped;
  
    if (writeToDB_) edm::LogInfo("TopInfo") << "data will be saved with tag and version=" << tag_ << ".version" << version_<< "\n";
    db_ = new EcalTPGDBApp(DBsid, DBuser, DBpass) ;

  writeToFiles_ =  pSet.getParameter<bool>("writeToFiles") ;
  if (writeToFiles_) {
    std::string outFile = pSet.getParameter<std::string>("outFile") ;
    out_file_ = new std::ofstream(outFile.c_str(), std::ios::out) ;  
  }
  geomFile_   = new std::ofstream("geomFile.txt", std::ios::out) ;  


  useTransverseEnergy_ = pSet.getParameter<bool>("useTransverseEnergy") ;
  
  Et_sat_EB_ = pSet.getParameter<double>("Et_sat_EB") ;
  Et_sat_EE_ = pSet.getParameter<double>("Et_sat_EE") ;
  sliding_ = pSet.getParameter<unsigned int>("sliding") ;
  weight_timeShift_ = pSet.getParameter<double>("weight_timeShift") ;
  sampleMax_ = pSet.getParameter<unsigned int>("weight_sampleMax") ;
  weight_unbias_recovery_ = pSet.getParameter<bool>("weight_unbias_recovery") ;

  forcedPedestalValue_ = pSet.getParameter<int>("forcedPedestalValue") ;
  forceEtaSlice_ = pSet.getParameter<bool>("forceEtaSlice") ;
    
  LUT_option_ = pSet.getParameter<std::string>("LUT_option") ;
  LUT_threshold_EB_ = pSet.getParameter<double>("LUT_threshold_EB") ;
  LUT_threshold_EE_ = pSet.getParameter<double>("LUT_threshold_EE") ;
  LUT_stochastic_EB_ = pSet.getParameter<double>("LUT_stochastic_EB") ;
  LUT_noise_EB_ =pSet.getParameter<double>("LUT_noise_EB") ;
  LUT_constant_EB_ =pSet.getParameter<double>("LUT_constant_EB") ;
  LUT_stochastic_EE_ = pSet.getParameter<double>("LUT_stochastic_EE") ;
  LUT_noise_EE_ =pSet.getParameter<double>("LUT_noise_EE") ;
  LUT_constant_EE_ =pSet.getParameter<double>("LUT_constant_EE") ;

  TTF_lowThreshold_EB_ = pSet.getParameter<double>("TTF_lowThreshold_EB") ;
  TTF_highThreshold_EB_ = pSet.getParameter<double>("TTF_highThreshold_EB") ;
  TTF_lowThreshold_EE_ = pSet.getParameter<double>("TTF_lowThreshold_EE") ;
  TTF_highThreshold_EE_ = pSet.getParameter<double>("TTF_highThreshold_EE") ;

  FG_lowThreshold_EB_ = pSet.getParameter<double>("FG_lowThreshold_EB") ;
  FG_highThreshold_EB_ = pSet.getParameter<double>("FG_highThreshold_EB") ;
  FG_lowRatio_EB_ = pSet.getParameter<double>("FG_lowRatio_EB") ;
  FG_highRatio_EB_ = pSet.getParameter<double>("FG_highRatio_EB") ;
  FG_lut_EB_ = pSet.getParameter<unsigned int>("FG_lut_EB") ;
  FG_Threshold_EE_ = pSet.getParameter<double>("FG_Threshold_EE") ;
  FG_lut_strip_EE_ = pSet.getParameter<unsigned int>("FG_lut_strip_EE") ;
  FG_lut_tower_EE_ = pSet.getParameter<unsigned int>("FG_lut_tower_EE") ;
  SFGVB_Threshold_ = pSet.getParameter<unsigned int>("SFGVB_Threshold") ;
  SFGVB_lut_ = pSet.getParameter<unsigned int>("SFGVB_lut") ;
  SFGVB_SpikeKillingThreshold_ = pSet.getParameter<int>("SFGVB_SpikeKillingThreshold") ; //modif-alex 21/01/11
  pedestal_offset_ = pSet.getParameter<unsigned int>("pedestal_offset") ;

  useInterCalibration_  = pSet.getParameter<bool>("useInterCalibration") ;
  H2_ = pSet.getUntrackedParameter<bool>("H2",false) ;

  useTransparencyCorr_ = false;
  useTransparencyCorr_ = pSet.getParameter<bool>("useTransparencyCorr") ; //modif-alex-25/04/2012
  Transparency_Corr_   = pSet.getParameter<std::string>("transparency_corrections");//modif-alex-30/01/2012 

  //modif-alex-23/02/2011
  //convert the spike killing first from GeV to ADC (10 bits)                                                                                                                                               
  //depending on the saturation scale: Et_sat_EB_                                                                                                                                                                
  if(SFGVB_SpikeKillingThreshold_ == -1 || (SFGVB_SpikeKillingThreshold_ > Et_sat_EB_)) 
    SFGVB_SpikeKillingThreshold_ = 1023; //nokilling
  else 
    SFGVB_SpikeKillingThreshold_ = int(SFGVB_SpikeKillingThreshold_ * 1024/Et_sat_EB_);    
  edm::LogInfo("TopInfo") << "INFO:SPIKE KILLING THRESHOLD (ADC)=" << SFGVB_SpikeKillingThreshold_ << "\n";
  
  //modif-alex-02/02/11
  //TIMING information 
  TimingDelays_EB_ = pSet.getParameter<std::string>("timing_delays_EB") ;
  TimingDelays_EE_ = pSet.getParameter<std::string>("timing_delays_EE") ;
  TimingPhases_EB_ = pSet.getParameter<std::string>("timing_phases_EB") ;
  TimingPhases_EE_ = pSet.getParameter<std::string>("timing_phases_EE") ;

  std::ostringstream ss; 
  //  edm::LogInfo("TopInfo") << "INFO: READING timing files" << "\n";
  ss << "INFO: READING timing files\n";
  std::ifstream delay_eb(TimingDelays_EB_.c_str()); 
  if(!delay_eb) edm::LogError("TopInfo") << "ERROR: File " << TimingDelays_EB_.c_str() << " could not be opened" << "\n";
  std::ifstream delay_ee(TimingDelays_EE_.c_str());
  if(!delay_ee) edm::LogError("TopInfo") << "ERROR: File " << TimingDelays_EE_.c_str() << " could not be opened" << "\n";
  std::ifstream phase_eb(TimingPhases_EB_.c_str());
  if(!phase_eb) edm::LogError("TopInfo") << "ERROR: File " << TimingPhases_EB_.c_str() << " could not be opened" << "\n";
  std::ifstream phase_ee(TimingPhases_EE_.c_str());
  if(!phase_ee) edm::LogError("TopInfo") << "ERROR: File " << TimingPhases_EE_.c_str() << " could not be opened" << "\n";
  
  char buf[1024];
  //READING DELAYS EB
  delay_eb.getline(buf,sizeof(buf),'\n');
  while( delay_eb ) {
    std::stringstream sin(buf);

    int tcc; 
    sin >> tcc;

    vector<int> vec_delays_eb;
    for(int ieb=0; ieb<68; ++ieb){
      int time_delay = -1;
      sin >> time_delay;
      vec_delays_eb.push_back(time_delay);
      if(time_delay==-1)  edm::LogError("TopInfo") << "ERROR:Barrel timing delay -1, check file" << "\n";
    }

    if(vec_delays_eb.size()!=68)
      edm::LogError("TopInfo") << "ERROR:Barrel timing delay wrong, not enough towers, check file" << "\n";

    if (delays_EB_.find(tcc) == delays_EB_.end())
      delays_EB_.insert(make_pair(tcc,vec_delays_eb));

    //    edm::LogInfo("TopInfo") << tcc << "\n"; 
    ss << tcc;
    for(unsigned int ieb=0; ieb<vec_delays_eb.size(); ++ieb)
      //      edm::LogInfo("TopInfo") << vec_delays_eb[ieb] << "\n";
      ss << " " << vec_delays_eb[ieb];
    delay_eb.getline(buf,sizeof(buf),'\n');
    ss << "\n";
  }//loop delay file EB
  delay_eb.close();

  //READING PHASES EB
  phase_eb.getline(buf,sizeof(buf),'\n');
  while( phase_eb ) {
    std::stringstream sin(buf);
    int tcc;
    sin >> tcc;

    vector<int> vec_phases_eb;
    for(unsigned int ieb=0; ieb<68; ++ieb){
      int time_phase = -1;
      sin >> time_phase;
      vec_phases_eb.push_back(time_phase);
      if(time_phase==-1)  edm::LogError("TopInfo") << "ERROR:Barrel timing phase -1, check file" << "\n";
    }

    if(vec_phases_eb.size()!=68)
      edm::LogError("TopInfo") << "ERROR:Barrel timing phase wrong, not enough towers, check file" << "\n";

    if (phases_EB_.find(tcc) == phases_EB_.end())
      phases_EB_.insert(make_pair(tcc,vec_phases_eb));

    //    edm::LogInfo("TopInfo") << tcc << "\n"; 
    ss << tcc;
    for(unsigned int ieb=0; ieb<vec_phases_eb.size(); ++ieb)
      //      edm::LogInfo("TopInfo") << vec_phases_eb[ieb] << "\n";
      ss << " " << vec_phases_eb[ieb];
    phase_eb.getline(buf,sizeof(buf),'\n');
    ss << "\n";
  }//loop phase file EB
  phase_eb.close();

  //READING DELAYS EE//------------------------------------------------
  delay_ee.getline(buf,sizeof(buf),'\n');
  while( delay_ee ) {
    std::stringstream sin(buf);
    int tcc; 
    sin >> tcc;

    vector<int> vec_delays_ee;
    for(unsigned int iee=0; iee<48; ++iee){
      int time_delay = -1;
      sin >> time_delay;
      vec_delays_ee.push_back(time_delay);
      if(time_delay==-1)  edm::LogError("TopInfo") << "ERROR:EE timing delay -1, check file" << "\n";
    }

    if(vec_delays_ee.size()!=48)
      edm::LogError("TopInfo") << "ERROR:EE timing delay wrong, not enough towers, check file" << "\n";

    if (delays_EE_.find(tcc) == delays_EE_.end())
      delays_EE_.insert(make_pair(tcc,vec_delays_ee));

    //    edm::LogInfo("TopInfo") << tcc << "\n"; 
    ss << tcc;
    for(unsigned int iee=0; iee<vec_delays_ee.size(); ++iee)
      //      edm::LogInfo("TopInfo") << vec_delays_ee[iee] << "\n";
      ss << " " << vec_delays_ee[iee];
    ss << "\n";
    delay_ee.getline(buf,sizeof(buf),'\n');
  }//loop delay file EE
  delay_ee.close();

  //READING PHASES EE
  phase_ee.getline(buf,sizeof(buf),'\n');
  while( phase_ee ) {
    std::stringstream sin(buf);
    int tcc;
    sin >> tcc;

    vector<int> vec_phases_ee;
    for(unsigned int iee=0; iee<48; ++iee){
      int time_phase = -1;
      sin >> time_phase;
      vec_phases_ee.push_back(time_phase);
      if(time_phase==-1)  edm::LogError("TopInfo") << "ERROR:EE timing phase -1, check file" << "\n";
    }

    if(vec_phases_ee.size()!=48)
      edm::LogError("TopInfo") << "ERROR:EE timing phase wrong, not enough towers, check file" << "\n";

    if (phases_EE_.find(tcc) == phases_EE_.end())
      phases_EE_.insert(make_pair(tcc,vec_phases_ee));
    //    edm::LogInfo("TopInfo") << tcc << "\n";
    ss << tcc;
    for(unsigned int iee = 0; iee < vec_phases_ee.size(); ++iee)
      //      edm::LogInfo("TopInfo") << vec_phases_ee[iee] << "\n";
      ss << " " << vec_phases_ee[iee];
    ss << "\n";
    phase_ee.getline(buf,sizeof(buf),'\n');
  }//loop phase file EE
  phase_ee.close();

  //  edm::LogInfo("TopInfo") << "INFO: DONE reading timing files for EB and EE" << "\n";
  ss << "INFO: DONE reading timing files for EB and EE\n";
  edm::LogInfo("TopInfo") << ss.str();
}

EcalTPGParamBuilder::~EcalTPGParamBuilder()
{ 
  if (writeToFiles_) {
    (*out_file_ )<<"EOF"<<std::endl ;
    out_file_->close() ;
    delete out_file_ ;
  }
}


bool EcalTPGParamBuilder::checkIfOK(EcalPedestals::Item item) 
{
  bool result=true;
  if( item.mean_x1 <150. || item.mean_x1 >250) result=false;
  if( item.mean_x6 <150. || item.mean_x6 >250) result=false;
  if( item.mean_x12 <150. || item.mean_x12 >250) result=false;
  if( item.rms_x1 <0 || item.rms_x1 > 2) result=false;
  if( item.rms_x6 <0 || item.rms_x6 > 3) result=false;
  if( item.rms_x12 <0 || item.rms_x12 > 5) result=false;
  return result; 
}

int EcalTPGParamBuilder::getEtaSlice(int tccId, int towerInTCC)
{
  int etaSlice = (towerInTCC-1)/4+1 ;
  // barrel
  if (tccId>36 || tccId<73) return etaSlice ;
  //endcap
  else {
    if (tccId >=1 && tccId <= 18) etaSlice += 21 ; // inner -
    if (tccId >=19 && tccId <= 36) etaSlice += 17 ; // outer -
    if (tccId >=91 && tccId <= 108) etaSlice += 21 ; // inner +
    if (tccId >=73 && tccId <= 90) etaSlice += 17 ; // outer +
  }
  return etaSlice ;
}

void EcalTPGParamBuilder::analyze(const edm::Event& evt, const edm::EventSetup& evtSetup) 
{
  using namespace edm;
  using namespace std;

  // geometry
  ESHandle<CaloGeometry> theGeometry;
  ESHandle<CaloSubdetectorGeometry> theEndcapGeometry_handle, theBarrelGeometry_handle;
  evtSetup.get<CaloGeometryRecord>().get( theGeometry );
  evtSetup.get<EcalEndcapGeometryRecord>().get("EcalEndcap",theEndcapGeometry_handle);
  evtSetup.get<EcalBarrelGeometryRecord>().get("EcalBarrel",theBarrelGeometry_handle);
  evtSetup.get<IdealGeometryRecord>().get(eTTmap_);
  theEndcapGeometry_ = theEndcapGeometry_handle.product();
  theBarrelGeometry_ = theBarrelGeometry_handle.product();

  // electronics mapping
  ESHandle< EcalElectronicsMapping > ecalmapping;
  evtSetup.get< EcalMappingRcd >().get(ecalmapping);
  theMapping_ = ecalmapping.product();

  // get record for alpha
  std::ostringstream ss; 
  ss << "EcalLaserDbAnalyzer::analyze\n";
  edm::ESHandle<EcalLaserAlphas> handle;
  evtSetup.get<EcalLaserAlphasRcd>().get(handle);
  ss << "EcalLaserDbAnalyzer::analyze-> got EcalLaserDbRecord: \n";
  const EcalLaserAlphaMap& laserAlphaMap = handle.product()->getMap(); // map of apdpns

  //modif-alex-27-july-2015-+ Jean june 2016 beg
  // use alpha to check
  EcalLaserAlphaMap::const_iterator italpha;
  int cnt = 0;
  // Barrel loop
  for(italpha = laserAlphaMap.barrelItems().begin(); italpha != laserAlphaMap.barrelItems().end(); ++italpha) {
    if (cnt %1000 == 0) {
      EBDetId ebdetid = EBDetId::unhashIndex(cnt);
      ss << " Barrel ALPHA = " << (*italpha) << " cmsswId " << ebdetid.rawId() << "\n";
    }
    cnt++;
  } 
  ss << "Number of barrel Alpha parameters : " << cnt << "\n";
  // Endcap loop
  cnt = 0;
  for(italpha = laserAlphaMap.endcapItems().begin(); italpha != laserAlphaMap.endcapItems().end(); ++italpha) {
    if (cnt %1000 == 0) {
      EEDetId eedetid = EEDetId::unhashIndex(cnt);
      ss  << "EndCap ALPHA = " << (*italpha) << " cmsswId " << eedetid.rawId() << "\n";
    }
    cnt++;
  }
  ss << "Number of Endcap Alpha parameters : " << cnt << "\n";
  edm::LogInfo("TopInfo") << ss.str();
  ss.str("");
  //modif-alex-27-july-2015 +-Jean june 2016-end
  
  //modif-alex-30/01/2012   displaced in analyze Jean 2018
  ss.str("");
  if(useTransparencyCorr_){
    if(Transparency_Corr_ != "tag") {       // Jean 2018
      edm::LogInfo("TopInfo") << "INFO: READING transparency correction files" << "\n";
      ss << "INFO: READING transparency correction files\n";
      std::ifstream transparency(Transparency_Corr_.c_str());
      if(!transparency) edm::LogError("TopInfo") << "ERROR: File " << Transparency_Corr_.c_str() << " could not be opened" << "\n";
    
      char buf[1024];
      transparency.getline(buf,sizeof(buf),'\n');
      int xtalcounter = 0;
      while( transparency ) {
    std::stringstream sin(buf);
      
    int raw_xtal_id;
    sin >> raw_xtal_id;
      
    double xtal_trans_corr;
    sin >> xtal_trans_corr; 
      
    //      edm::LogInfo("TopInfo") << raw_xtal_id << " " << xtal_trans_corr << "\n";
    ss << raw_xtal_id << " " << xtal_trans_corr << "\n";
      
    Transparency_Correction_.insert(make_pair(raw_xtal_id,xtal_trans_corr));
      
    xtalcounter++;
    transparency.getline(buf,sizeof(buf),'\n');
      }//loop transparency
      transparency.close();
      ss << "INFO: DONE transparency correction files " << xtalcounter << "\n";
      edm::LogInfo("TopInfo") << ss.str();
      ss.str("");
    }   // file
    else {       // Jean 2018
      edm::LogInfo("TopInfo") << "INFO: READING transparency correction tag" << "\n";
      //      std::cout << "new feature, read a tag" << std::endl;
      ESHandle<EcalLaserAPDPNRatios> pAPDPNRatios;          
      evtSetup.get<EcalLaserAPDPNRatiosRcd>().get(pAPDPNRatios);  
      const EcalLaserAPDPNRatios* lratio = pAPDPNRatios.product();
      //      std::cout << " laser map size " << lratio->getLaserMap().size() << std::endl;

      EcalLaserAPDPNRatios::EcalLaserAPDPNRatiosMap::const_iterator itratio; 
      // Barrel loop
      for(int ib = 0; ib < 61200; ib++) {
    EBDetId ebdetid = EBDetId::unhashIndex(ib);
    itratio = lratio->getLaserMap().find(ebdetid.rawId());
    //  std::cout  << ib << " " << ebdetid.rawId() << " Barrel APDPNRatios = " << (*itratio).p1 << " " << (*itratio).p2 << " " << (*itratio).p3 << std::endl;
    Transparency_Correction_.insert(make_pair(ebdetid.rawId(), (*itratio).p2));
      }
      // Endcap loop
      for(int ie = 0; ie < 14648; ie++) {
    EEDetId eedetid = EEDetId::unhashIndex(ie);
    itratio = lratio->getLaserMap().find(eedetid.rawId());
    //  std::cout  << ie << " " << eedetid.rawId() << " Endcap APDPNRatios = " << (*itratio).p1 << " " << (*itratio).p2 << " " << (*itratio).p3 << std::endl;
    Transparency_Correction_.insert(make_pair(eedetid.rawId(), (*itratio).p2));
      }
      ss << "INFO: DONE transparency correction tag\n";
      edm::LogInfo("TopInfo") << ss.str();
      ss.str("");
    }       // Jean 2018
  }//if transparency

  // histo
  TFile saving ("EcalTPGParam.root","recreate") ;
  saving.cd () ;
  TH2F * ICEB = new TH2F("ICEB", "IC: Barrel", 360, 1, 361, 172, -86, 86) ;
  ICEB->GetYaxis()->SetTitle("eta index") ;
  ICEB->GetXaxis()->SetTitle("phi index") ;  
  TH2F * tpgFactorEB = new TH2F("tpgFactorEB", "tpgFactor: Barrel", 360, 1, 361, 172, -86, 86) ;
  tpgFactorEB->GetYaxis()->SetTitle("eta index") ;
  tpgFactorEB->GetXaxis()->SetTitle("phi index") ;  

  TH2F * ICEEPlus = new TH2F("ICEEPlus", "IC: Plus Endcap", 120, -9, 111, 120, -9, 111) ;
  ICEEPlus->GetYaxis()->SetTitle("y index") ;
  ICEEPlus->GetXaxis()->SetTitle("x index") ;
  TH2F * tpgFactorEEPlus = new TH2F("tpgFactorEEPlus", "tpgFactor: Plus Endcap", 120, -9, 111, 120, -9, 111) ;
  tpgFactorEEPlus->GetYaxis()->SetTitle("y index") ;
  tpgFactorEEPlus->GetXaxis()->SetTitle("x index") ;
  TH2F * ICEEMinus = new TH2F("ICEEMinus", "IC: Minus Endcap", 120, -9, 111, 120, -9, 111) ;
  ICEEMinus->GetYaxis()->SetTitle("y index") ;
  ICEEMinus->GetXaxis()->SetTitle("x index") ;
  TH2F * tpgFactorEEMinus = new TH2F("tpgFactorEEMinus", "tpgFactor: Minus Endcap", 120, -9, 111, 120, -9, 111) ;
  tpgFactorEEMinus->GetYaxis()->SetTitle("y index") ;
  tpgFactorEEMinus->GetXaxis()->SetTitle("x index") ;

  TH2F * IC = new TH2F("IC", "IC", 720, -acos(-1.), acos(-1.), 600, -3., 3.) ;
  IC->GetYaxis()->SetTitle("eta") ;
  IC->GetXaxis()->SetTitle("phi") ;  
  TH2F * tpgFactor = new TH2F("tpgFactor", "tpgFactor", 720, -acos(-1.), acos(-1.), 600, -3., 3.) ;
  tpgFactor->GetYaxis()->SetTitle("eta") ;
  tpgFactor->GetXaxis()->SetTitle("phi") ;  

  TH1F * hshapeEB = new TH1F("shapeEB", "shapeEB", 250, 0., 10.) ;
  TH1F * hshapeEE = new TH1F("shapeEE", "shapeEE", 250, 0., 10.) ;

  TTree * ntuple = new TTree("tpgmap","TPG geometry map") ;
  const std::string branchFloat[26] = {"fed","tcc","tower","stripInTower","xtalInStrip",
                       "CCU","VFE","xtalInVFE","xtalInCCU","ieta","iphi",
                       "ix","iy","iz","hashedId","ic","cmsswId","dbId","ietaTT","iphiTT",
                       "TCCch","TCCslot","SLBch","SLBslot","ietaGCT","iphiGCT"} ;
  ntupleInts_ = new Int_t[26] ;
  for (int i=0 ; i<26 ; i++) ntuple->Branch(branchFloat[i].c_str(),&ntupleInts_[i],(branchFloat[i]+string("/I")).c_str()) ;
  ntuple->Branch("det",ntupleDet_,"det/C") ;
  ntuple->Branch("crate",ntupleCrate_,"crate/C") ;


  TNtuple *ntupleSpike = new TNtuple("spikeParam","Spike parameters","gainId:theta:G:g:ped:pedLin") ;




  // Initialization section //
  list<uint32_t> towerListEB ;
  list<uint32_t> stripListEB ;
  list<uint32_t> towerListEE ;
  list<uint32_t> stripListEE ;
  list<uint32_t>::iterator itList ;
  
  map<int, uint32_t> stripMapEB ; // <EcalLogicId.hashed, strip elec id>
  map<uint32_t, uint32_t> stripMapEBsintheta ; // <strip elec id, sintheta>


  // Pedestals

  EcalPedestalsMap  pedMap ;

  if(m_write_ped == 1) {
    ss <<"Getting the pedestals from offline DB...\n";

    ESHandle<EcalPedestals> pedHandle;
    evtSetup.get<EcalPedestalsRcd>().get( pedHandle );
    pedMap = pedHandle.product()->getMap() ;


    EcalPedestalsMapIterator pedIter ; 
    int nPed = 0 ;
    for (pedIter = pedMap.begin() ; pedIter != pedMap.end() && nPed<10 ; ++pedIter, nPed++) {
      EcalPedestals::Item aped = (*pedIter);
      ss <<aped.mean_x12<<", "<<aped.mean_x6<<", "<<aped.mean_x1<< "\n" ;
    }
  } else if(m_write_ped==0) {
    ss <<"Getting the pedestals from previous configuration\n";
    
    EcalPedestals peds ;

    FEConfigMainInfo fe_main_info;
    fe_main_info.setConfigTag(tag_);
    ss << "trying to read previous tag if it exists tag="<< tag_<< ".version"<<version_<< "\n";
      db_-> fetchConfigSet(&fe_main_info);
      if(fe_main_info.getPedId()>0 ) ped_conf_id_=fe_main_info.getPedId();
      
      FEConfigPedInfo fe_ped_info;
      fe_ped_info.setId(ped_conf_id_);
      db_-> fetchConfigSet(&fe_ped_info);
      std::map<EcalLogicID, FEConfigPedDat> dataset_TpgPed;
      db_->fetchDataSet(&dataset_TpgPed, &fe_ped_info);
      
      typedef std::map<EcalLogicID, FEConfigPedDat>::const_iterator CIfeped;
      EcalLogicID ecid_xt;
      FEConfigPedDat  rd_ped;
      int icells=0;
      for (CIfeped p = dataset_TpgPed.begin(); p != dataset_TpgPed.end(); p++) 
    {
      ecid_xt = p->first;
      rd_ped  = p->second;
      
      std::string ecid_name=ecid_xt.getName();
      
      // EB data
      if (ecid_name=="EB_crystal_number") {


        int sm_num=ecid_xt.getID1();
        int xt_num=ecid_xt.getID2();
        
        EBDetId ebdetid(sm_num,xt_num,EBDetId::SMCRYSTALMODE);
        EcalPedestals::Item item;
        item.mean_x1  =rd_ped.getPedMeanG1() ;
        item.mean_x6  =rd_ped.getPedMeanG6();
        item.mean_x12 =rd_ped.getPedMeanG12();
        item.rms_x1  =0.5 ;
        item.rms_x6  =1. ;
        item.rms_x12  =1.2 ;

        if(icells<10) ss << " copy the EB data " << " ped = "  << item.mean_x12<< "\n";

        peds.insert(std::make_pair(ebdetid.rawId(),item));

        ++icells;
      }else if (ecid_name=="EE_crystal_number"){
              
        // EE data
        int z=ecid_xt.getID1();
        int x=ecid_xt.getID2();
        int y=ecid_xt.getID3();
        EEDetId eedetid(x,y,z,EEDetId::XYMODE);
        EcalPedestals::Item item;
        item.mean_x1  =rd_ped.getPedMeanG1() ;
        item.mean_x6  =rd_ped.getPedMeanG6();
        item.mean_x12 =rd_ped.getPedMeanG12();
        item.rms_x1  =0.5 ;
        item.rms_x6  =1. ;
        item.rms_x12  =1.2 ;
        
        peds.insert(std::make_pair(eedetid.rawId(),item));
        ++icells;
      }
    }
      
      pedMap = peds.getMap() ;
      
      EcalPedestalsMapIterator pedIter ; 
      int nPed = 0 ;
      for (pedIter = pedMap.begin() ; pedIter != pedMap.end() && nPed<10 ; ++pedIter, nPed++) {
    EcalPedestals::Item aped = (*pedIter);
    ss << aped.mean_x12 <<", " << aped.mean_x6 << ", " << aped.mean_x1 << "\n" ;
      }
      
      

      

  } else if(m_write_ped >1) {
    ss <<"Getting the pedestals from configuration number"<< m_write_ped << "\n";



    EcalPedestals peds ;

      

      FEConfigPedInfo fe_ped_info;
      fe_ped_info.setId(m_write_ped);
      db_-> fetchConfigSet(&fe_ped_info);
      std::map<EcalLogicID, FEConfigPedDat> dataset_TpgPed;
      db_->fetchDataSet(&dataset_TpgPed, &fe_ped_info);
        
      typedef std::map<EcalLogicID, FEConfigPedDat>::const_iterator CIfeped;
      EcalLogicID ecid_xt;
      FEConfigPedDat  rd_ped;
      int icells=0;
      for (CIfeped p = dataset_TpgPed.begin(); p != dataset_TpgPed.end(); p++) 
        {
          ecid_xt = p->first;
          rd_ped  = p->second;
              
          std::string ecid_name=ecid_xt.getName();
              
          // EB data
          if (ecid_name=="EB_crystal_number") {
        if(icells<10) edm::LogInfo("TopInfo") << " copy the EB data " << " icells = " << icells << "\n";
        int sm_num=ecid_xt.getID1();
        int xt_num=ecid_xt.getID2();
              
        EBDetId ebdetid(sm_num,xt_num,EBDetId::SMCRYSTALMODE);
        EcalPedestals::Item item;
        item.mean_x1  =(unsigned int)rd_ped.getPedMeanG1() ;
        item.mean_x6  =(unsigned int)rd_ped.getPedMeanG6();
        item.mean_x12 =(unsigned int)rd_ped.getPedMeanG12();
        item.rms_x1  =0.5 ;
        item.rms_x6  =1. ;
        item.rms_x12  =1.2 ;
              
        peds.insert(std::make_pair(ebdetid.rawId(),item));
        ++icells;
          }else if (ecid_name=="EE_crystal_number"){
              
        // EE data
        int z=ecid_xt.getID1();
        int x=ecid_xt.getID2();
        int y=ecid_xt.getID3();
        EEDetId eedetid(x,y,z,EEDetId::XYMODE);
        EcalPedestals::Item item;
        item.mean_x1  =(unsigned int)rd_ped.getPedMeanG1();
        item.mean_x6  =(unsigned int)rd_ped.getPedMeanG6();
        item.mean_x12 =(unsigned int)rd_ped.getPedMeanG12();
        item.rms_x1  =0.5 ;
        item.rms_x6  =1. ;
        item.rms_x12  =1.2 ;
              
        peds.insert(std::make_pair(eedetid.rawId(),item));
        ++icells;
          }
        }

      pedMap = peds.getMap() ;

      EcalPedestalsMapIterator pedIter ; 
      int nPed = 0 ;
      for (pedIter = pedMap.begin() ; pedIter != pedMap.end() && nPed<10 ; ++pedIter, nPed++) {
        EcalPedestals::Item aped = (*pedIter);
        ss << aped.mean_x12 << ", " << aped.mean_x6 << ", " << aped.mean_x1 << "\n" ;
      }
  }

  edm::LogInfo("TopInfo") << ss.str();
  ss.str("");

  // Intercalib constants
  ss <<"Getting intercalib from offline DB...\n";
  ESHandle<EcalIntercalibConstants> pIntercalib ;
  evtSetup.get<EcalIntercalibConstantsRcd>().get(pIntercalib) ;
  const EcalIntercalibConstants * intercalib = pIntercalib.product() ;
  const EcalIntercalibConstantMap & calibMap = intercalib->getMap() ;
  EcalIntercalibConstantMap::const_iterator calIter ;
  int nCal = 0 ;
  for (calIter = calibMap.begin() ; calIter != calibMap.end() && nCal<10 ; ++calIter, nCal++) {
    ss << (*calIter)<< "\n" ;
  }  
  edm::LogInfo("TopInfo") << ss.str();
  ss.str("");
  float calibvec[1700] ;
  if (H2_) {
    edm::LogInfo("TopInfo") <<"H2: overwriting IC coef with file"<< "\n" ;
    std::ifstream calibfile("calib_sm36.txt", std::ios::out) ;  
    int idata, icry ; 
    float fdata, fcali ;
    std::string strdata ;
    if (calibfile.is_open()) {
      calibfile >> strdata >>  strdata >>  strdata >>  strdata >>  strdata ;
      while(!calibfile.eof()) {
    calibfile >> idata >> icry >> fcali >> fdata >> fdata ;
    calibvec[icry-1] = fcali ;
    if(calibfile.eof()){
      break ; // avoid last line duplication
    }
      }
    }
  }    

  // Gain Ratios
  ss << "Getting the gain ratios from offline DB...\n";
  ESHandle<EcalGainRatios> pRatio;
  evtSetup.get<EcalGainRatiosRcd>().get(pRatio);
  const EcalGainRatioMap & gainMap = pRatio.product()->getMap();
  EcalGainRatioMap::const_iterator gainIter ;
  int nGain = 0 ;
  for (gainIter = gainMap.begin() ; gainIter != gainMap.end() && nGain<10 ; ++gainIter, nGain++) {
    const EcalMGPAGainRatio & aGain = (*gainIter) ;
    ss << aGain.gain12Over6() << ", " << aGain.gain6Over1() * aGain.gain12Over6()<< "\n" ;
  }
  edm::LogInfo("TopInfo") << ss.str();
  ss.str("");

  // ADCtoGeV
  ss << "Getting the ADC to GeV from offline DB...\n";
  ESHandle<EcalADCToGeVConstant> pADCToGeV ;
  evtSetup.get<EcalADCToGeVConstantRcd>().get(pADCToGeV) ;
  const EcalADCToGeVConstant * ADCToGeV = pADCToGeV.product() ;
  xtal_LSB_EB_ = ADCToGeV->getEBValue() ;
  xtal_LSB_EE_ = ADCToGeV->getEEValue() ;
  ss <<"xtal_LSB_EB_ = "<<xtal_LSB_EB_<< "\n" ;
  ss <<"xtal_LSB_EE_ = "<<xtal_LSB_EE_<< "\n" ;
 
  vector<EcalLogicID> my_EcalLogicId;
  vector<EcalLogicID> my_TTEcalLogicId;
  vector<EcalLogicID> my_StripEcalLogicId;
  EcalLogicID my_EcalLogicId_EB;
    // Endcap identifiers
  EcalLogicID my_EcalLogicId_EE;
  vector<EcalLogicID> my_TTEcalLogicId_EE;
  vector<EcalLogicID> my_RTEcalLogicId_EE;
  vector<EcalLogicID> my_StripEcalLogicId1_EE;
  vector<EcalLogicID> my_StripEcalLogicId2_EE;
  vector<EcalLogicID> my_CrystalEcalLogicId_EE;
  vector<EcalLogicID> my_TTEcalLogicId_EB_by_TCC;
  vector<EcalLogicID> my_StripEcalLogicId_EE_strips_by_TCC;

  my_EcalLogicId_EB = db_->getEcalLogicID( "EB",EcalLogicID::NULLID,EcalLogicID::NULLID,EcalLogicID::NULLID,"EB");
  my_EcalLogicId_EE = db_->getEcalLogicID( "EE",EcalLogicID::NULLID,EcalLogicID::NULLID,EcalLogicID::NULLID,"EE");
  my_EcalLogicId = db_->getEcalLogicIDSetOrdered( "EB_crystal_number",
                          1, 36,
                          1, 1700,
                          EcalLogicID::NULLID,EcalLogicID::NULLID,
                          "EB_crystal_number",12 );
  my_TTEcalLogicId = db_->getEcalLogicIDSetOrdered( "EB_trigger_tower",
                            1, 36,
                            1, 68,
                            EcalLogicID::NULLID,EcalLogicID::NULLID,
                            "EB_trigger_tower",12 );
  my_StripEcalLogicId = db_->getEcalLogicIDSetOrdered( "EB_VFE",   1, 36,   1, 68,   1,5 ,  "EB_VFE",123 ); //last digi means ordered 1st by SM,then TT, then strip 
  ss <<"got the 3 ecal barrel logic id set\n";

  // EE crystals identifiers
  my_CrystalEcalLogicId_EE = db_->getEcalLogicIDSetOrdered("EE_crystal_number",
                               -1, 1,
                               0, 200,
                               0, 200,
                               "EE_crystal_number",123 );
  
  // EE Strip identifiers
  // DCC=601-609 TT = ~40 EEstrip = 5
  my_StripEcalLogicId1_EE = db_->getEcalLogicIDSetOrdered( "ECAL_readout_strip",   
                               601, 609,   
                               1, 100,   
                               0,5 ,  
                               "ECAL_readout_strip",123 );
  // EE Strip identifiers
  // DCC=646-654 TT = ~40 EEstrip = 5
  my_StripEcalLogicId2_EE = db_->getEcalLogicIDSetOrdered( "ECAL_readout_strip",   
                               646, 654,   
                               1, 100,   
                               0,5 ,  
                               "ECAL_readout_strip",123 );
  // ----> modif here 31/1/2011   
  // EE Strip identifiers by TCC tower strip
  // TCC=1-108 TT = ~40 EEstrip = 1-5
  my_StripEcalLogicId_EE_strips_by_TCC = db_->getEcalLogicIDSetOrdered( "EE_trigger_strip",   
                               1, 108,   
                               1, 100,   
                               1,5 ,  
                               "EE_trigger_strip",123 );

  
  // ----> modif here 31/1/2011
  // ECAL Barrel trigger towers by TCC/tower identifiers
  // TTC=38-72 TT = 1-68
  my_TTEcalLogicId_EB_by_TCC = db_->getEcalLogicIDSetOrdered( "ECAL_trigger_tower",
                               1, 108,
                               1, 68,
                               EcalLogicID::NULLID,EcalLogicID::NULLID,
                               "ECAL_trigger_tower",12 );
  
  // EE TT identifiers
  // TTC=72 TT = 1440
  my_TTEcalLogicId_EE = db_->getEcalLogicIDSetOrdered( "EE_trigger_tower",
                               1, 108,
                               1, 40,
                               EcalLogicID::NULLID,EcalLogicID::NULLID,
                               "EE_trigger_tower",12 );
  
  // EE TT identifiers    
  // TTC=72 TT = 1440
  my_RTEcalLogicId_EE = db_->getEcalLogicIDSetOrdered( "EE_readout_tower",
                               1, 1000,
                               1, 100,
                               EcalLogicID::NULLID,EcalLogicID::NULLID,
                               "EE_readout_tower",12 );
  edm::LogInfo("TopInfo") << ss.str();
  ss.str("");

  if (writeToDB_) {
    ss <<"Getting the latest ids for this tag (latest version) "<< "\n";

    FEConfigMainInfo fe_main_info;
    fe_main_info.setConfigTag(tag_);
    
      ss << "trying to read previous tag if it exists tag="<< tag_<< ".version"<<version_<< "\n";

      db_-> fetchConfigSet(&fe_main_info);
      if(fe_main_info.getPedId()>0 && ped_conf_id_ ==0 ) ped_conf_id_=fe_main_info.getPedId();
      lin_conf_id_=fe_main_info.getLinId();
      lut_conf_id_=fe_main_info.getLUTId();
      wei_conf_id_=fe_main_info.getWeiId();
      fgr_conf_id_=fe_main_info.getFgrId();
      sli_conf_id_=fe_main_info.getSliId();
      spi_conf_id_=fe_main_info.getSpiId(); //modif-alex 21/01/11
      del_conf_id_=fe_main_info.getTimId(); //modif-alex 21/01/11
      if(fe_main_info.getBxtId()>0 && bxt_conf_id_==0 ) bxt_conf_id_=fe_main_info.getBxtId();
      if(fe_main_info.getBttId()>0 && btt_conf_id_==0 ) btt_conf_id_=fe_main_info.getBttId();
      if(fe_main_info.getBstId()>0 && bst_conf_id_==0 ) bst_conf_id_=fe_main_info.getBstId();
      // those that are not written specifically in this program are propagated
      // from the previous record with the same tag and the highest version

  }

  // Compute linearization coeff section //

  map<EcalLogicID, FEConfigPedDat> pedset ;
  map<EcalLogicID, FEConfigLinDat> linset ;
  map<EcalLogicID, FEConfigLinParamDat> linparamset ;
  map<EcalLogicID, FEConfigLUTParamDat> lutparamset ;
  map<EcalLogicID, FEConfigFgrParamDat> fgrparamset ;

  map<int, linStruc> linEtaSlice ;
  map< vector<int>, linStruc > linMap ;

  // count number of strip per tower
  int NbOfStripPerTCC[108][68] ; 
  for (int i=0 ; i<108 ; i++)
    for (int j=0 ; j<68 ; j++) 
      NbOfStripPerTCC[i][j] = 0 ;
  const std::vector<DetId> & ebCells = theBarrelGeometry_->getValidDetIds(DetId::Ecal, EcalBarrel);
  const std::vector<DetId> & eeCells = theEndcapGeometry_->getValidDetIds(DetId::Ecal, EcalEndcap);
  for (vector<DetId>::const_iterator it = ebCells.begin(); it != ebCells.end(); ++it) {
    EBDetId id(*it) ;
    const EcalTrigTowerDetId towid= id.tower();
    const EcalTriggerElectronicsId elId = theMapping_->getTriggerElectronicsId(id) ;
    int tccNb = theMapping_->TCCid(towid) ;
    int towerInTCC = theMapping_->iTT(towid) ;
    int stripInTower = elId.pseudoStripId() ;  
    if (stripInTower>NbOfStripPerTCC[tccNb-1][towerInTCC-1]) NbOfStripPerTCC[tccNb-1][towerInTCC-1] = stripInTower ;
  }
  for (vector<DetId>::const_iterator it = eeCells.begin(); it != eeCells.end(); ++it) {
    EEDetId id(*it) ;
    const EcalTrigTowerDetId towid= (*eTTmap_).towerOf(id) ;
    const EcalTriggerElectronicsId elId = theMapping_->getTriggerElectronicsId(id) ;
    int tccNb = theMapping_->TCCid(towid) ;
    int towerInTCC = theMapping_->iTT(towid) ; 
    int stripInTower = elId.pseudoStripId() ;
    if (stripInTower>NbOfStripPerTCC[tccNb-1][towerInTCC-1]) NbOfStripPerTCC[tccNb-1][towerInTCC-1] = stripInTower ;
  }



  // loop on EB xtals
  if (writeToFiles_) (*out_file_)<<"COMMENT ====== barrel crystals ====== "<<std::endl ;

  // special case of eta slices
  for (vector<DetId>::const_iterator it = ebCells.begin(); it != ebCells.end(); ++it) {
    EBDetId id(*it) ;
    double theta = theBarrelGeometry_->getGeometry(id)->getPosition().theta() ;
    if (!useTransverseEnergy_) theta = acos(0.) ;
    const EcalTrigTowerDetId towid= id.tower();
    const EcalTriggerElectronicsId elId = theMapping_->getTriggerElectronicsId(id) ;
    int dccNb = theMapping_->DCCid(towid) ;
    int tccNb = theMapping_->TCCid(towid) ;
    int towerInTCC = theMapping_->iTT(towid) ; // from 1 to 68 (EB)
    int stripInTower = elId.pseudoStripId() ;  // from 1 to 5
    int xtalInStrip = elId.channelId() ;       // from 1 to 5
    const EcalElectronicsId Id = theMapping_->getElectronicsId(id) ;
    int CCUid = Id.towerId() ;
    int VFEid = Id.stripId() ;
    int xtalInVFE = Id.xtalId() ;
    int xtalWithinCCUid = 5*(VFEid-1) + xtalInVFE -1 ; // Evgueni expects [0,24]
         
    (*geomFile_)<<"dccNb = "<<dccNb<<" tccNb = "<<tccNb<<" towerInTCC = "<<towerInTCC
        <<" stripInTower = "<<stripInTower<<" xtalInStrip = "<<xtalInStrip
        <<" CCUid = "<<CCUid<<" VFEid = "<<VFEid<<" xtalInVFE = "<<xtalInVFE
        <<" xtalWithinCCUid = "<<xtalWithinCCUid<<" ieta = "<<id.ieta()<<" iphi = "<<id.iphi()
        <<" xtalhashedId = "<<id.hashedIndex()<<" xtalNb = "<<id.ic()
        <<" ietaTT = "<<towid.ieta()<<" iphiTT = "<<towid.iphi()<<endl ;

    int TCCch = towerInTCC ;
    int SLBslot = int((towerInTCC-1)/8.) + 1 ;
    int SLBch = (towerInTCC-1)%8 + 1 ;
    int cmsswId = id.rawId() ;
    int ixtal=(id.ism()-1)*1700+(id.ic()-1);
    EcalLogicID logicId =my_EcalLogicId[ixtal];
    int dbId = logicId.getLogicID() ;
    int val[] = {dccNb+600,tccNb,towerInTCC,stripInTower,
         xtalInStrip,CCUid,VFEid,xtalInVFE,xtalWithinCCUid,id.ieta(),id.iphi(),
         -999,-999,towid.ieta()/abs(towid.ieta()),id.hashedIndex(),
         id.ic(), cmsswId, dbId, 
         towid.ieta(),towid.iphi(), TCCch, getCrate(tccNb).second, SLBch, SLBslot, 
         getGCTRegionEta(towid.ieta()),getGCTRegionPhi(towid.iphi())} ;
    for (int i=0 ; i<26 ; i++) ntupleInts_[i] = val[i] ;
    
    strcpy(ntupleDet_,getDet(tccNb).c_str()) ;
    strcpy(ntupleCrate_,getCrate(tccNb).first.c_str()) ;
    ntuple->Fill() ;
    
    
    if (tccNb == 37 && stripInTower == 3 && xtalInStrip == 3 && (towerInTCC-1)%4==0) {
      int etaSlice = towid.ietaAbs() ;
      coeffStruc coeff ;
      //getCoeff(coeff, calibMap, id.rawId()) ;
      //modif-alex-27-july-2015
      string str;
      getCoeff(coeff, calibMap, laserAlphaMap, id.rawId(), str) ;
      ss << str;
      getCoeff(coeff, gainMap, id.rawId()) ;
      getCoeff(coeff, pedMap, id.rawId()) ;
      linStruc lin ;
      for (int i=0 ; i<3 ; i++) {
    int mult, shift ;
    bool ok = computeLinearizerParam(theta, coeff.gainRatio_[i], coeff.calibCoeff_, "EB", mult , shift) ;
    if (!ok) edm::LogError("TopInfo") << "unable to compute the parameters for SM="<< id.ism()<<" xt="<< id.ic()<<" " <<dec<<id.rawId()<< "\n" ;  
    else {
      lin.pedestal_[i] = coeff.pedestals_[i] ;
      lin.mult_[i] = mult ;
      lin.shift_[i] = shift ;
    }
      }

      bool ok(true) ;
      if (forcedPedestalValue_ == -2) ok = realignBaseline(lin, 0) ;
      if (!ok) ss << "SM="<< id.ism()<<" xt="<< id.ic()<<" " <<dec<<id.rawId()<< "\n" ; 
      linEtaSlice[etaSlice] = lin ; 
    }
  }

  // general case
  for (vector<DetId>::const_iterator it = ebCells.begin(); it != ebCells.end(); ++it) {
    EBDetId id(*it) ;
    double theta = theBarrelGeometry_->getGeometry(id)->getPosition().theta() ;
    if (!useTransverseEnergy_) theta = acos(0.) ;
    const EcalTrigTowerDetId towid= id.tower();
    towerListEB.push_back(towid.rawId()) ;
    const EcalTriggerElectronicsId elId = theMapping_->getTriggerElectronicsId(id) ;
    stripListEB.push_back(elId.rawId() & 0xfffffff8) ;
    int dccNb = theMapping_->DCCid(towid) ;
    //int tccNb = theMapping_->TCCid(towid) ;
    int towerInTCC = theMapping_->iTT(towid) ; // from 1 to 68 (EB)
    //int stripInTower = elId.pseudoStripId() ;  // from 1 to 5
    //int xtalInStrip = elId.channelId() ;       // from 1 to 5
    const EcalElectronicsId Id = theMapping_->getElectronicsId(id) ;
    int CCUid = Id.towerId() ;
    int VFEid = Id.stripId() ;
    int xtalInVFE = Id.xtalId() ;
    int xtalWithinCCUid = 5*(VFEid-1) + xtalInVFE -1 ; // Evgueni expects [0,24]
    int etaSlice = towid.ietaAbs() ;

    // hashed index of strip EcalLogicID:
    int hashedStripLogicID = 68*5*(id.ism()-1) + 5*(towerInTCC-1) +  (VFEid-1) ;
    stripMapEB[hashedStripLogicID] = elId.rawId() & 0xfffffff8 ;
    stripMapEBsintheta[elId.rawId() & 0xfffffff8] = SFGVB_Threshold_ + abs(int(sin(theta)*pedestal_offset_)) ;

    //modif-debug
    /*
    FEConfigFgrEEStripDat stripdebug;
    EcalLogicID thestrip_debug = my_StripEcalLogicId[hashedStripLogicID] ;
    if(towid.ieta() == 6 && towid.iphi() == 7){
      std::cout << "xtal info=" << id << " VFE=" << VFEid << std::endl;
      std::cout << "TOWER DEBUG ieta=" << towid.ieta() << " iphi=" << towid.iphi() << " SFGVB=" << SFGVB_Threshold_ + abs(int(sin(theta)*pedestal_offset_)) 
        << " dbId=" << (elId.rawId() & 0xfffffff8) << " " << hashedStripLogicID << " " << thestrip_debug.getLogicID() << std::endl;  //modif-debug
    }//EB+3 TT24
    */
    //std::cout<<std::dec<<SFGVB_Threshold_ + abs(int(sin(theta)*pedestal_offset_))<<" "<<SFGVB_Threshold_<<" "<<abs(int(sin(theta)*pedestal_offset_))<<std::endl ;

    FEConfigPedDat pedDB ;
    FEConfigLinDat linDB ;
    if (writeToFiles_) (*out_file_)<<"CRYSTAL "<<dec<<id.rawId()<<std::endl ;
    //  if (writeToDB_) logicId = db_->getEcalLogicID ("EB_crystal_number", id.ism(), id.ic()) ;

    coeffStruc coeff ;
    //getCoeff(coeff, calibMap, id.rawId()) ; 
    //modif-alex-27-july-2015
    std::string str;
    getCoeff(coeff, calibMap, laserAlphaMap, id.rawId(), str) ;
    ss << str;

    if (H2_) coeff.calibCoeff_ = calibvec[id.ic()-1] ;
    getCoeff(coeff, gainMap, id.rawId()) ;
    getCoeff(coeff, pedMap, id.rawId()) ;
    ICEB->Fill(id.iphi(), id.ieta(), coeff.calibCoeff_) ;  
    IC->Fill(theBarrelGeometry_->getGeometry(id)->getPosition().phi(), theBarrelGeometry_->getGeometry(id)->getPosition().eta(), coeff.calibCoeff_) ;  

    vector<int> xtalCCU ;
    xtalCCU.push_back(dccNb+600) ; 
    xtalCCU.push_back(CCUid) ; 
    xtalCCU.push_back(xtalWithinCCUid) ;
    xtalCCU.push_back(id.rawId()) ;
    
    // compute and fill linearization parameters
    // case of eta slice
    if (forceEtaSlice_) {
      map<int, linStruc>::const_iterator itLin = linEtaSlice.find(etaSlice);
      if (itLin != linEtaSlice.end()) {
    linMap[xtalCCU] = itLin->second ;
    if (writeToFiles_) {
      for (int i=0 ; i<3 ; i++) 
        (*out_file_) << hex <<" 0x"<<itLin->second.pedestal_[i]<<" 0x"<<itLin->second.mult_[i]<<" 0x"<<itLin->second.shift_[i]<<std::endl;
    }
    if (writeToDB_) {
      for (int i=0 ; i<3 ; i++) {
        if (i==0)  {pedDB.setPedMeanG12(itLin->second.pedestal_[i]) ; linDB.setMultX12(itLin->second.mult_[i]) ; linDB.setShift12(itLin->second.shift_[i]) ; } 
        if (i==1)  {pedDB.setPedMeanG6(itLin->second.pedestal_[i]) ; linDB.setMultX6(itLin->second.mult_[i]) ; linDB.setShift6(itLin->second.shift_[i]) ; } 
        if (i==2)  {pedDB.setPedMeanG1(itLin->second.pedestal_[i]) ; linDB.setMultX1(itLin->second.mult_[i]) ; linDB.setShift1(itLin->second.shift_[i]) ; } 
      }
    }
    float factor = float(itLin->second.mult_[0])*pow(2.,-itLin->second.shift_[0])/xtal_LSB_EB_ ;
    tpgFactorEB->Fill(id.iphi(), id.ieta(), factor) ;
    tpgFactor->Fill(theBarrelGeometry_->getGeometry(id)->getPosition().phi(), 
            theBarrelGeometry_->getGeometry(id)->getPosition().eta(), factor) ;
      }
      else ss <<"current EtaSlice = "<<etaSlice<<" not found in the EtaSlice map"<< "\n" ;
    }
    else {
      // general case
      linStruc lin ;
      int forceBase12 = 0 ;
      for (int i=0 ; i<3 ; i++) {
    int mult, shift ;
    bool ok = computeLinearizerParam(theta, coeff.gainRatio_[i], coeff.calibCoeff_, "EB", mult , shift) ;
    if (!ok) edm::LogError("TopInfo") << "unable to compute the parameters for SM="<< id.ism()<<" xt="<< id.ic()<<" " <<dec<<id.rawId()<< "\n" ;  
    else {
      //PP begin
      //  mult = 0 ; shift = 0 ;
      //  if (CCUid==1 && xtalWithinCCUid==21) {
      //    if (i==0) {mult = 0x80 ; shift = 0x3 ;}
      //    if (i==1) {mult = 0x80 ; shift = 0x2 ;}
      //    if (i==2) {mult = 0xc0 ; shift = 0x0 ;}
      //  } 
      //PP end
      double base = coeff.pedestals_[i] ;
      if (forcedPedestalValue_ == -3 && i==0) {
        double G = mult*pow(2.0,-(shift+2)) ;
        double g = G/sin(theta) ;
        // int pedestal = coeff.pedestals_[i] ;
        base = double(coeff.pedestals_[i]) - pedestal_offset_/g ;
        if (base<0.) base = 0 ;
        forceBase12 = int(base) ;
      }
      lin.pedestal_[i] = coeff.pedestals_[i] ;
      lin.mult_[i] = mult ;
      lin.shift_[i] = shift ;

//    if (xtalWithinCCUid != 14) {
//      forceBase12 = 0 ;
//      lin.pedestal_[i] = 0 ; 
//      lin.mult_[i] = 0 ;
//      lin.shift_[i] = 0 ;
//    }

    }
      }

      bool ok(true) ;
      if (forcedPedestalValue_ == -2) ok = realignBaseline(lin, 0) ;
      if (forcedPedestalValue_ == -3) ok = realignBaseline(lin, forceBase12) ;
      if (!ok) ss << "SM="<< id.ism()<<" xt="<< id.ic()<<" " <<dec<<id.rawId()<< "\n" ; 
      
      for (int i=0 ; i<3 ; i++) {      
    if (writeToFiles_) (*out_file_) << hex <<" 0x"<<lin.pedestal_[i]<<" 0x"<<lin.mult_[i]<<" 0x"<<lin.shift_[i]<<std::endl; 
    if (writeToDB_) {
      if (i==0)  {pedDB.setPedMeanG12(lin.pedestal_[i]) ; linDB.setMultX12(lin.mult_[i]) ; linDB.setShift12(lin.shift_[i]) ; } 
      if (i==1)  {pedDB.setPedMeanG6(lin.pedestal_[i]) ; linDB.setMultX6(lin.mult_[i]) ; linDB.setShift6(lin.shift_[i]) ; } 
      if (i==2)  {pedDB.setPedMeanG1(lin.pedestal_[i]) ; linDB.setMultX1(lin.mult_[i]) ; linDB.setShift1(lin.shift_[i]) ; }
    }
    if (i==0) {
      float factor = float(lin.mult_[i])*pow(2.,-lin.shift_[i])/xtal_LSB_EB_ ;
      tpgFactorEB->Fill(id.iphi(), id.ieta(), factor) ;
      tpgFactor->Fill(theBarrelGeometry_->getGeometry(id)->getPosition().phi(), 
              theBarrelGeometry_->getGeometry(id)->getPosition().eta(), factor) ;               
    }
    double G = lin.mult_[i]*pow(2.0,-(lin.shift_[i]+2)) ;
    double g = G/sin(theta) ;
    float val[] = { float(i),
            float(theta),
            float(G),
            float(g),
            float(coeff.pedestals_[i]),
            float(lin.pedestal_[i])};// first arg = gainId (0 means gain12)
    ntupleSpike->Fill(val) ;
      }
      linMap[xtalCCU] = lin ;
    }
    if (writeToDB_) {
    // 1700 crystals/SM in the ECAL barrel
      int ixtal=(id.ism()-1)*1700+(id.ic()-1);
      EcalLogicID logicId =my_EcalLogicId[ixtal];
      pedset[logicId] = pedDB ;
      linset[logicId] = linDB ; 
    }
  } //ebCells

  if (writeToDB_) {
    // EcalLogicID  of the whole barrel is: my_EcalLogicId_EB
    FEConfigLinParamDat linparam ;
    linparam.setETSat(Et_sat_EB_); 
    linparamset[my_EcalLogicId_EB] = linparam ;

    FEConfigFgrParamDat fgrparam ;
    fgrparam.setFGlowthresh(FG_lowThreshold_EB_); 
    fgrparam.setFGhighthresh(FG_highThreshold_EB_); 
    fgrparam.setFGlowratio(FG_lowRatio_EB_); 
    fgrparam.setFGhighratio(FG_highRatio_EB_);
    fgrparamset[my_EcalLogicId_EB] = fgrparam ;


    FEConfigLUTParamDat lutparam ;
    lutparam.setETSat(Et_sat_EB_); 
    lutparam.setTTThreshlow(TTF_lowThreshold_EB_); 
    lutparam.setTTThreshhigh(TTF_highThreshold_EB_); 
    lutparamset[my_EcalLogicId_EB] = lutparam ;

  }


  // loop on EE xtals
  if (writeToFiles_) (*out_file_)<<"COMMENT ====== endcap crystals ====== "<<std::endl ;
  
  // special case of eta slices
  for (vector<DetId>::const_iterator it = eeCells.begin(); it != eeCells.end(); ++it) {
    EEDetId id(*it) ;
    double theta = theEndcapGeometry_->getGeometry(id)->getPosition().theta() ;
    if (!useTransverseEnergy_) theta = acos(0.) ;
    const EcalTrigTowerDetId towid= (*eTTmap_).towerOf(id) ;
    const EcalTriggerElectronicsId elId = theMapping_->getTriggerElectronicsId(id) ;
    const EcalElectronicsId Id = theMapping_->getElectronicsId(id) ;
    int dccNb = Id.dccId() ;
    int tccNb = theMapping_->TCCid(towid) ;
    int towerInTCC = theMapping_->iTT(towid) ; 
    int stripInTower = elId.pseudoStripId() ;
    int xtalInStrip = elId.channelId() ;
    int CCUid = Id.towerId() ;
    int VFEid = Id.stripId() ;
    int xtalInVFE = Id.xtalId() ;
    int xtalWithinCCUid = 5*(VFEid-1) + xtalInVFE -1 ; // Evgueni expects [0,24]

    (*geomFile_)<<"dccNb = "<<dccNb<<" tccNb = "<<tccNb<<" towerInTCC = "<<towerInTCC
        <<" stripInTower = "<<stripInTower<<" xtalInStrip = "<<xtalInStrip
        <<" CCUid = "<<CCUid<<" VFEid = "<<VFEid<<" xtalInVFE = "<<xtalInVFE
        <<" xtalWithinCCUid = "<<xtalWithinCCUid<<" ix = "<<id.ix()<<" iy = "<<id.iy()
        <<" xtalhashedId = "<<id.hashedIndex()<<" xtalNb = "<<id.isc()
        <<" ietaTT = "<<towid.ieta()<<" iphiTT = "<<towid.iphi()<<endl ;

    int TCCch = stripInTower ;
    int SLBslot, SLBch ;
    if (towerInTCC<5) {
      SLBslot = 1 ;
      SLBch = 4 + towerInTCC ;
    }
    else {
      SLBslot = int((towerInTCC-5)/8.) + 2 ;
      SLBch = (towerInTCC-5)%8 + 1 ;
    } 
    for (int j=0 ; j<towerInTCC-1 ; j++) TCCch += NbOfStripPerTCC[tccNb-1][j] ;

    int cmsswId = id.rawId() ;
    EcalLogicID logicId ;
    int iz = id.positiveZ() ;
    if (iz ==0) iz = -1 ;
    for (int k=0; k<(int)my_CrystalEcalLogicId_EE.size(); k++) {
      int z= my_CrystalEcalLogicId_EE[k].getID1() ;
      int x= my_CrystalEcalLogicId_EE[k].getID2() ;
      int y= my_CrystalEcalLogicId_EE[k].getID3() ;
      if (id.ix()==x && id.iy()==y && iz==z) logicId=my_CrystalEcalLogicId_EE[k];
    }
    int dbId = logicId.getLogicID() ;
    
    int val[] = {dccNb+600,tccNb,towerInTCC,stripInTower,
         xtalInStrip,CCUid,VFEid,xtalInVFE,xtalWithinCCUid,-999,-999,
         id.ix(),id.iy(),towid.ieta()/abs(towid.ieta()),id.hashedIndex(),
         id.ic(),cmsswId, dbId, 
         towid.ieta(),towid.iphi(),TCCch, getCrate(tccNb).second, SLBch, SLBslot, 
         getGCTRegionEta(towid.ieta()),getGCTRegionPhi(towid.iphi())} ;
    for (int i=0 ; i<26 ; i++) ntupleInts_[i] = val[i] ;
    strcpy(ntupleDet_,getDet(tccNb).c_str()) ;
    strcpy(ntupleCrate_,getCrate(tccNb).first.c_str()) ;
    ntuple->Fill() ;
     
    if ((tccNb == 76 || tccNb == 94) && stripInTower == 1 && xtalInStrip == 3 && (towerInTCC-1)%4==0) {
      int etaSlice = towid.ietaAbs() ;
      coeffStruc coeff ;
      //getCoeff(coeff, calibMap, id.rawId()) ;
      //modif-alex-27-july-2015
      std::string str;
      getCoeff(coeff, calibMap, laserAlphaMap, id.rawId(), str) ;
      ss << str;
      getCoeff(coeff, gainMap, id.rawId()) ;
      getCoeff(coeff, pedMap, id.rawId()) ;
      linStruc lin ;
      for (int i=0 ; i<3 ; i++) {
    int mult, shift ;
    bool ok = computeLinearizerParam(theta, coeff.gainRatio_[i], coeff.calibCoeff_, "EE", mult , shift) ;
    if (!ok) edm::LogError("TopInfo") << "unable to compute the parameters for Quadrant="<< id.iquadrant()<<" xt="<< id.ic()<<" " <<dec<<id.rawId()<< "\n" ;  
    else {
      lin.pedestal_[i] = coeff.pedestals_[i] ;
      lin.mult_[i] = mult ;
      lin.shift_[i] = shift ;
    }
      }

      bool ok(true) ;
      if (forcedPedestalValue_ == -2 || forcedPedestalValue_ == -3) ok = realignBaseline(lin, 0) ;
      if (!ok) ss <<"Quadrant="<< id.iquadrant()<<" xt="<< id.ic()<<" " <<dec<<id.rawId()<< "\n" ; 

      linEtaSlice[etaSlice] = lin ;
    }
  }

  // general case
  for (vector<DetId>::const_iterator it = eeCells.begin(); it != eeCells.end(); ++it) {
    EEDetId id(*it);
    double theta = theEndcapGeometry_->getGeometry(id)->getPosition().theta() ;
    if (!useTransverseEnergy_) theta = acos(0.) ;
    const EcalTrigTowerDetId towid= (*eTTmap_).towerOf(id) ;
    const EcalTriggerElectronicsId elId = theMapping_->getTriggerElectronicsId(id) ;
    const EcalElectronicsId Id = theMapping_->getElectronicsId(id) ;
    towerListEE.push_back(towid.rawId()) ;
    // special case of towers in inner rings of EE
    if (towid.ietaAbs() == 27 || towid.ietaAbs() == 28) {
      EcalTrigTowerDetId additionalTower(towid.zside(), towid.subDet(), towid.ietaAbs(), towid.iphi()+1) ;
      towerListEE.push_back(additionalTower.rawId()) ;
    }
    stripListEE.push_back(elId.rawId() & 0xfffffff8) ;
    int dccNb = Id.dccId() ;
    //int tccNb = theMapping_->TCCid(towid) ;
    //int towerInTCC = theMapping_->iTT(towid) ;
    //int stripInTower = elId.pseudoStripId() ;
    //int xtalInStrip = elId.channelId() ;
    int CCUid = Id.towerId() ;
    int VFEid = Id.stripId() ;
    int xtalInVFE = Id.xtalId() ;
    int xtalWithinCCUid = 5*(VFEid-1) + xtalInVFE -1 ; // Evgueni expects [0,24]
    int etaSlice = towid.ietaAbs() ;

    EcalLogicID logicId ;
    FEConfigPedDat pedDB ;
    FEConfigLinDat linDB ;
    if (writeToFiles_) (*out_file_)<<"CRYSTAL "<<dec<<id.rawId()<<std::endl ;
    if (writeToDB_ && DBEE_) {
      int iz = id.positiveZ() ;
      if (iz ==0) iz = -1 ;      
      for (int k=0; k<(int)my_CrystalEcalLogicId_EE.size(); k++) {
    int z= my_CrystalEcalLogicId_EE[k].getID1() ;
    int x= my_CrystalEcalLogicId_EE[k].getID2() ;
    int y= my_CrystalEcalLogicId_EE[k].getID3() ;   
    if (id.ix()==x && id.iy()==y && iz==z) logicId=my_CrystalEcalLogicId_EE[k];
      }
    }
    
    coeffStruc coeff ;
    //getCoeff(coeff, calibMap, id.rawId()) ;
    //modif-alex-27-july-2015
    string str;
    getCoeff(coeff, calibMap, laserAlphaMap, id.rawId(), str) ;
    ss << str;
    getCoeff(coeff, gainMap, id.rawId()) ;
    getCoeff(coeff, pedMap, id.rawId()) ;
    if (id.zside()>0) ICEEPlus->Fill(id.ix(), id.iy(), coeff.calibCoeff_) ;  
    else ICEEMinus->Fill(id.ix(), id.iy(), coeff.calibCoeff_) ;  
    IC->Fill(theEndcapGeometry_->getGeometry(id)->getPosition().phi(), theEndcapGeometry_->getGeometry(id)->getPosition().eta(), coeff.calibCoeff_) ;  
  
    vector<int> xtalCCU ;
    xtalCCU.push_back(dccNb+600) ; 
    xtalCCU.push_back(CCUid) ; 
    xtalCCU.push_back(xtalWithinCCUid) ;
    xtalCCU.push_back(id.rawId()) ;

    // compute and fill linearization parameters
    // case of eta slice
    if (forceEtaSlice_) {
      map<int, linStruc>::const_iterator itLin = linEtaSlice.find(etaSlice);
      if (itLin != linEtaSlice.end()) {
    linMap[xtalCCU] = itLin->second ;
    if (writeToFiles_) {
      for (int i=0 ; i<3 ; i++) 
        (*out_file_) << hex <<" 0x"<<itLin->second.pedestal_[i]<<" 0x"<<itLin->second.mult_[i]<<" 0x"<<itLin->second.shift_[i]<<std::endl;
    }
    if (writeToDB_ && DBEE_) {
      for (int i=0 ; i<3 ; i++) {
        if (i==0)  {pedDB.setPedMeanG12(itLin->second.pedestal_[i]) ; linDB.setMultX12(itLin->second.mult_[i]) ; linDB.setShift12(itLin->second.shift_[i]) ; } 
        if (i==1)  {pedDB.setPedMeanG6(itLin->second.pedestal_[i]) ; linDB.setMultX6(itLin->second.mult_[i]) ; linDB.setShift6(itLin->second.shift_[i]) ; } 
        if (i==2)  {pedDB.setPedMeanG1(itLin->second.pedestal_[i]) ; linDB.setMultX1(itLin->second.mult_[i]) ; linDB.setShift1(itLin->second.shift_[i]) ; } 
      }
    }
    float factor = float(itLin->second.mult_[0])*pow(2.,-itLin->second.shift_[0])/xtal_LSB_EE_ ;
    if (id.zside()>0) tpgFactorEEPlus->Fill(id.ix(), id.iy(), factor) ;
    else tpgFactorEEMinus->Fill(id.ix(), id.iy(), factor) ;
    tpgFactor->Fill(theEndcapGeometry_->getGeometry(id)->getPosition().phi(), 
            theEndcapGeometry_->getGeometry(id)->getPosition().eta(), factor) ;
      }
      else ss <<"current EtaSlice = "<<etaSlice<<" not found in the EtaSlice map"<< "\n";      
    }
    else {
      // general case
      linStruc lin ;
      for (int i=0 ; i<3 ; i++) {
    int mult, shift ;
    bool ok = computeLinearizerParam(theta, coeff.gainRatio_[i], coeff.calibCoeff_, "EE", mult , shift) ;
    if (!ok) edm::LogError("TopInfo") << "unable to compute the parameters for "<<dec<<id.rawId()<< "\n" ;  
    else {
      lin.pedestal_[i] = coeff.pedestals_[i] ;
      lin.mult_[i] = mult ;
      lin.shift_[i] = shift ;
    }
      }

      bool ok(true) ;
      if (forcedPedestalValue_ == -2 || forcedPedestalValue_ == -3) ok = realignBaseline(lin, 0) ;
      if (!ok) ss <<"Quadrant="<< id.iquadrant()<<" xt="<< id.ic()<<" " <<dec<<id.rawId()<< "\n" ; 

      for (int i=0 ; i<3 ; i++) {
    if (writeToFiles_) (*out_file_) << hex <<" 0x"<<lin.pedestal_[i]<<" 0x"<<lin.mult_[i]<<" 0x"<<lin.shift_[i]<<std::endl; 
    if (writeToDB_ && DBEE_) {
      if (i==0)  {pedDB.setPedMeanG12(lin.pedestal_[i]) ; linDB.setMultX12(lin.mult_[i]) ; linDB.setShift12(lin.shift_[i]) ; } 
      if (i==1)  {pedDB.setPedMeanG6(lin.pedestal_[i]) ; linDB.setMultX6(lin.mult_[i]) ; linDB.setShift6(lin.shift_[i]) ; } 
      if (i==2)  {pedDB.setPedMeanG1(lin.pedestal_[i]) ; linDB.setMultX1(lin.mult_[i]) ; linDB.setShift1(lin.shift_[i]) ; }
    }
    if (i==0) {
      float factor = float(lin.mult_[i])*pow(2.,-lin.shift_[i])/xtal_LSB_EE_ ;
      if (id.zside()>0) tpgFactorEEPlus->Fill(id.ix(), id.iy(), factor) ;
      else tpgFactorEEMinus->Fill(id.ix(), id.iy(), factor) ;       
      tpgFactor->Fill(theEndcapGeometry_->getGeometry(id)->getPosition().phi(), 
              theEndcapGeometry_->getGeometry(id)->getPosition().eta(), factor) ;                   
    }
      }
      linMap[xtalCCU] = lin ;
    }
    if (writeToDB_ && DBEE_) {
      pedset[logicId] = pedDB ;
      linset[logicId] = linDB ;
    }
  } //eeCells
  edm::LogInfo("TopInfo") << ss.str();
  ss.str("");

  if (writeToDB_ ) {
    // EcalLogicID  of the whole barrel is: my_EcalLogicId_EB
    FEConfigLinParamDat linparam ;
    linparam.setETSat(Et_sat_EE_);
    linparamset[my_EcalLogicId_EE] = linparam ;

    FEConfigLUTParamDat lutparam ;
    lutparam.setETSat(Et_sat_EE_);
    lutparam.setTTThreshlow(TTF_lowThreshold_EE_); 
    lutparam.setTTThreshhigh(TTF_highThreshold_EE_);
    lutparamset[my_EcalLogicId_EE] = lutparam ;

 
    FEConfigFgrParamDat fgrparam ;
    fgrparam.setFGlowthresh(FG_Threshold_EE_); 
    fgrparam.setFGhighthresh(FG_Threshold_EE_); 
    fgrparamset[my_EcalLogicId_EE] = fgrparam ;

  }

  if (writeToDB_) {

    ostringstream ltag;
    ltag.str("EB_"); ltag<<Et_sat_EB_<<"_EE_"<<Et_sat_EE_;
    std::string lin_tag=ltag.str();
    ss << " LIN tag "<<lin_tag<< "\n";

    if(m_write_ped==1) {
      ped_conf_id_=db_->writeToConfDB_TPGPedestals(pedset, 1, "from_OfflineDB") ;
    } else {
      ss << "the ped id ="<<ped_conf_id_<<" will be used for the pedestals "<< "\n"; 
    }

    if(m_write_lin==1) lin_conf_id_=db_->writeToConfDB_TPGLinearCoef(linset,linparamset, 1, lin_tag) ;
  }

  // Evgueni interface
  std::ofstream evgueni("TPG_hardcoded.hh", std::ios::out) ;  
  evgueni<<"void getLinParamTPG_hardcoded(int fed, int ccu, int xtal,"<<endl ;
  evgueni<<"                        int & mult12, int & shift12, int & base12,"<<endl ;
  evgueni<<"                        int & mult6, int & shift6, int & base6,"<<endl ;
  evgueni<<"                        int & mult1, int & shift1, int & base1)"<<endl ;
  evgueni<<"{"<<endl;
  evgueni<<"  mult12 = 0 ; shift12 = 0 ; base12 = 0 ; mult6 = 0 ; shift6 = 0 ; base6 = 0 ; mult1 = 0 ; shift1 = 0 ; base1 = 0 ;"<<endl ;
  map< vector<int>, linStruc>::const_iterator itLinMap ;
  for (itLinMap = linMap.begin() ; itLinMap != linMap.end() ; itLinMap++) {
    vector<int> xtalInCCU = itLinMap->first ;
    evgueni<<"  if (fed=="<<xtalInCCU[0]<<" && ccu=="<<xtalInCCU[1]<<" && xtal=="<<xtalInCCU[2]<<") {" ;
    evgueni<<"  mult12 = "<<itLinMap->second.mult_[0]<<" ; shift12 = "<<itLinMap->second.shift_[0]<<" ; base12 = "<<itLinMap->second.pedestal_[0]<<" ; " ;
    evgueni<<"  mult6 = "<<itLinMap->second.mult_[1]<<" ; shift6 = "<<itLinMap->second.shift_[1]<<" ; base6 = "<<itLinMap->second.pedestal_[1]<<" ; " ;
    evgueni<<"  mult1 = "<<itLinMap->second.mult_[2]<<" ; shift1 = "<<itLinMap->second.shift_[2]<<" ; base1 = "<<itLinMap->second.pedestal_[2]<<" ; " ;
    evgueni<<"  return ;}" <<endl ;
  }
  evgueni<<"}" <<endl ;
  evgueni.close() ;



  // Compute weights section //

  const int NWEIGROUPS = 2 ; 
  std::vector<unsigned int> weights[NWEIGROUPS] ;

  bool useDBShape = useDBShape_;
  EBShape shapeEB(useDBShape) ; shapeEB.setEventSetup(evtSetup); // EBShape, EEShape are fetched now from DB (2018.05.22 K. Theofilatos)
  EEShape shapeEE(useDBShape) ; shapeEE.setEventSetup(evtSetup); //
  weights[0] = computeWeights(shapeEB, hshapeEB) ;
  weights[1] = computeWeights(shapeEE, hshapeEE) ;

  map<EcalLogicID, FEConfigWeightGroupDat> dataset;

  for (int igrp=0 ; igrp<NWEIGROUPS ; igrp++) {

    if (weights[igrp].size() == 5) {

      if (writeToFiles_) {
    (*out_file_) <<std::endl ;
    (*out_file_) <<"WEIGHT "<<igrp<<endl ;
    for (unsigned int sample=0 ; sample<5 ; sample++) (*out_file_) << "0x" <<hex<<weights[igrp][sample]<<" " ;
    (*out_file_)<<std::endl ; 
    (*out_file_) <<std::endl ;
      }
      if (writeToDB_) {
    ss <<"going to write the weights for groupe:"<<igrp<< "\n";
    FEConfigWeightGroupDat gut;
    gut.setWeightGroupId(igrp);
    //PP WARNING: weights order is reverted when stored in the DB 
    gut.setWeight0(weights[igrp][4] );
    gut.setWeight1(weights[igrp][3]+ 0x80); //0x80 to identify the max of the pulse in the FENIX (doesn't exist in emulator)
    gut.setWeight2(weights[igrp][2] );
    gut.setWeight3(weights[igrp][1] );
    gut.setWeight4(weights[igrp][0] );
    EcalLogicID ecid = EcalLogicID( "DUMMY", igrp,igrp); //1 dummy ID per group
    // Fill the dataset
    dataset[ecid] = gut;
      }
    }
  }

  if (writeToDB_) {

    // now we store in the DB the correspondence btw channels and groups 
    map<EcalLogicID, FEConfigWeightDat> dataset2;

    // EB loop
    for (int ich=0; ich<(int)my_StripEcalLogicId.size() ; ich++){
      FEConfigWeightDat wut;
      int igroup = 0; // this group is for EB
      wut.setWeightGroupId(igroup);
      dataset2[my_StripEcalLogicId[ich]] = wut;
    }
    
    // EE loop
    for (int ich=0; ich<(int)my_StripEcalLogicId1_EE.size() ; ich++){
      FEConfigWeightDat wut;
      int igroup = 1; // this group is for EE
      wut.setWeightGroupId(igroup);
      // Fill the dataset
      dataset2[my_StripEcalLogicId1_EE[ich]] = wut;
    }
    // EE loop 2 (we had to split the ids of EE in 2 vectors to avoid crash!)
    for (int ich=0; ich<(int)my_StripEcalLogicId2_EE.size() ; ich++){
      FEConfigWeightDat wut;
      int igroup = 1; // this group is for EE
      wut.setWeightGroupId(igroup);
      // Fill the dataset
      dataset2[my_StripEcalLogicId2_EE[ich]] = wut;
    }

    // Insert the datasets
    ostringstream wtag;
    wtag.str(""); wtag<<"Shape_NGroups_"<<NWEIGROUPS;
    std::string weight_tag=wtag.str();
    ss << " weight tag "<<weight_tag<< "\n"; 
    if (m_write_wei==1) wei_conf_id_=db_->writeToConfDB_TPGWeight(dataset, dataset2, NWEIGROUPS, weight_tag) ;
  }      
  edm::LogInfo("TopInfo") << ss.str();
  ss.str("");

  // Compute FG section //

  // barrel
  unsigned int lowRatio, highRatio, lowThreshold, highThreshold, lutFG ;
  computeFineGrainEBParameters(lowRatio, highRatio, lowThreshold, highThreshold, lutFG) ;
  if (writeToFiles_) {
    (*out_file_) <<std::endl ;
    (*out_file_) <<"FG 0"<<std::endl ;
    (*out_file_)<<hex<<"0x"<<lowThreshold<<" 0x"<<highThreshold
          <<" 0x"<<lowRatio<<" 0x"<<highRatio<<" 0x"<<lutFG
          <<std::endl ;
  }

  // endcap
  unsigned int threshold, lut_tower ;
  unsigned int lut_strip;
  computeFineGrainEEParameters(threshold, lut_strip, lut_tower) ; 

  // and here we store the fgr part

  
  if (writeToDB_) {
    ss <<"going to write the fgr "<< "\n";
    map<EcalLogicID, FEConfigFgrGroupDat> dataset;
      // we create 1 group
      int NFGRGROUPS =1; 
      for (int ich=0; ich<NFGRGROUPS; ich++){
    FEConfigFgrGroupDat gut;
    gut.setFgrGroupId(ich);
    gut.setThreshLow(lowRatio );
    gut.setThreshHigh(highRatio);
    gut.setRatioLow(lowThreshold);
    gut.setRatioHigh(highThreshold);
    gut.setLUTValue(lutFG);
    EcalLogicID ecid = EcalLogicID( "DUMMY", ich,ich);
    // Fill the dataset
    dataset[ecid] = gut; // we use any logic id but different, because it is in any case ignored... 
      }
      
      // now we store in the DB the correspondence btw channels and groups 
      map<EcalLogicID, FEConfigFgrDat> dataset2;
      // in this case I decide in a stupid way which channel belongs to which group 
      for (int ich=0; ich<(int)my_TTEcalLogicId.size() ; ich++){
    FEConfigFgrDat wut;
    int igroup=0;
    wut.setFgrGroupId(igroup);
    // Fill the dataset
    // the logic ids are ordered by SM (1,...36) and TT (1,...68)  
    // you have to calculate the right index here 
    dataset2[my_TTEcalLogicId[ich]] = wut;
      }

      // endcap loop
      for (int ich=0; ich<(int)my_RTEcalLogicId_EE.size() ; ich++){
    //  std::cout << " endcap FGR " << std::endl;
    FEConfigFgrDat wut;
    int igroup=0;
    wut.setFgrGroupId(igroup);
    // Fill the dataset
    // the logic ids are ordered by .... ?  
    // you have to calculate the right index here 
    dataset2[my_RTEcalLogicId_EE[ich]] = wut;
      }

      // endcap TT loop for the FEfgr EE Tower
      map<EcalLogicID, FEConfigFgrEETowerDat> dataset3;
      for (int ich=0; ich<(int)my_TTEcalLogicId_EE.size() ; ich++){
    FEConfigFgrEETowerDat fgreett;
    fgreett.setLutValue(lut_tower);
    dataset3[my_TTEcalLogicId_EE[ich]]=fgreett;
      }

      // endcap strip loop for the FEfgr EE strip
      // and barrel strip loop for the spike parameters (same structure than EE FGr)
      map<EcalLogicID, FEConfigFgrEEStripDat> dataset4;
      for (int ich=0; ich<(int)my_StripEcalLogicId1_EE.size() ; ich++){
    FEConfigFgrEEStripDat zut;
    zut.setThreshold(threshold);
    zut.setLutFgr(lut_strip);
    dataset4[my_StripEcalLogicId1_EE[ich]] = zut;
      }
      for (int ich=0; ich<(int)my_StripEcalLogicId2_EE.size() ; ich++){
    FEConfigFgrEEStripDat zut;
    zut.setThreshold(threshold);
    zut.setLutFgr(lut_strip);
    // Fill the dataset
    dataset4[my_StripEcalLogicId2_EE[ich]] = zut;
      }
      for (int ich=0; ich<(int)my_StripEcalLogicId.size() ; ich++){
    // EB
    FEConfigFgrEEStripDat zut;
    EcalLogicID thestrip = my_StripEcalLogicId[ich] ;
    uint32_t elStripId = stripMapEB[ich] ;
    map<uint32_t, uint32_t>::const_iterator it = stripMapEBsintheta.find(elStripId) ;
    if (it != stripMapEBsintheta.end()) zut.setThreshold(it->second);
    else {
      edm::LogError("TopInfo") <<"ERROR: strip SFGVB threshold parameter not found for that strip:"<<thestrip.getID1()<<" "<<thestrip.getID3()<<" "<<thestrip.getID3()<< "\n" ; 
      edm::LogError("TopInfo") <<" using value = "<<SFGVB_Threshold_+pedestal_offset_<< "\n" ;
      zut.setThreshold(SFGVB_Threshold_+pedestal_offset_) ;
    }
    zut.setLutFgr(SFGVB_lut_);
    // Fill the dataset
    dataset4[thestrip] = zut;
      }

      // Insert the dataset
      ostringstream wtag;
      wtag.str(""); wtag<<"FGR_"<<lutFG<<"_N_"<<NFGRGROUPS<<"_eb_"<<FG_lowThreshold_EB_<<"_EB_"<<FG_highThreshold_EB_;
      std::string weight_tag=wtag.str();
      ss << " weight tag "<<weight_tag<< "\n"; 
      if(m_write_fgr==1) fgr_conf_id_=db_->writeToConfDB_TPGFgr(dataset, dataset2,  fgrparamset,dataset3, dataset4, NFGRGROUPS, weight_tag) ;

      //modif-alex 21/01/11
      map<EcalLogicID, FEConfigSpikeDat> datasetspike; //loob EB TT
      for (int ich=0; ich<(int)my_TTEcalLogicId.size() ; ich++){
        FEConfigSpikeDat spiketh;
        spiketh.setSpikeThreshold(SFGVB_SpikeKillingThreshold_);
        datasetspike[my_TTEcalLogicId[ich]] = spiketh;
      }//loop EB TT towers    

      //modif-alex 21/01/11
      ostringstream stag;
      stag.str(""); stag<<"SpikeTh"<<SFGVB_SpikeKillingThreshold_;
      std::string spike_tag=stag.str();
      ss << " spike tag "<<spike_tag<< "\n";
      if(m_write_spi==1) spi_conf_id_=db_->writeToConfDB_Spike(datasetspike, spike_tag) ; //modif-alex 21/01/11

      //modif-alex 31/01/11
      //DELAYS EB
      map<EcalLogicID, FEConfigTimingDat> datasetdelay; // the loop goes from TCC 38 to 72 and throught the towers from 1 to 68 
      for (int ich=0; ich<(int)my_TTEcalLogicId_EB_by_TCC.size() ; ich++){
        FEConfigTimingDat delay;

    EcalLogicID logiciddelay = my_TTEcalLogicId_EB_by_TCC[ich] ;
    int id1_tcc=my_TTEcalLogicId_EB_by_TCC[ich].getID1(); // the TCC
    int id2_tt =my_TTEcalLogicId_EB_by_TCC[ich].getID2(); // the tower
    std::map<int, vector<int> >::const_iterator ittEB  = delays_EB_.find(id1_tcc); 
    std::vector<int> TimingDelaysEB = ittEB->second;

    if (ittEB != delays_EB_.end()){
      if(TimingDelaysEB[id2_tt-1] == -1){
          edm::LogError("TopInfo") << "ERROR: Barrel timing delay not specified, check file, putting default value 1" << "\n";
          delay.setTimingPar1(1); 
        }
      else delay.setTimingPar1(TimingDelaysEB[id2_tt-1]);
    } else {
      edm::LogError("TopInfo") << "ERROR:Barrel Could not find delay parameter for that trigger tower " << "\n";
      edm::LogError("TopInfo") << "Using default value = 1" << "\n";
      delay.setTimingPar1(1);
    }

    std::map<int, vector<int> >::const_iterator ittpEB  = phases_EB_.find(id1_tcc);
    std::vector<int> TimingPhasesEB = ittpEB->second;

        if (ittpEB != phases_EB_.end()){
          if(TimingPhasesEB[id2_tt-1] == -1){
        edm::LogError("TopInfo") << "ERROR: Barrel timing phase not specified, check file, putting default value 0" << "\n";
        delay.setTimingPar2(0);
      }
          else delay.setTimingPar2(TimingPhasesEB[id2_tt-1]);
        } else {
      edm::LogError("TopInfo") << "ERROR:Barrel Could not find phase parameter for that trigger tower " << "\n";
      edm::LogError("TopInfo") << "Using default value = 0" << "\n";
          delay.setTimingPar2(0);
        }

    ss << ich << " tcc=" << id1_tcc << " TT=" << id2_tt << " logicId=" << logiciddelay.getLogicID() 
          << " delay=" << TimingDelaysEB[id2_tt-1] << " phase=" << TimingPhasesEB[id2_tt-1] <<  "\n"; 
    
        //delay.setTimingPar1(1);
        //delay.setTimingPar2(2);
        datasetdelay[my_TTEcalLogicId_EB_by_TCC[ich]] = delay;
      }//loop EB TT towers

      //DELAYS EE
      int stripindex = 0;
      int tccin = 1;
      for (int ich=0; ich<(int)my_StripEcalLogicId_EE_strips_by_TCC.size() ; ich++){
        FEConfigTimingDat delay;
    //int id1_strip=my_StripEcalLogicId_EE_strips_by_TCC[ich].getID1(); // the TCC
    //int id2_strip=my_StripEcalLogicId_EE_strips_by_TCC[ich].getID2(); // the Tower
    //int id3_strip=my_StripEcalLogicId_EE_strips_by_TCC[ich].getID3(); // the strip

        EcalLogicID logiciddelay = my_StripEcalLogicId_EE_strips_by_TCC[ich] ;
        int id1_tcc=my_StripEcalLogicId_EE_strips_by_TCC[ich].getID1(); // the TCC
        int id2_tt =my_StripEcalLogicId_EE_strips_by_TCC[ich].getID2(); // the tower
        int id3_st =my_StripEcalLogicId_EE_strips_by_TCC[ich].getID3(); // the strip

    //reset strip counter
    if(id1_tcc != tccin) {tccin = id1_tcc; stripindex=0;}

    std::map<int, vector<int> >::const_iterator ittEE  = delays_EE_.find(id1_tcc);
    std::vector<int> TimingDelaysEE = ittEE->second;

        if (ittEE != delays_EE_.end()){
          if(TimingDelaysEE[stripindex] == -1){
        edm::LogError("TopInfo") << "ERROR: Endcap timing delay not specified, check file, putting default value 1" << "\n";
        delay.setTimingPar1(1);
      }
          else delay.setTimingPar1(TimingDelaysEE[stripindex]);
        } else {
      edm::LogError("TopInfo") << "ERROR:Endcap Could not find delay parameter for that trigger tower " << "\n";
      edm::LogError("TopInfo") << "Using default value = 1" << "\n";
          delay.setTimingPar1(1);
        }

    std::map<int, vector<int> >::const_iterator ittpEE  = phases_EE_.find(id1_tcc);
    std::vector<int> TimingPhasesEE = ittpEE->second;

    if (ittpEE != phases_EE_.end()){
          if(TimingPhasesEE[stripindex] == -1){
        edm::LogError("TopInfo") << "ERROR: Endcap timing phase not specified, check file, putting default value 0" << "\n";
            delay.setTimingPar2(0);
          }
          else delay.setTimingPar2(TimingPhasesEE[stripindex]);
        } else {
      edm::LogError("TopInfo") << "ERROR:Endcap Could not find phase parameter for that trigger tower " << "\n";
      edm::LogError("TopInfo") << "Using default value = 0" << "\n";
          delay.setTimingPar2(0);
        }

    ss << ich << " stripindex=" << stripindex << " tcc=" << id1_tcc << " TT=" << id2_tt << " id3_st=" << id3_st 
          << " logicId=" << logiciddelay.getLogicID()
                  << " delay=" << TimingDelaysEE[stripindex] << " phase=" << TimingPhasesEE[stripindex] <<  "\n";


        //delay.setTimingPar1(1);
        //delay.setTimingPar2(2);
        datasetdelay[my_StripEcalLogicId_EE_strips_by_TCC[ich]] = delay;
    stripindex++;
      }//loop EE strip towers

      ostringstream de_tag;
      de_tag.str(""); de_tag<<"DelaysFromFile";
      std::string delay_tag=de_tag.str();
      ss << " delay tag "<<delay_tag<< "\n";
      if(m_write_del==1) del_conf_id_=db_->writeToConfDB_Delay(datasetdelay, delay_tag) ; //modif-alex 31/01/11   


  } //write to DB

  if (writeToDB_) {
    ss <<"going to write the sliding "<< "\n";
    map<EcalLogicID, FEConfigSlidingDat> dataset;
      // in this case I decide in a stupid way which channel belongs to which group 
      for (int ich=0; ich<(int)my_StripEcalLogicId.size() ; ich++){
    FEConfigSlidingDat wut;
    wut.setSliding(sliding_);
    // Fill the dataset
    // the logic ids are ordered by SM (1,...36) , TT (1,...68) and strip (1..5) 
    // you have to calculate the right index here 
    dataset[my_StripEcalLogicId[ich]] = wut;
      }

      // endcap loop
      for (int ich=0; ich<(int)my_StripEcalLogicId1_EE.size() ; ich++){
    FEConfigSlidingDat wut;
    wut.setSliding(sliding_);
    // Fill the dataset
    // the logic ids are ordered by fed tower strip
    // you have to calculate the right index here 
    dataset[my_StripEcalLogicId1_EE[ich]] = wut;
      }
      for (int ich=0; ich<(int)my_StripEcalLogicId2_EE.size() ; ich++){
    FEConfigSlidingDat wut;
    wut.setSliding(sliding_);
    // Fill the dataset
    // the logic ids are ordered by ... ? 
    // you have to calculate the right index here 
    dataset[my_StripEcalLogicId2_EE[ich]] = wut;
      }
      

      // Insert the dataset
      ostringstream wtag;
      wtag.str(""); wtag<<"Sliding_"<<sliding_;
      std::string justatag=wtag.str();
      ss << " sliding tag "<<justatag<< "\n";
      int iov_id=0; // just a parameter ... 
      if(m_write_sli==1) sli_conf_id_=db_->writeToConfDB_TPGSliding(dataset,iov_id, justatag) ;
  }

  



  // Compute LUT section //

  int lut_EB[1024], lut_EE[1024] ;

  // barrel
  computeLUT(lut_EB, "EB") ; 
  if (writeToFiles_) {
    (*out_file_) <<std::endl ;
    (*out_file_) <<"LUT 0"<<std::endl ;
    for (int i=0 ; i<1024 ; i++) (*out_file_)<<"0x"<<hex<<lut_EB[i]<<endl ;
    (*out_file_)<<endl ;
  }
  
  // endcap
  computeLUT(lut_EE, "EE") ;
  // check first if lut_EB and lut_EE are the same
  bool newLUT(false) ;
  for (int i=0 ; i<1024 ; i++) if (lut_EE[i] != lut_EB[i]) newLUT = true ;
  if (newLUT && writeToFiles_) { 
    (*out_file_) <<std::endl ;
    (*out_file_) <<"LUT 1"<<std::endl ;
    for (int i=0 ; i<1024 ; i++) (*out_file_)<<"0x"<<hex<<lut_EE[i]<<endl ;
    (*out_file_)<<endl ;
  }



  if (writeToDB_) {



    map<EcalLogicID, FEConfigLUTGroupDat> dataset;
    // we create 1 LUT group
    int NLUTGROUPS =0; 
    int ich=0;
      FEConfigLUTGroupDat lut;
      lut.setLUTGroupId(ich);
      for (int i=0; i<1024; i++){
    lut.setLUTValue(i, lut_EB[i] );
      }
      EcalLogicID ecid = EcalLogicID( "DUMMY", ich,ich);
      // Fill the dataset
      dataset[ecid] = lut; // we use any logic id but different, because it is in any case ignored... 
      
      ich++;
      
      FEConfigLUTGroupDat lute;
      lute.setLUTGroupId(ich);
      for (int i=0; i<1024; i++){
    lute.setLUTValue(i, lut_EE[i] );
      }
      EcalLogicID ecide = EcalLogicID( "DUMMY", ich,ich);
      // Fill the dataset
      dataset[ecide] = lute; // we use any logic id but different, because it is in any case ignored...

      ich++;

      NLUTGROUPS=ich;

    // now we store in the DB the correspondence btw channels and LUT groups 
    map<EcalLogicID, FEConfigLUTDat> dataset2;
    // in this case I decide in a stupid way which channel belongs to which group 
    for (int ich=0; ich<(int)my_TTEcalLogicId.size() ; ich++){
      FEConfigLUTDat lut;
      int igroup=0;
      lut.setLUTGroupId(igroup);
      // calculate the right TT - in the vector they are ordered by SM and by TT 
      // Fill the dataset
      dataset2[my_TTEcalLogicId[ich]] = lut;
    }

    // endcap loop 
    for (int ich=0; ich<(int)my_TTEcalLogicId_EE.size() ; ich++){
      FEConfigLUTDat lut;
      int igroup=1;
      lut.setLUTGroupId(igroup);
      // calculate the right TT  
      // Fill the dataset
      dataset2[my_TTEcalLogicId_EE[ich]] = lut;
    }    

    // Insert the dataset
    ostringstream ltag;
    ltag.str(""); ltag<<LUT_option_<<"_NGroups_"<<NLUTGROUPS;
    std::string lut_tag=ltag.str();
    ss << " LUT tag "<<lut_tag<< "\n"; 
    if(m_write_lut==1) lut_conf_id_=db_->writeToConfDB_TPGLUT(dataset, dataset2,lutparamset, NLUTGROUPS, lut_tag) ;

  }

  // last we insert the FE_CONFIG_MAIN table 
 if (writeToDB_) {
   
   //int conf_id_=db_->writeToConfDB_TPGMain(ped_conf_id_,lin_conf_id_, lut_conf_id_, fgr_conf_id_, 
   //               sli_conf_id_, wei_conf_id_, bxt_conf_id_, btt_conf_id_, tag_, version_) ;
   int conf_id_=db_->writeToConfDB_TPGMain(ped_conf_id_,lin_conf_id_, lut_conf_id_, fgr_conf_id_, 
                       sli_conf_id_, wei_conf_id_, spi_conf_id_, del_conf_id_, bxt_conf_id_, btt_conf_id_, bst_conf_id_, tag_, version_) ; //modif-alex 21/01/11
   
   ss << "\n Conf ID = " << conf_id_ << "\n";

 }

  // loop on strips and associate them with  values //

  // Barrel
  stripListEB.sort() ;
  stripListEB.unique() ;
  ss <<"Number of EB strips="<<dec<<stripListEB.size()<< "\n" ;
  if (writeToFiles_) {
    (*out_file_) <<std::endl ;
    for (itList = stripListEB.begin(); itList != stripListEB.end(); itList++ ) {
      (*out_file_) <<"STRIP_EB "<<dec<<(*itList)<<endl ;
      (*out_file_) << hex << "0x" <<sliding_<<std::endl ;
      (*out_file_) <<"0" <<std::endl ;
      (*out_file_) <<"0x"<<stripMapEBsintheta[(*itList)] << " 0x" << SFGVB_lut_ <<std::endl ;
    }
  }

  // Endcap
  stripListEE.sort() ;
  stripListEE.unique() ;
  ss <<"Number of EE strips="<<dec<<stripListEE.size()<< "\n" ;
  if (writeToFiles_) {
    (*out_file_) <<std::endl ;
    for (itList = stripListEE.begin(); itList != stripListEE.end(); itList++ ) {
      (*out_file_) <<"STRIP_EE "<<dec<<(*itList)<<endl ;
      (*out_file_) << hex << "0x" <<sliding_<<std::endl ;
      //(*out_file_) <<" 0" << std::endl ;
      (*out_file_) <<" 1" << std::endl ; //modif-debug to get the correct EE TPG
      (*out_file_)<<hex<<"0x"<<threshold<<" 0x"<<lut_strip<<std::endl ;  
    }
  }
  edm::LogInfo("TopInfo") << ss.str();
  ss.str("");


  // loop on towers and associate them with default values //

  // Barrel
  towerListEB.sort() ;
  towerListEB.unique() ;
  ss <<"Number of EB towers="<<dec<<towerListEB.size()<< "\n" ;
  if (writeToFiles_) {
    (*out_file_) <<std::endl ;
    for (itList = towerListEB.begin(); itList != towerListEB.end(); itList++ ) {
      (*out_file_) <<"TOWER_EB "<<dec<<(*itList)<<endl ;
      (*out_file_) <<" 0\n 0\n" ;
      (*out_file_) <<" " <<  SFGVB_SpikeKillingThreshold_ << std::endl; //modif-alex
    }
  }

  // Endcap
  towerListEE.sort() ;
  towerListEE.unique() ;
  ss <<"Number of EE towers="<<dec<<towerListEE.size()<< "\n" ;
  if (writeToFiles_) {
    (*out_file_) <<std::endl ;
    for (itList = towerListEE.begin(); itList != towerListEE.end(); itList++ ) {
      (*out_file_) <<"TOWER_EE "<<dec<<(*itList)<<endl ;
      if (newLUT) (*out_file_) <<" 1\n" ;
      else  (*out_file_) <<" 0\n" ;
      (*out_file_)<<hex<<"0x"<<lut_tower<<std::endl ;
    }
  }
  edm::LogInfo("TopInfo") << ss.str();
  ss.str("");



  // store control histos //
  ICEB->Write() ;
  tpgFactorEB->Write() ;
  ICEEPlus->Write() ;
  tpgFactorEEPlus->Write() ;  
  ICEEMinus->Write() ;
  tpgFactorEEMinus->Write() ;
  IC->Write() ;
  tpgFactor->Write() ;
  hshapeEB->Write() ;
  hshapeEE->Write() ;
  ntuple->Write() ;
  ntupleSpike->Write() ;
  saving.Close () ;
}

void EcalTPGParamBuilder::beginJob()
{
  using namespace edm;
  using namespace std;

  edm::LogInfo("TopInfo") <<"we are in beginJob\n";

  create_header() ; 

  DetId eb(DetId::Ecal,EcalBarrel) ;
  DetId ee(DetId::Ecal,EcalEndcap) ;

  if (writeToFiles_) {
    (*out_file_)<<"PHYSICS_EB "<<dec<<eb.rawId()<<std::endl ;
    (*out_file_)<<Et_sat_EB_<<" "<<TTF_lowThreshold_EB_<<" "<<TTF_highThreshold_EB_<<std::endl ;
    (*out_file_)<<FG_lowThreshold_EB_<<" "<<FG_highThreshold_EB_<<" "
          <<FG_lowRatio_EB_<<" "<<FG_highRatio_EB_<<std::endl ;
    //(*out_file_) << SFGVB_SpikeKillingThreshold_ << std::endl; //modif-alex02/02/2011
    (*out_file_) <<std::endl ;

    (*out_file_)<<"PHYSICS_EE "<<dec<<ee.rawId()<<std::endl ;
    (*out_file_)<<Et_sat_EE_<<" "<<TTF_lowThreshold_EE_<<" "<<TTF_highThreshold_EE_<<std::endl ;
    (*out_file_)<<FG_Threshold_EE_<<" "<<-1<<" "
          <<-1<<" "<<-1<<std::endl ;
    (*out_file_) <<std::endl ;
  }

}



bool EcalTPGParamBuilder::computeLinearizerParam(double theta, double gainRatio, double calibCoeff, std::string subdet, int & mult , int & shift) 
{
  /*
    Linearization coefficient are determined in order to satisfy:
    tpg(ADC_sat) = 1024
    where: 
    tpg() is a model of the linearized tpg response on 10b 
    ADC_sat is the number of ADC count corresponding the Et_sat, the maximum scale of the transverse energy
    
    Since we have:
    Et_sat = xtal_LSB * ADC_sat * gainRatio * calibCoeff * sin(theta)
    and a simple model of tpg() being given by:
    tpg(X) = [ (X*mult) >> (shift+2) ] >> (sliding+shiftDet) 
    we must satisfy:
    [ (Et_sat/(xtal_LSB * gainRatio * calibCoeff * sin(theta)) * mult) >> (shift+2) ] >> (sliding+shiftDet) = 1024 
    that is:
    mult = 1024/Et_sat * xtal_LSB * gainRatio * calibCoeff * sin(theta) * 2^-(sliding+shiftDet+2) * 2^-shift
    mult = factor * 2^-shift
  */

  // case barrel:
  int shiftDet = 2 ; //fixed, due to FE FENIX TCP format
  double ratio = xtal_LSB_EB_/Et_sat_EB_ ;
  // case endcap:
  if (subdet=="EE") {
    shiftDet = 2 ; //applied in TCC-EE and not in FE FENIX TCP... This parameters is setable in the TCC-EE
    //shiftDet = 0 ; //was like this before with FE bug
    ratio = xtal_LSB_EE_/Et_sat_EE_ ;
  }

  //modif-alex-30/01/2012
  //std::cout << "calibCoeff="<<calibCoeff<<endl; 

  double factor = 1024 * ratio * gainRatio * calibCoeff * sin(theta) * (1 << (sliding_ + shiftDet + 2)) ;
  // Let's try first with shift = 0 (trivial solution)
  mult = (int)(factor+0.5) ; 
  for (shift = 0 ; shift<15 ; shift++) {
    if (mult>=128  && mult<256) return true ;
    factor *= 2 ; 
    mult = (int)(factor+0.5) ;
  }
  edm::LogError("TopInfo") << "too bad we did not manage to calculate the factor for calib=" << calibCoeff << "\n";
  return false ;
}

void EcalTPGParamBuilder::create_header() 
{
  if (!writeToFiles_) return ;
  (*out_file_) <<"COMMENT put your comments here"<<std::endl ; 

  (*out_file_) <<"COMMENT ================================="<<std::endl ;
  (*out_file_) <<"COMMENT           physics EB structure"<<std::endl ;
  (*out_file_) <<"COMMENT"<<std::endl ;
  (*out_file_) <<"COMMENT  EtSaturation (GeV), ttf_threshold_Low (GeV), ttf_threshold_High (GeV)"<<std::endl ;
  (*out_file_) <<"COMMENT  FG_lowThreshold (GeV), FG_highThreshold (GeV), FG_lowRatio, FG_highRatio"<<std::endl ;
  //(*out_file_) <<"COMMENT  SFGVB_SpikeKillingThreshold (GeV)"<<std::endl ; //modif-alex-02/02/2011
  (*out_file_) <<"COMMENT ================================="<<std::endl ;
  (*out_file_) <<"COMMENT"<<std::endl ;

  (*out_file_) <<"COMMENT ================================="<<std::endl ;
  (*out_file_) <<"COMMENT           physics EE structure"<<std::endl ;
  (*out_file_) <<"COMMENT"<<std::endl ;
  (*out_file_) <<"COMMENT  EtSaturation (GeV), ttf_threshold_Low (GeV), ttf_threshold_High (GeV)"<<std::endl ;
  (*out_file_) <<"COMMENT  FG_Threshold (GeV), dummy, dummy, dummy"<<std::endl ;
  (*out_file_) <<"COMMENT ================================="<<std::endl ;
  (*out_file_) <<"COMMENT"<<std::endl ;

  (*out_file_) <<"COMMENT ================================="<<std::endl ;
  (*out_file_) <<"COMMENT           crystal structure (same for EB and EE)"<<std::endl ;
  (*out_file_) <<"COMMENT"<<std::endl ;
  (*out_file_) <<"COMMENT  ped, mult, shift [gain12]"<<std::endl ;
  (*out_file_) <<"COMMENT  ped, mult, shift [gain6]"<<std::endl ;
  (*out_file_) <<"COMMENT  ped, mult, shift [gain1]"<<std::endl ;
  (*out_file_) <<"COMMENT ================================="<<std::endl ;
  (*out_file_) <<"COMMENT"<<std::endl ;

  (*out_file_) <<"COMMENT ================================="<<std::endl ;
  (*out_file_) <<"COMMENT           strip EB structure"<<std::endl ;
  (*out_file_) <<"COMMENT"<<std::endl ;
  (*out_file_) <<"COMMENT  sliding_window"<<std::endl ;
  (*out_file_) <<"COMMENT  weightGroupId"<<std::endl ;
  (*out_file_) <<"COMMENT  threshold_sfg lut_sfg"<<std::endl ;
  (*out_file_) <<"COMMENT ================================="<<std::endl ;
  (*out_file_) <<"COMMENT"<<std::endl ;

  (*out_file_) <<"COMMENT ================================="<<std::endl ;
  (*out_file_) <<"COMMENT           strip EE structure"<<std::endl ;
  (*out_file_) <<"COMMENT"<<std::endl ;
  (*out_file_) <<"COMMENT  sliding_window"<<std::endl ;
  (*out_file_) <<"COMMENT  weightGroupId"<<std::endl ;
  (*out_file_) <<"COMMENT  threshold_fg lut_fg"<<std::endl ;
  (*out_file_) <<"COMMENT ================================="<<std::endl ;
  (*out_file_) <<"COMMENT"<<std::endl ;

  (*out_file_) <<"COMMENT ================================="<<std::endl ;
  (*out_file_) <<"COMMENT           tower EB structure"<<std::endl ;
  (*out_file_) <<"COMMENT"<<std::endl ;
  (*out_file_) <<"COMMENT  LUTGroupId"<<std::endl ;
  (*out_file_) <<"COMMENT  FgGroupId"<<std::endl ;
  (*out_file_) <<"COMMENT  spike_killing_threshold"<<std::endl ;//modif alex
  (*out_file_) <<"COMMENT ================================="<<std::endl ;
  (*out_file_) <<"COMMENT"<<std::endl ;

  (*out_file_) <<"COMMENT ================================="<<std::endl ;
  (*out_file_) <<"COMMENT           tower EE structure"<<std::endl ;
  (*out_file_) <<"COMMENT"<<std::endl ;
  (*out_file_) <<"COMMENT  LUTGroupId"<<std::endl ;
  (*out_file_) <<"COMMENT  tower_lut_fg"<<std::endl ;
  (*out_file_) <<"COMMENT ================================="<<std::endl ;
  (*out_file_) <<"COMMENT"<<std::endl ;
  
  (*out_file_) <<"COMMENT ================================="<<std::endl ;
  (*out_file_) <<"COMMENT           Weight structure"<<std::endl ;
  (*out_file_) <<"COMMENT"<<std::endl ;
  (*out_file_) <<"COMMENT  weightGroupId"<<std::endl ;
  (*out_file_) <<"COMMENT  w0, w1, w2, w3, w4"<<std::endl ;
  (*out_file_) <<"COMMENT ================================="<<std::endl ;
  (*out_file_) <<"COMMENT"<<std::endl ;

  (*out_file_) <<"COMMENT ================================="<<std::endl ;
  (*out_file_) <<"COMMENT           lut structure"<<std::endl ;
  (*out_file_) <<"COMMENT"<<std::endl ;
  (*out_file_) <<"COMMENT  LUTGroupId"<<std::endl ;
  (*out_file_) <<"COMMENT  LUT[1-1024]"<<std::endl ;
  (*out_file_) <<"COMMENT ================================="<<std::endl ;
  (*out_file_) <<"COMMENT"<<std::endl ;

  (*out_file_) <<"COMMENT ================================="<<std::endl ;
  (*out_file_) <<"COMMENT           fg EB structure"<<std::endl ;
  (*out_file_) <<"COMMENT"<<std::endl ;
  (*out_file_) <<"COMMENT  FgGroupId"<<std::endl ;
  (*out_file_) <<"COMMENT  el, eh, tl, th, lut_fg"<<std::endl ;
  (*out_file_) <<"COMMENT ================================="<<std::endl ;
  (*out_file_) <<"COMMENT"<<std::endl ;

  (*out_file_) <<std::endl ;
}


int EcalTPGParamBuilder::uncodeWeight(double weight, int complement2)
{
  int iweight ;
  unsigned int max = (unsigned int)(pow(2.,complement2)-1) ;
  if (weight>0) iweight=int((1<<6)*weight+0.5) ; // +0.5 for rounding pb
  else iweight= max - int(-weight*(1<<6)+0.5) +1 ;
  iweight = iweight & max ;
  return iweight ;
}

double EcalTPGParamBuilder::uncodeWeight(int iweight, int complement2) {
  double weight = double(iweight)/pow(2., 6.) ;
  // test if negative weight:
  if ( (iweight & (1<<(complement2-1))) != 0) weight = (double(iweight)-pow(2., complement2))/pow(2., 6.) ;
  return weight ;
}

std::vector<unsigned int> EcalTPGParamBuilder::computeWeights(EcalShapeBase & shape, TH1F * histo) {
  std::ostringstream ss; 
  ss <<"Computing Weights...\n" ;
  double timeMax = shape.timeOfMax() - shape.timeOfThr() ; // timeMax w.r.t begining of pulse
  double max = shape(timeMax) ;

  double sumf = 0. ;
  double sumf2 = 0. ;
  for (unsigned int sample = 0 ; sample<nSample_ ; sample++) {
    double time = timeMax - ((double)sampleMax_-(double)sample)*25. ;
    time -= weight_timeShift_ ;
    sumf += shape(time)/max ;
    sumf2 += shape(time)/max * shape(time)/max ;
    for (int subtime = 0 ; subtime<25 ; subtime++) histo->Fill(float(sample*25. + subtime)/25., shape(time+subtime)) ;
  }
  double lambda = 1./(sumf2-sumf*sumf/nSample_) ;
  double gamma = -lambda*sumf/nSample_ ;
  double * weight = new double[nSample_] ;
  for (unsigned int sample = 0 ; sample<nSample_ ; sample++) {
    double time = timeMax - ((double)sampleMax_-(double)sample)*25. ;
    time -= weight_timeShift_ ;
    weight[sample] = lambda*shape(time)/max + gamma ;
  }


  int * iweight = new int[nSample_] ;
  for (unsigned int sample = 0 ; sample<nSample_ ; sample++)   iweight[sample] = uncodeWeight(weight[sample], complement2_) ;

  // Let's check:  
  int isumw  = 0 ;  
  for (unsigned int sample = 0 ; sample<nSample_ ; sample++) isumw  += iweight[sample] ;
  unsigned int imax = (unsigned int)(pow(2.,int(complement2_))-1) ;
  isumw = (isumw & imax ) ;

  double ampl = 0. ;
  double sumw = 0. ;
  for (unsigned int sample = 0 ; sample<nSample_ ; sample++) {
     double time = timeMax - ((double)sampleMax_-(double)sample)*25. ;
     time -= weight_timeShift_ ;
     ampl += weight[sample]*shape(time) ;
     sumw += weight[sample] ;
     ss <<"weight="<<weight[sample]<<" shape="<<shape(time)<< "\n" ;
  }
  ss <<"Weights: sum="<<isumw<<" in float ="<<uncodeWeight(isumw, complement2_)<<" sum of floats ="<<sumw<< "\n" ;
  ss <<"Weights: sum (weight*shape) = "<<ampl<< "\n" ;



  // Let's correct for bias if any
  if (weight_unbias_recovery_) {
    int count = 0 ;
    while (isumw != 0 && count<10) {
      double min = 99. ;
      unsigned int index = 0 ;
      if ( (isumw & (1<<(complement2_-1))) != 0) {
    // add 1:
    ss << "Correcting for bias: adding 1\n" ;
    for (unsigned int sample = 0 ; sample<nSample_ ; sample++) {
      int new_iweight = iweight[sample]+1 ; 
      double new_weight = uncodeWeight(new_iweight, complement2_) ;
      if (fabs(new_weight-weight[sample])<min) {
        min = fabs(new_weight-weight[sample]) ;
        index = sample ;
      }
    }
    iweight[index] ++ ; 
      } else {
    // Sub 1:
    ss <<"Correcting for bias: subtracting 1\n" ;
    for (unsigned int sample = 0 ; sample<nSample_ ; sample++) {
      int new_iweight = iweight[sample]-1 ;    
      double new_weight = uncodeWeight(new_iweight, complement2_) ;
      if (fabs(new_weight-weight[sample])<min) {
        min = fabs(new_weight-weight[sample]) ;
        index = sample ;
      }
    }
    iweight[index] -- ; 
      } 
      isumw  = 0 ;  
      for (unsigned int sample = 0 ; sample<nSample_ ; sample++) isumw  += iweight[sample] ; 
      imax = (unsigned int)(pow(2.,int(complement2_))-1) ;
      isumw = (isumw & imax ) ;
      ss <<"Correcting weight number: "<< index << " sum weights = " << isumw<< "\n" ;
      count ++ ;
    }
  }
  
  // let's check again
  isumw  = 0 ;  
  for (unsigned int sample = 0 ; sample<nSample_ ; sample++) isumw  += iweight[sample] ;
  imax = (unsigned int)(pow(2.,int(complement2_))-1) ;
  isumw = (isumw & imax ) ;
  ampl = 0. ;
  for (unsigned int sample = 0 ; sample<nSample_ ; sample++) {
     double time = timeMax - ((double)sampleMax_-(double)sample)*25. ;
     time -= weight_timeShift_ ;
     double new_weight = uncodeWeight(iweight[sample], complement2_) ;
     sumw += uncodeWeight(iweight[sample], complement2_) ;
     ampl += new_weight*shape(time) ;
     ss <<"weight unbiased after integer conversion="<<new_weight<<" shape="<<shape(time)<< "\n" ;
  }
  ss <<"Weights: sum="<<isumw<<" in float ="<<uncodeWeight(isumw, complement2_)<<" sum of floats ="<<sumw<< "\n" ;
  ss <<"Weights: sum (weight*shape) = "<<ampl<< "\n" ;
  edm::LogInfo("TopInfo") << ss.str();

  std::vector<unsigned int> theWeights ;
  for (unsigned int sample = 0 ; sample<nSample_ ; sample++) theWeights.push_back(iweight[sample]) ;

  delete[] weight ;
  delete[] iweight ;
  return theWeights ;
}

void EcalTPGParamBuilder::computeLUT(int * lut, std::string det) 
{
  double Et_sat = Et_sat_EB_ ;
  double LUT_threshold = LUT_threshold_EB_ ;
  double LUT_stochastic = LUT_stochastic_EB_ ;
  double LUT_noise = LUT_noise_EB_ ;
  double LUT_constant = LUT_constant_EB_ ;
  double TTF_lowThreshold = TTF_lowThreshold_EB_ ;
  double TTF_highThreshold = TTF_highThreshold_EB_ ;
  if (det == "EE") {
    Et_sat = Et_sat_EE_ ;
    LUT_threshold = LUT_threshold_EE_ ;
    LUT_stochastic = LUT_stochastic_EE_ ;
    LUT_noise = LUT_noise_EE_ ;
    LUT_constant = LUT_constant_EE_ ;
    TTF_lowThreshold = TTF_lowThreshold_EE_ ;
    TTF_highThreshold = TTF_highThreshold_EE_ ;
  }

  // initialisation with identity
  for (int i=0 ; i<1024 ; i++) {
    lut[i] = i ;
    if (lut[i]>0xff) lut[i] = 0xff ;
  }

  // case linear LUT
  if (LUT_option_ == "Linear") {
    int mylut = 0 ;
    for (int i=0 ; i<1024 ; i++) {
      lut[i] = mylut ;
      if ((i+1)%4 == 0 ) mylut++ ; 
      //if ((i+1)%8 == 0 ) mylut++ ;//modif-alex 16/12/2010 LSB==500MeV ONLY USED FOR BEAMV4 key 
    }
  }

  // case LUT following Ecal resolution
  if (LUT_option_ == "EcalResolution") {
    TF1 * func = new TF1("func",oneOverEtResolEt, 0., Et_sat,3) ;
    func->SetParameters(LUT_stochastic, LUT_noise, LUT_constant) ;
    double norm = func->Integral(0., Et_sat) ;
    for (int i=0 ; i<1024 ; i++) {   
      double Et = i*Et_sat/1024. ;
      lut[i] =  int(0xff*func->Integral(0., Et)/norm + 0.5) ;
    }
  }

  // Now, add TTF thresholds to LUT and apply LUT threshold if needed
  for (int j=0 ; j<1024 ; j++) {
    double Et_GeV = Et_sat/1024*(j+0.5) ;
    if (Et_GeV <= LUT_threshold) lut[j] = 0 ; // LUT threshold
    int ttf = 0x0 ;    
    if (Et_GeV >= TTF_highThreshold) ttf = 3 ;
    if (Et_GeV >= TTF_lowThreshold && Et_GeV < TTF_highThreshold) ttf = 1 ;
    ttf = ttf << 8 ;
    lut[j] += ttf ;
  }

}

void EcalTPGParamBuilder::getCoeff(coeffStruc & coeff, const EcalIntercalibConstantMap & calibMap, const EcalLaserAlphaMap& laserAlphaMap,  uint rawId, std::string & st)
{
  // get current intercalibration coeff
  coeff.calibCoeff_ = 1. ;
  if (!useInterCalibration_) return ;
  EcalIntercalibConstantMap::const_iterator icalit = calibMap.find(rawId);
  
  //modif-alex-30/01/2012
  std::map<int, double >::const_iterator itCorr  = Transparency_Correction_.find(rawId);
  double icorr = 1.0;
  double alpha_factor = 1.0;

  if(useTransparencyCorr_){
    icorr = itCorr->second;
    if( itCorr != Transparency_Correction_.end() )  {
      stringstream ss;
      ss << rawId;
      stringstream ss1;
      ss1 << icorr;
      stringstream ss2;
      ss2 << (*icalit);
      st = "Transparency correction found for xtal " + ss.str() + " corr=" + ss1.str() + " intercalib=" + ss2.str() + "\n";
    }
    else
      edm::LogError("TopInfo") << "ERROR = Transparency correction not found for xtal " << rawId << "\n";

    //modif-alex-27-july-2015
    DetId ECALdetid(rawId);
    //    ss << "DETID=" << ECALdetid.subdetId() << "\n";
    stringstream ss;
    ss << ECALdetid.subdetId();
    st += "DETID=" + ss.str() + "\n";
    if(ECALdetid.subdetId() == 1){ //ECAL BARREL 
      EBDetId barrel_detid(rawId);
      EcalLaserAlphaMap::const_iterator italpha = laserAlphaMap.find(barrel_detid); 
      if ( italpha != laserAlphaMap.end() ) alpha_factor = (*italpha);
      else edm::LogError("TopInfo") << "ERROR:LaserAlphe parameter note found!!" << "\n";
    }
    if(ECALdetid.subdetId() == 1){ //ECAL ENDCAP
      EEDetId endcap_detid(rawId);
      EcalLaserAlphaMap::const_iterator italpha = laserAlphaMap.find(endcap_detid); 
      if ( italpha != laserAlphaMap.end() ) alpha_factor = (*italpha);
      else edm::LogError("TopInfo") << "ERROR:LaserAlphe parameter note found!!" << "\n";
    }

  }//transparency corrections applied

  //if( icalit != calibMap.end() ) coeff.calibCoeff_ = (*icalit) ;
  //if( icalit != calibMap.end() ) coeff.calibCoeff_ = (*icalit)/icorr; //modif-alex-30/01/2010 tansparency corrections
  //  ss << "rawId " << (*icalit) << " " << icorr << " " << alpha_factor << "\n"; 
  stringstream ss;
  ss << (*icalit);
  stringstream ss1;
  ss1 << icorr;
  stringstream ss2;
  ss2 << alpha_factor;
  st += "rawId " + ss.str() + " " + ss1.str() + " " + ss2.str() + "\n"; 
  if( icalit != calibMap.end() ) coeff.calibCoeff_ = (*icalit)/std::pow(icorr,alpha_factor); //modif-alex-27/07/2015 tansparency corrections with alpha parameters

  else edm::LogError("TopInfo") <<"getCoeff: "<<rawId<<" not found in EcalIntercalibConstantMap"<< "\n" ;
}

void EcalTPGParamBuilder::getCoeff(coeffStruc & coeff, const EcalGainRatioMap & gainMap, unsigned int rawId)
{
  // get current gain ratio
  coeff.gainRatio_[0]  = 1. ;
  coeff.gainRatio_[1]  = 2. ;
  coeff.gainRatio_[2]  = 12. ;
  EcalGainRatioMap::const_iterator gainIter = gainMap.find(rawId);
  if (gainIter != gainMap.end()) {
    const EcalMGPAGainRatio & aGain = (*gainIter) ;
    coeff.gainRatio_[1] = aGain.gain12Over6() ;
    coeff.gainRatio_[2] = aGain.gain6Over1() * aGain.gain12Over6() ;
  }
  else edm::LogError("TopInfo") <<"getCoeff: "<<rawId<<" not found in EcalGainRatioMap"<< "\n" ;
}

void EcalTPGParamBuilder::getCoeff(coeffStruc & coeff, const EcalPedestalsMap & pedMap, unsigned int rawId)
{
  coeff.pedestals_[0] = 0 ;
  coeff.pedestals_[1] = 0 ;
  coeff.pedestals_[2] = 0 ;

  if (forcedPedestalValue_ >= 0) {
    coeff.pedestals_[0] = forcedPedestalValue_ ;
    coeff.pedestals_[1] = forcedPedestalValue_ ;
    coeff.pedestals_[2] = forcedPedestalValue_ ;  
    return ;
  }

  // get current pedestal
  EcalPedestalsMapIterator pedIter = pedMap.find(rawId);
  if (pedIter != pedMap.end()) {
    EcalPedestals::Item aped = (*pedIter);
    coeff.pedestals_[0] = int(aped.mean_x12 + 0.5) ; 
    coeff.pedestals_[1] = int(aped.mean_x6 + 0.5) ;
    coeff.pedestals_[2] = int(aped.mean_x1 + 0.5) ;
  }
  else edm::LogError("TopInfo") << "getCoeff: " << rawId << " not found in EcalPedestalsMap\n" ;
}

void EcalTPGParamBuilder::getCoeff(coeffStruc & coeff, const map<EcalLogicID, MonPedestalsDat> & pedMap, const EcalLogicID & logicId)
{
  // get current pedestal
  coeff.pedestals_[0] = 0 ;
  coeff.pedestals_[1] = 0 ;
  coeff.pedestals_[2] = 0 ;

  map<EcalLogicID, MonPedestalsDat>::const_iterator it =  pedMap.find(logicId);
  if (it != pedMap.end()) {
    MonPedestalsDat ped = it->second ;
    coeff.pedestals_[0] = int(ped.getPedMeanG12() + 0.5) ; 
    coeff.pedestals_[1] = int(ped.getPedMeanG6() + 0.5) ; 
    coeff.pedestals_[2] = int(ped.getPedMeanG1() + 0.5) ; 
  } 
  else edm::LogError("TopInfo") <<"getCoeff: "<<logicId.getID1()<<", "<<logicId.getID2()<<", "<<logicId.getID3()
        <<" not found in map<EcalLogicID, MonPedestalsDat\n" ;
}

void EcalTPGParamBuilder::computeFineGrainEBParameters(unsigned int & lowRatio, unsigned int & highRatio,
                               unsigned int & lowThreshold, unsigned int & highThreshold, unsigned int & lut)
{
  lowRatio = int(0x80*FG_lowRatio_EB_ + 0.5) ;
  if (lowRatio>0x7f) lowRatio = 0x7f ;
  highRatio = int(0x80*FG_highRatio_EB_ + 0.5) ;
  if (highRatio>0x7f) highRatio = 0x7f ;
  
  // lsb at the stage of the FG calculation is:
  double lsb_FG = Et_sat_EB_/1024./4 ;
  lowThreshold = int(FG_lowThreshold_EB_/lsb_FG+0.5) ;
  if (lowThreshold>0xff) lowThreshold = 0xff ;
  highThreshold = int(FG_highThreshold_EB_/lsb_FG+0.5) ;
  if (highThreshold>0xff) highThreshold = 0xff ;

  // FG lut: FGVB response is LUT(adress) where adress is: 
  // bit3: maxof2/ET >= lowRatio, bit2: maxof2/ET >= highRatio, bit1: ET >= lowThreshold, bit0: ET >= highThreshold
  // FGVB =1 if jet-like (veto active), =0 if E.M.-like
  // the condition for jet-like is: ET>Threshold and  maxof2/ET < Ratio (only TT with enough energy are vetoed)

  // With the following lut, what matters is only max(TLow, Thigh) and max(Elow, Ehigh)
  // So, jet-like if maxof2/ettot<max(TLow, Thigh) && ettot >= max(Elow, Ehigh)
  if (FG_lut_EB_ == 0) lut = 0x0888 ; 
  else lut = FG_lut_EB_ ; // let's use the users value (hope he/she knows what he/she does!)
}

void EcalTPGParamBuilder::computeFineGrainEEParameters(unsigned int & threshold, unsigned int & lut_strip, unsigned int & lut_tower) 
{
  // lsb for EE:
  double lsb_FG = Et_sat_EE_/1024. ; // FIXME is it true????
  threshold = int(FG_Threshold_EE_/lsb_FG+0.5) ;
  lut_strip = FG_lut_strip_EE_  ;
  lut_tower = FG_lut_tower_EE_  ;
}

bool EcalTPGParamBuilder::realignBaseline(linStruc & lin, float forceBase12)
{
  bool ok(true) ;
  float base[3] = {forceBase12, float(lin.pedestal_[1]), float(lin.pedestal_[2])} ;
  for (int i=1 ; i<3 ; i++) {
    if (lin.mult_[i]>0) {
      base[i] = float(lin.pedestal_[i]) - 
    float(lin.mult_[0])/float(lin.mult_[i])*pow(2., -(lin.shift_[0]-lin.shift_[i]))*(lin.pedestal_[0]-base[0]) ;
    } else base[i] = 0 ;
  }

  for (int i=0 ; i<3 ; i++) {
    lin.pedestal_[i] = base[i] ;
    //cout<<lin.pedestal_[i]<<" "<<base[i]<<endl ;
    if (base[i]<0 || lin.pedestal_[i]>1000) {
      edm::LogError("TopInfo") <<"WARNING: base= "<<base[i]<<", "<<lin.pedestal_[i]<<" for gainId[0-2]="<<i<<" ==> forcing at 0"<< "\n" ;
      lin.pedestal_[i] = 0 ; 
      ok = false ;
    }
  }
  return ok ;

}

int EcalTPGParamBuilder::getGCTRegionPhi(int ttphi){
  int gctphi=0;
  gctphi = (ttphi+1)/4;
  if(ttphi<=2) gctphi=0;
  if(ttphi>=71) gctphi=0;
  return gctphi;
}

int EcalTPGParamBuilder::getGCTRegionEta(int tteta){
  int gcteta = 0;
  if(tteta>0) gcteta = (tteta-1)/4 + 11;
  else if(tteta<0) gcteta = (tteta+1)/4 + 10;
  return gcteta;
}

std::string EcalTPGParamBuilder::getDet(int tcc){
  std::stringstream sdet ;

  if (tcc>36 && tcc<55) sdet<<"EB-"<<tcc-36 ;
  else if (tcc>=55 && tcc<73) sdet<<"EB+"<<tcc-54 ;
  else if (tcc<=36) sdet<<"EE-" ;
  else sdet<<"EE+" ;
  
  if (tcc<=36 || tcc>=73) {
    if (tcc>=73) tcc-=72 ;
    if (tcc==1 || tcc==18 || tcc==19 || tcc==36)  sdet<<7 ;
    else if (tcc==2 || tcc==3 || tcc==20 || tcc==21)  sdet<<8 ;
    else if (tcc==4 || tcc==5 || tcc==22 || tcc==23)  sdet<<9 ;
    else if (tcc==6 || tcc==7 || tcc==24 || tcc==25)  sdet<<1 ;
    else if (tcc==8 || tcc==9 || tcc==26 || tcc==27)  sdet<<2 ;
    else if (tcc==10 || tcc==11 || tcc==28 || tcc==29)  sdet<<3 ;
    else if (tcc==12 || tcc==13 || tcc==30 || tcc==31)  sdet<<4 ;
    else if (tcc==14 || tcc==15 || tcc==32 || tcc==33)  sdet<<5 ;
    else if (tcc==16 || tcc==17 || tcc==34 || tcc==35)  sdet<<6 ;
  }

  return sdet.str() ;
}

std::pair < std::string, int >  EcalTPGParamBuilder::getCrate(int tcc){
  std::stringstream crate ;
  std::string pos ;
  int slot = 0 ;

  crate<<"S2D" ;  
  if (tcc>=40 && tcc<=42) {crate<<"02d" ; slot = 5 + (tcc-40)*6 ;}
  if (tcc>=43 && tcc<=45) {crate<<"03d" ; slot = 5 + (tcc-43)*6 ;}
  if (tcc>=37 && tcc<=39) {crate<<"04d" ; slot = 5 + (tcc-37)*6 ;}
  if (tcc>=52 && tcc<=54) {crate<<"06d" ; slot = 5 + (tcc-52)*6 ;}
  if (tcc>=46 && tcc<=48) {crate<<"07d" ; slot = 5 + (tcc-46)*6 ;}
  if (tcc>=49 && tcc<=51) {crate<<"08d" ; slot = 5 + (tcc-49)*6 ;}
  if (tcc>=58 && tcc<=60) {crate<<"02h" ; slot = 5 + (tcc-58)*6 ;}
  if (tcc>=61 && tcc<=63) {crate<<"03h" ; slot = 5 + (tcc-61)*6 ;}
  if (tcc>=55 && tcc<=57) {crate<<"04h" ; slot = 5 + (tcc-55)*6 ;}
  if (tcc>=70 && tcc<=72) {crate<<"06h" ; slot = 5 + (tcc-70)*6 ;}
  if (tcc>=64 && tcc<=66) {crate<<"07h" ; slot = 5 + (tcc-64)*6 ;}
  if (tcc>=67 && tcc<=69) {crate<<"08h" ; slot = 5 + (tcc-67)*6 ;}

  if (tcc>=76 && tcc<=81) {
    crate<<"02l" ; 
    if (tcc%2==0) slot = 2 + (tcc-76)*3 ;
    else slot = 4 + (tcc-77)*3 ;
  }
  if (tcc>=94 && tcc<=99) {
    crate<<"02l" ; 
    if (tcc%2==0) slot = 3 + (tcc-94)*3 ;
    else slot = 5 + (tcc-95)*3 ;
  }

  if (tcc>=22 && tcc<=27) {
    crate<<"03l" ; 
    if (tcc%2==0) slot = 2 + (tcc-22)*3 ;
    else slot = 4 + (tcc-23)*3 ;
  }
  if (tcc>=4 && tcc<=9) {
    crate<<"03l" ; 
    if (tcc%2==0) slot = 3 + (tcc-4)*3 ;
    else slot = 5 + (tcc-5)*3 ;
  }

  if (tcc>=82 && tcc<=87) {
    crate<<"07l" ; 
    if (tcc%2==0) slot = 2 + (tcc-82)*3 ;
    else slot = 4 + (tcc-83)*3 ;
  }
  if (tcc>=100 && tcc<=105) {
    crate<<"07l" ; 
    if (tcc%2==0) slot = 3 + (tcc-100)*3 ;
    else slot = 5 + (tcc-101)*3 ;
  }

  if (tcc>=28 && tcc<=33) {
    crate<<"08l" ; 
    if (tcc%2==0) slot = 2 + (tcc-28)*3 ;
    else slot = 4 + (tcc-29)*3 ;
  }
  if (tcc>=10 && tcc<=15) {
    crate<<"08l" ; 
    if (tcc%2==0) slot = 3 + (tcc-10)*3 ;
    else slot = 5 + (tcc-11)*3 ;
  }

  if (tcc==34) {crate<<"04l" ; slot = 2 ;}
  if (tcc==16) {crate<<"04l" ; slot = 3 ;}
  if (tcc==35) {crate<<"04l" ; slot = 4 ;}
  if (tcc==17) {crate<<"04l" ; slot = 5 ;}
  if (tcc==36) {crate<<"04l" ; slot = 8 ;}
  if (tcc==18) {crate<<"04l" ; slot = 9 ;}
  if (tcc==19) {crate<<"04l" ; slot = 10 ;}
  if (tcc==1)  {crate<<"04l" ; slot = 11 ;}
  if (tcc==20) {crate<<"04l" ; slot = 14 ;}
  if (tcc==2)  {crate<<"04l" ; slot = 15 ;}
  if (tcc==21) {crate<<"04l" ; slot = 16 ;}
  if (tcc==3)  {crate<<"04l" ; slot = 17 ;}

  if (tcc==88)  {crate<<"06l" ; slot = 2 ;}
  if (tcc==106) {crate<<"06l" ; slot = 3 ;}
  if (tcc==89)  {crate<<"06l" ; slot = 4 ;}
  if (tcc==107) {crate<<"06l" ; slot = 5 ;}
  if (tcc==90)  {crate<<"06l" ; slot = 8 ;}
  if (tcc==108) {crate<<"06l" ; slot = 9 ;}
  if (tcc==73)  {crate<<"06l" ; slot = 10 ;}
  if (tcc==91)  {crate<<"06l" ; slot = 11 ;}
  if (tcc==74)  {crate<<"06l" ; slot = 14 ;}
  if (tcc==92)  {crate<<"06l" ; slot = 15 ;}
  if (tcc==75)  {crate<<"06l" ; slot = 16 ;}
  if (tcc==93)  {crate<<"06l" ; slot = 17 ;}

  return std::pair< std::string, int > (crate.str(),slot) ;
}