CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_9_patch3/src/CalibCalorimetry/EcalTPGTools/plugins/EcalTPGParamBuilder.cc

Go to the documentation of this file.
00001 #include "EcalTPGParamBuilder.h"
00002 
00003 #include "CalibCalorimetry/EcalTPGTools/plugins/EcalTPGDBApp.h"
00004 
00005 
00006 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
00007 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00008 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
00009 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00010 #include "Geometry/EcalMapping/interface/EcalElectronicsMapping.h"
00011 #include "Geometry/EcalMapping/interface/EcalMappingRcd.h"
00012 
00013 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
00014 #include "DataFormats/EcalDetId/interface/EBDetId.h"
00015 #include "DataFormats/EcalDetId/interface/EEDetId.h"
00016 
00017 #include "CondFormats/DataRecord/interface/EcalIntercalibConstantsRcd.h"
00018 #include "CondFormats/EcalObjects/interface/EcalADCToGeVConstant.h"
00019 #include "CondFormats/DataRecord/interface/EcalADCToGeVConstantRcd.h"
00020 #include "CondFormats/EcalObjects/interface/EcalMGPAGainRatio.h"
00021 #include "CondFormats/DataRecord/interface/EcalGainRatiosRcd.h"
00022 #include "CondFormats/DataRecord/interface/EcalPedestalsRcd.h"
00023 #include "CondFormats/EcalObjects/interface/EcalTPGPedestals.h"
00024 #include "CondFormats/DataRecord/interface/EcalTPGPedestalsRcd.h"
00025 
00026 #include "SimCalorimetry/EcalSimAlgos/interface/EcalSimParameterMap.h"
00027 #if (CMSSW_VERSION>=340)
00028 #include "SimCalorimetry/EcalSimAlgos/interface/EBShape.h"
00029 #include "SimCalorimetry/EcalSimAlgos/interface/EEShape.h"
00030 #endif
00031 
00032 #include <iostream>
00033 #include <string>
00034 #include <sstream>
00035 #include <vector>
00036 #include <time.h>
00037 
00038 #include <TF1.h>
00039 #include <TH2F.h>
00040 #include <TFile.h>
00041 #include <TNtuple.h>
00042 #include <iomanip>
00043 #include <fstream>
00044 
00045 using namespace std;
00046 
00047 double oneOverEtResolEt(double *x, double *par) { 
00048   double Et = x[0] ;
00049   if (Et<1e-6) return 1./par[1] ; // to avoid division by 0.
00050   double resolEt_overEt = sqrt( (par[0]/sqrt(Et))*(par[0]/sqrt(Et)) + (par[1]/Et)*(par[1]/Et) + par[2]*par[2] ) ;
00051   return 1./(Et*resolEt_overEt) ;
00052 }
00053 
00054 EcalTPGParamBuilder::EcalTPGParamBuilder(edm::ParameterSet const& pSet)
00055   : xtal_LSB_EB_(0), xtal_LSB_EE_(0), nSample_(5), complement2_(7)
00056 {
00057   ped_conf_id_=0;
00058   lin_conf_id_=0;
00059   lut_conf_id_=0;
00060   wei_conf_id_=0;
00061   fgr_conf_id_=0;
00062   sli_conf_id_=0;
00063   spi_conf_id_=0; //modif-alex 21/01/11
00064   del_conf_id_=0; //modif-alex 21/01/11
00065   bxt_conf_id_=0;
00066   btt_conf_id_=0;
00067   bst_conf_id_=0;
00068   tag_="";
00069   version_=0;
00070 
00071   m_write_ped=1;
00072   m_write_lin=1;
00073   m_write_lut=1;
00074   m_write_wei=1;
00075   m_write_fgr=1;
00076   m_write_sli=1;
00077   m_write_spi=1; //modif-alex 21/01/11
00078   m_write_del=1; //modif-alex 21/01/11
00079   m_write_bxt=1;
00080   m_write_btt=1;
00081   m_write_bst=1;
00082 
00083   writeToDB_  = pSet.getParameter<bool>("writeToDB") ;
00084   DBEE_ = pSet.getParameter<bool>("allowDBEE") ;
00085   string DBsid    = pSet.getParameter<std::string>("DBsid") ;
00086   string DBuser   = pSet.getParameter<std::string>("DBuser") ;
00087   string DBpass   = pSet.getParameter<std::string>("DBpass") ;
00088   //uint32_t DBport = pSet.getParameter<unsigned int>("DBport") ;
00089   
00090   tag_   = pSet.getParameter<std::string>("TPGtag") ;
00091   version_ = pSet.getParameter<unsigned int>("TPGversion") ;
00092 
00093   m_write_ped = pSet.getParameter<unsigned int>("TPGWritePed") ;
00094   m_write_lin = pSet.getParameter<unsigned int>("TPGWriteLin") ;
00095   m_write_lut = pSet.getParameter<unsigned int>("TPGWriteLut") ;
00096   m_write_wei = pSet.getParameter<unsigned int>("TPGWriteWei") ;
00097   m_write_fgr = pSet.getParameter<unsigned int>("TPGWriteFgr") ;
00098   m_write_sli = pSet.getParameter<unsigned int>("TPGWriteSli") ;
00099   m_write_spi = pSet.getParameter<unsigned int>("TPGWriteSpi") ; //modif-alex 21/01/11
00100   m_write_del = pSet.getParameter<unsigned int>("TPGWriteDel") ; //modif-alex 21/01/11
00101   m_write_bxt = pSet.getParameter<unsigned int>("TPGWriteBxt") ;
00102   m_write_btt = pSet.getParameter<unsigned int>("TPGWriteBtt") ;
00103   m_write_bst = pSet.getParameter<unsigned int>("TPGWriteBst") ;
00104 
00105   btt_conf_id_=m_write_btt;
00106   bxt_conf_id_=m_write_bxt;
00107   bst_conf_id_=m_write_bst;
00108 
00109   if(m_write_ped != 0 && m_write_ped != 1 ) ped_conf_id_=m_write_ped;
00110   
00111   try {
00112     if (writeToDB_) std::cout << "data will be saved with tag and version="<< tag_<< ".version"<<version_<< endl;
00113     db_ = new EcalTPGDBApp(DBsid, DBuser, DBpass) ;
00114   } catch (exception &e) {
00115     cout << "ERROR:  " << e.what() << endl;
00116   } catch (...) {
00117     cout << "Unknown error caught" << endl;
00118   }
00119 
00120   writeToFiles_ =  pSet.getParameter<bool>("writeToFiles") ;
00121   if (writeToFiles_) {
00122     std::string outFile = pSet.getParameter<std::string>("outFile") ;
00123     out_file_ = new std::ofstream(outFile.c_str(), std::ios::out) ;  
00124   }
00125   geomFile_   = new std::ofstream("geomFile.txt", std::ios::out) ;  
00126 
00127 
00128   useTransverseEnergy_ = pSet.getParameter<bool>("useTransverseEnergy") ;
00129   
00130   Et_sat_EB_ = pSet.getParameter<double>("Et_sat_EB") ;
00131   Et_sat_EE_ = pSet.getParameter<double>("Et_sat_EE") ;
00132   sliding_ = pSet.getParameter<unsigned int>("sliding") ;
00133   weight_timeShift_ = pSet.getParameter<double>("weight_timeShift") ;
00134   sampleMax_ = pSet.getParameter<unsigned int>("weight_sampleMax") ;
00135   weight_unbias_recovery_ = pSet.getParameter<bool>("weight_unbias_recovery") ;
00136 
00137   forcedPedestalValue_ = pSet.getParameter<int>("forcedPedestalValue") ;
00138   forceEtaSlice_ = pSet.getParameter<bool>("forceEtaSlice") ;
00139     
00140   LUT_option_ = pSet.getParameter<std::string>("LUT_option") ;
00141   LUT_threshold_EB_ = pSet.getParameter<double>("LUT_threshold_EB") ;
00142   LUT_threshold_EE_ = pSet.getParameter<double>("LUT_threshold_EE") ;
00143   LUT_stochastic_EB_ = pSet.getParameter<double>("LUT_stochastic_EB") ;
00144   LUT_noise_EB_ =pSet.getParameter<double>("LUT_noise_EB") ;
00145   LUT_constant_EB_ =pSet.getParameter<double>("LUT_constant_EB") ;
00146   LUT_stochastic_EE_ = pSet.getParameter<double>("LUT_stochastic_EE") ;
00147   LUT_noise_EE_ =pSet.getParameter<double>("LUT_noise_EE") ;
00148   LUT_constant_EE_ =pSet.getParameter<double>("LUT_constant_EE") ;
00149 
00150   TTF_lowThreshold_EB_ = pSet.getParameter<double>("TTF_lowThreshold_EB") ;
00151   TTF_highThreshold_EB_ = pSet.getParameter<double>("TTF_highThreshold_EB") ;
00152   TTF_lowThreshold_EE_ = pSet.getParameter<double>("TTF_lowThreshold_EE") ;
00153   TTF_highThreshold_EE_ = pSet.getParameter<double>("TTF_highThreshold_EE") ;
00154 
00155   FG_lowThreshold_EB_ = pSet.getParameter<double>("FG_lowThreshold_EB") ;
00156   FG_highThreshold_EB_ = pSet.getParameter<double>("FG_highThreshold_EB") ;
00157   FG_lowRatio_EB_ = pSet.getParameter<double>("FG_lowRatio_EB") ;
00158   FG_highRatio_EB_ = pSet.getParameter<double>("FG_highRatio_EB") ;
00159   FG_lut_EB_ = pSet.getParameter<unsigned int>("FG_lut_EB") ;
00160   FG_Threshold_EE_ = pSet.getParameter<double>("FG_Threshold_EE") ;
00161   FG_lut_strip_EE_ = pSet.getParameter<unsigned int>("FG_lut_strip_EE") ;
00162   FG_lut_tower_EE_ = pSet.getParameter<unsigned int>("FG_lut_tower_EE") ;
00163   SFGVB_Threshold_ = pSet.getParameter<unsigned int>("SFGVB_Threshold") ;
00164   SFGVB_lut_ = pSet.getParameter<unsigned int>("SFGVB_lut") ;
00165   SFGVB_SpikeKillingThreshold_ = pSet.getParameter<int>("SFGVB_SpikeKillingThreshold") ; //modif-alex 21/01/11
00166   pedestal_offset_ = pSet.getParameter<unsigned int>("pedestal_offset") ;
00167 
00168   useInterCalibration_  = pSet.getParameter<bool>("useInterCalibration") ;
00169   H2_ = pSet.getUntrackedParameter<bool>("H2",false) ;
00170 
00171   //modif-alex-23/02/2011
00172   //convert the spike killing first from GeV to ADC (10 bits)                                                                                                                                                    
00173   //depending on the saturation scale: Et_sat_EB_                                                                                                                                                                
00174   if(SFGVB_SpikeKillingThreshold_ == -1 || (SFGVB_SpikeKillingThreshold_ > Et_sat_EB_)) 
00175     SFGVB_SpikeKillingThreshold_ = 1023; //nokilling
00176   else
00177     SFGVB_SpikeKillingThreshold_ = int(SFGVB_SpikeKillingThreshold_ * 1024/Et_sat_EB_);    
00178   std::cout << "INFO:SPIKE KILLING THRESHOLD (ADC)=" << SFGVB_SpikeKillingThreshold_ << std::endl;
00179 
00180   //modif-alex-02/02/11
00181   //TIMING information 
00182   TimingDelays_EB_ = pSet.getParameter<std::string>("timing_delays_EB") ;
00183   TimingDelays_EE_ = pSet.getParameter<std::string>("timing_delays_EE") ;
00184   TimingPhases_EB_ = pSet.getParameter<std::string>("timing_phases_EB") ;
00185   TimingPhases_EE_ = pSet.getParameter<std::string>("timing_phases_EE") ;
00186 
00187   std::cout << "INFO: READING timing files" << std::endl;
00188   std::ifstream delay_eb(TimingDelays_EB_.c_str()); 
00189   if(!delay_eb) std::cout << "ERROR: File " << TimingDelays_EB_.c_str() << " could not be opened" << std::endl;
00190   std::ifstream delay_ee(TimingDelays_EE_.c_str());
00191   if(!delay_ee) std::cout << "ERROR: File " << TimingDelays_EE_.c_str() << " could not be opened" << std::endl;
00192   std::ifstream phase_eb(TimingPhases_EB_.c_str());
00193   if(!phase_eb) std::cout << "ERROR: File " << TimingPhases_EB_.c_str() << " could not be opened" << std::endl;
00194   std::ifstream phase_ee(TimingPhases_EE_.c_str());
00195   if(!phase_ee) std::cout << "ERROR: File " << TimingPhases_EE_.c_str() << " could not be opened" << std::endl;
00196   
00197   char buf[1024];
00198   //READING DELAYS EB
00199   delay_eb.getline(buf,sizeof(buf),'\n');
00200   while( delay_eb ) {
00201     std::stringstream sin(buf);
00202 
00203     int tcc; 
00204     sin >> tcc;
00205 
00206     vector<int> vec_delays_eb;
00207     for(int ieb=0; ieb<68; ++ieb){
00208       int time_delay = -1;
00209       sin >> time_delay;
00210       vec_delays_eb.push_back(time_delay);
00211       if(time_delay==-1)  std::cout << "ERROR:Barrel timing delay -1, check file" << std::endl;
00212     }
00213 
00214     if(vec_delays_eb.size()!=68)
00215       std::cout << "ERROR:Barrel timing delay wrong, not enough towers, check file" << std::endl;
00216 
00217     if (delays_EB_.find(tcc) == delays_EB_.end())
00218       delays_EB_.insert(make_pair(tcc,vec_delays_eb));
00219 
00220     cout << tcc << " "; 
00221     for(unsigned int ieb=0; ieb<vec_delays_eb.size(); ++ieb)
00222       cout << vec_delays_eb[ieb] << " ";
00223     cout << endl;
00224 
00225     delay_eb.getline(buf,sizeof(buf),'\n');
00226   }//loop delay file EB
00227   delay_eb.close();
00228 
00229   //READING PHASES EB
00230   phase_eb.getline(buf,sizeof(buf),'\n');
00231   while( phase_eb ) {
00232     std::stringstream sin(buf);
00233     int tcc;
00234     sin >> tcc;
00235 
00236     vector<int> vec_phases_eb;
00237     for(unsigned int ieb=0; ieb<68; ++ieb){
00238       int time_phase = -1;
00239       sin >> time_phase;
00240       vec_phases_eb.push_back(time_phase);
00241       if(time_phase==-1)  std::cout << "ERROR:Barrel timing phase -1, check file" << std::endl;
00242     }
00243 
00244     if(vec_phases_eb.size()!=68)
00245       std::cout << "ERROR:Barrel timing phase wrong, not enough towers, check file" << std::endl;
00246 
00247     if (phases_EB_.find(tcc) == phases_EB_.end())
00248       phases_EB_.insert(make_pair(tcc,vec_phases_eb));
00249 
00250     cout << tcc << " ";
00251     for(unsigned int ieb=0; ieb<vec_phases_eb.size(); ++ieb)
00252       cout << vec_phases_eb[ieb] << " ";
00253     cout << endl;
00254 
00255     phase_eb.getline(buf,sizeof(buf),'\n');
00256   }//loop phase file EB
00257   phase_eb.close();
00258 
00259   //READING DELAYS EE//------------------------------------------------
00260   delay_ee.getline(buf,sizeof(buf),'\n');
00261   while( delay_ee ) {
00262     std::stringstream sin(buf);
00263     int tcc; 
00264     sin >> tcc;
00265 
00266     vector<int> vec_delays_ee;
00267     for(unsigned int iee=0; iee<48; ++iee){
00268       int time_delay = -1;
00269       sin >> time_delay;
00270       vec_delays_ee.push_back(time_delay);
00271       if(time_delay==-1)  std::cout << "ERROR:EE timing delay -1, check file" << std::endl;
00272     }
00273 
00274     if(vec_delays_ee.size()!=48)
00275       std::cout << "ERROR:EE timing delay wrong, not enough towers, check file" << std::endl;
00276 
00277     if (delays_EE_.find(tcc) == delays_EE_.end())
00278       delays_EE_.insert(make_pair(tcc,vec_delays_ee));
00279 
00280     cout << tcc << " "; 
00281     for(unsigned int iee=0; iee<vec_delays_ee.size(); ++iee)
00282       cout << vec_delays_ee[iee] << " ";
00283     cout << endl;
00284 
00285     delay_ee.getline(buf,sizeof(buf),'\n');
00286   }//loop delay file EE
00287   delay_ee.close();
00288 
00289 
00290   //READING PHASES EE
00291   phase_ee.getline(buf,sizeof(buf),'\n');
00292   while( phase_ee ) {
00293     std::stringstream sin(buf);
00294     int tcc;
00295     sin >> tcc;
00296 
00297     vector<int> vec_phases_ee;
00298     for(unsigned int iee=0; iee<48; ++iee){
00299       int time_phase = -1;
00300       sin >> time_phase;
00301       vec_phases_ee.push_back(time_phase);
00302       if(time_phase==-1)  std::cout << "ERROR:EE timing phase -1, check file" << std::endl;
00303     }
00304 
00305     if(vec_phases_ee.size()!=48)
00306       std::cout << "ERROR:EE timing phase wrong, not enough towers, check file" << std::endl;
00307 
00308     if (phases_EE_.find(tcc) == phases_EE_.end())
00309       phases_EE_.insert(make_pair(tcc,vec_phases_ee));
00310 
00311     cout << tcc << " ";
00312     for(unsigned int iee=0; iee<vec_phases_ee.size(); ++iee)
00313       cout << vec_phases_ee[iee] << " ";
00314     cout << endl;
00315 
00316     phase_ee.getline(buf,sizeof(buf),'\n');
00317   }//loop phase file EE
00318   phase_ee.close();
00319 
00320   std::cout << "INFO: DONE reading timing files for EB and EE" << std::endl;
00321 
00322 }
00323 
00324 EcalTPGParamBuilder::~EcalTPGParamBuilder()
00325 { 
00326   if (writeToFiles_) {
00327     (*out_file_ )<<"EOF"<<std::endl ;
00328     out_file_->close() ;
00329     delete out_file_ ;
00330   }
00331 }
00332 
00333 
00334 bool EcalTPGParamBuilder::checkIfOK(EcalPedestals::Item item) 
00335 {
00336   bool result=true;
00337   if( item.mean_x1 <150. || item.mean_x1 >250) result=false;
00338   if( item.mean_x6 <150. || item.mean_x6 >250) result=false;
00339   if( item.mean_x12 <150. || item.mean_x12 >250) result=false;
00340   if( item.rms_x1 <0 || item.rms_x1 > 2) result=false;
00341   if( item.rms_x6 <0 || item.rms_x6 > 3) result=false;
00342   if( item.rms_x12 <0 || item.rms_x12 > 5) result=false;
00343   return result; 
00344 }
00345 
00346 int EcalTPGParamBuilder::getEtaSlice(int tccId, int towerInTCC)
00347 {
00348   int etaSlice = (towerInTCC-1)/4+1 ;
00349   // barrel
00350   if (tccId>36 || tccId<73) return etaSlice ;
00351   //endcap
00352   else {
00353     if (tccId >=1 && tccId <= 18) etaSlice += 21 ; // inner -
00354     if (tccId >=19 && tccId <= 36) etaSlice += 17 ; // outer -
00355     if (tccId >=91 && tccId <= 108) etaSlice += 21 ; // inner +
00356     if (tccId >=73 && tccId <= 90) etaSlice += 17 ; // outer +
00357   }
00358   return etaSlice ;
00359 }
00360 
00361 void EcalTPGParamBuilder::analyze(const edm::Event& evt, const edm::EventSetup& evtSetup) 
00362 {
00363   using namespace edm;
00364   using namespace std;
00365 
00366   // geometry
00367   ESHandle<CaloGeometry> theGeometry;
00368   ESHandle<CaloSubdetectorGeometry> theEndcapGeometry_handle, theBarrelGeometry_handle;
00369   evtSetup.get<CaloGeometryRecord>().get( theGeometry );
00370   evtSetup.get<EcalEndcapGeometryRecord>().get("EcalEndcap",theEndcapGeometry_handle);
00371   evtSetup.get<EcalBarrelGeometryRecord>().get("EcalBarrel",theBarrelGeometry_handle);
00372   evtSetup.get<IdealGeometryRecord>().get(eTTmap_);
00373   theEndcapGeometry_ = &(*theEndcapGeometry_handle);
00374   theBarrelGeometry_ = &(*theBarrelGeometry_handle);
00375 
00376   // electronics mapping
00377   ESHandle< EcalElectronicsMapping > ecalmapping;
00378   evtSetup.get< EcalMappingRcd >().get(ecalmapping);
00379   theMapping_ = ecalmapping.product();
00380 
00381   
00382   // histo
00383   TFile saving ("EcalTPGParam.root","recreate") ;
00384   saving.cd () ;
00385   TH2F * ICEB = new TH2F("ICEB", "IC: Barrel", 360, 1, 361, 172, -86, 86) ;
00386   ICEB->GetYaxis()->SetTitle("eta index") ;
00387   ICEB->GetXaxis()->SetTitle("phi index") ;  
00388   TH2F * tpgFactorEB = new TH2F("tpgFactorEB", "tpgFactor: Barrel", 360, 1, 361, 172, -86, 86) ;
00389   tpgFactorEB->GetYaxis()->SetTitle("eta index") ;
00390   tpgFactorEB->GetXaxis()->SetTitle("phi index") ;  
00391 
00392   TH2F * ICEEPlus = new TH2F("ICEEPlus", "IC: Plus Endcap", 120, -9, 111, 120, -9, 111) ;
00393   ICEEPlus->GetYaxis()->SetTitle("y index") ;
00394   ICEEPlus->GetXaxis()->SetTitle("x index") ;
00395   TH2F * tpgFactorEEPlus = new TH2F("tpgFactorEEPlus", "tpgFactor: Plus Endcap", 120, -9, 111, 120, -9, 111) ;
00396   tpgFactorEEPlus->GetYaxis()->SetTitle("y index") ;
00397   tpgFactorEEPlus->GetXaxis()->SetTitle("x index") ;
00398   TH2F * ICEEMinus = new TH2F("ICEEMinus", "IC: Minus Endcap", 120, -9, 111, 120, -9, 111) ;
00399   ICEEMinus->GetYaxis()->SetTitle("y index") ;
00400   ICEEMinus->GetXaxis()->SetTitle("x index") ;
00401   TH2F * tpgFactorEEMinus = new TH2F("tpgFactorEEMinus", "tpgFactor: Minus Endcap", 120, -9, 111, 120, -9, 111) ;
00402   tpgFactorEEMinus->GetYaxis()->SetTitle("y index") ;
00403   tpgFactorEEMinus->GetXaxis()->SetTitle("x index") ;
00404 
00405   TH2F * IC = new TH2F("IC", "IC", 720, -acos(-1.), acos(-1.), 600, -3., 3.) ;
00406   IC->GetYaxis()->SetTitle("eta") ;
00407   IC->GetXaxis()->SetTitle("phi") ;  
00408   TH2F * tpgFactor = new TH2F("tpgFactor", "tpgFactor", 720, -acos(-1.), acos(-1.), 600, -3., 3.) ;
00409   tpgFactor->GetYaxis()->SetTitle("eta") ;
00410   tpgFactor->GetXaxis()->SetTitle("phi") ;  
00411 
00412   TH1F * hshapeEB = new TH1F("shapeEB", "shapeEB", 250, 0., 10.) ;
00413   TH1F * hshapeEE = new TH1F("shapeEE", "shapeEE", 250, 0., 10.) ;
00414 
00415   TTree * ntuple = new TTree("tpgmap","TPG geometry map") ;
00416   const std::string branchFloat[26] = {"fed","tcc","tower","stripInTower","xtalInStrip",
00417                                        "CCU","VFE","xtalInVFE","xtalInCCU","ieta","iphi",
00418                                        "ix","iy","iz","hashedId","ic","cmsswId","dbId","ietaTT","iphiTT",
00419                                        "TCCch","TCCslot","SLBch","SLBslot","ietaGCT","iphiGCT"} ;
00420   ntupleInts_ = new Int_t[26] ;
00421   for (int i=0 ; i<26 ; i++) ntuple->Branch(branchFloat[i].c_str(),&ntupleInts_[i],(branchFloat[i]+string("/I")).c_str()) ;
00422   ntuple->Branch("det",ntupleDet_,"det/C") ;
00423   ntuple->Branch("crate",ntupleCrate_,"crate/C") ;
00424 
00425 
00426   TNtuple *ntupleSpike = new TNtuple("spikeParam","Spike parameters","gainId:theta:G:g:ped:pedLin") ;
00427 
00428 
00429 
00430 
00432   // Initialization section //
00434   list<uint32_t> towerListEB ;
00435   list<uint32_t> stripListEB ;
00436   list<uint32_t> towerListEE ;
00437   list<uint32_t> stripListEE ;
00438   list<uint32_t>::iterator itList ;
00439   
00440   map<int, uint32_t> stripMapEB ; // <EcalLogicId.hashed, strip elec id>
00441   map<uint32_t, uint32_t> stripMapEBsintheta ; // <strip elec id, sintheta>
00442 
00443 
00444   // Pedestals
00445 
00446   EcalPedestalsMap  pedMap ;
00447 
00448   if(m_write_ped == 1) {
00449     std::cout <<"Getting the pedestals from offline DB..."<<endl;
00450 
00451     ESHandle<EcalPedestals> pedHandle;
00452     evtSetup.get<EcalPedestalsRcd>().get( pedHandle );
00453     pedMap = pedHandle.product()->getMap() ;
00454 
00455 
00456     EcalPedestalsMapIterator pedIter ; 
00457     int nPed = 0 ;
00458     for (pedIter = pedMap.begin() ; pedIter != pedMap.end() && nPed<10 ; ++pedIter, nPed++) {
00459       EcalPedestals::Item aped = (*pedIter);
00460       std::cout<<aped.mean_x12<<", "<<aped.mean_x6<<", "<<aped.mean_x1<<std::endl ;
00461     }
00462   } else if(m_write_ped==0) {
00463     std::cout <<"Getting the pedestals from previous configuration"<<std::endl;
00464     
00465     EcalPedestals peds ;
00466 
00467     FEConfigMainInfo fe_main_info;
00468     fe_main_info.setConfigTag(tag_);
00469     try {
00470       std::cout << "trying to read previous tag if it exists tag="<< tag_<< ".version"<<version_<< endl;
00471       db_-> fetchConfigSet(&fe_main_info);
00472       if(fe_main_info.getPedId()>0 ) ped_conf_id_=fe_main_info.getPedId();
00473       
00474       FEConfigPedInfo fe_ped_info;
00475       fe_ped_info.setId(ped_conf_id_);
00476       db_-> fetchConfigSet(&fe_ped_info);
00477       std::map<EcalLogicID, FEConfigPedDat> dataset_TpgPed;
00478       db_->fetchDataSet(&dataset_TpgPed, &fe_ped_info);
00479       
00480       typedef std::map<EcalLogicID, FEConfigPedDat>::const_iterator CIfeped;
00481       EcalLogicID ecid_xt;
00482       FEConfigPedDat  rd_ped;
00483       int icells=0;
00484       for (CIfeped p = dataset_TpgPed.begin(); p != dataset_TpgPed.end(); p++) 
00485         {
00486           ecid_xt = p->first;
00487           rd_ped  = p->second;
00488           
00489           std::string ecid_name=ecid_xt.getName();
00490           
00491           // EB data
00492           if (ecid_name=="EB_crystal_number") {
00493 
00494 
00495             int sm_num=ecid_xt.getID1();
00496             int xt_num=ecid_xt.getID2();
00497             
00498             EBDetId ebdetid(sm_num,xt_num,EBDetId::SMCRYSTALMODE);
00499             EcalPedestals::Item item;
00500             item.mean_x1  =rd_ped.getPedMeanG1() ;
00501             item.mean_x6  =rd_ped.getPedMeanG6();
00502             item.mean_x12 =rd_ped.getPedMeanG12();
00503             item.rms_x1  =0.5 ;
00504             item.rms_x6  =1. ;
00505             item.rms_x12  =1.2 ;
00506 
00507             if(icells<10) std::cout << " copy the EB data " << " ped = "  << item.mean_x12<< std::endl;
00508 
00509             peds.insert(std::make_pair(ebdetid.rawId(),item));
00510 
00511             ++icells;
00512           }else if (ecid_name=="EE_crystal_number"){
00513                       
00514             // EE data
00515             int z=ecid_xt.getID1();
00516             int x=ecid_xt.getID2();
00517             int y=ecid_xt.getID3();
00518             EEDetId eedetid(x,y,z,EEDetId::XYMODE);
00519             EcalPedestals::Item item;
00520             item.mean_x1  =rd_ped.getPedMeanG1() ;
00521             item.mean_x6  =rd_ped.getPedMeanG6();
00522             item.mean_x12 =rd_ped.getPedMeanG12();
00523             item.rms_x1  =0.5 ;
00524             item.rms_x6  =1. ;
00525             item.rms_x12  =1.2 ;
00526             
00527             peds.insert(std::make_pair(eedetid.rawId(),item));
00528             ++icells;
00529           }
00530         }
00531       
00532       pedMap = peds.getMap() ;
00533       
00534       EcalPedestalsMapIterator pedIter ; 
00535       int nPed = 0 ;
00536       for (pedIter = pedMap.begin() ; pedIter != pedMap.end() && nPed<10 ; ++pedIter, nPed++) {
00537         EcalPedestals::Item aped = (*pedIter);
00538         std::cout<<aped.mean_x12<<", "<<aped.mean_x6<<", "<<aped.mean_x1<<std::endl ;
00539       }
00540       
00541       
00542     } catch (exception &e) {
00543       cout << " error reading previous pedestals " << endl;
00544     } catch (...) {
00545       cout << "Unknown error reading previous pedestals " << endl;
00546     }
00547     
00548 
00549       
00550 
00551   } else if(m_write_ped >1) {
00552     std::cout <<"Getting the pedestals from configuration number"<< m_write_ped <<std::endl;
00553 
00554 
00555 
00556     EcalPedestals peds ;
00557 
00558       try {
00559 
00560           FEConfigPedInfo fe_ped_info;
00561           fe_ped_info.setId(m_write_ped);
00562           db_-> fetchConfigSet(&fe_ped_info);
00563           std::map<EcalLogicID, FEConfigPedDat> dataset_TpgPed;
00564           db_->fetchDataSet(&dataset_TpgPed, &fe_ped_info);
00565             
00566           typedef std::map<EcalLogicID, FEConfigPedDat>::const_iterator CIfeped;
00567           EcalLogicID ecid_xt;
00568           FEConfigPedDat  rd_ped;
00569           int icells=0;
00570           for (CIfeped p = dataset_TpgPed.begin(); p != dataset_TpgPed.end(); p++) 
00571             {
00572               ecid_xt = p->first;
00573               rd_ped  = p->second;
00574                   
00575               std::string ecid_name=ecid_xt.getName();
00576                   
00577               // EB data
00578               if (ecid_name=="EB_crystal_number") {
00579                 if(icells<10) std::cout << " copy the EB data " << " icells = " << icells << std::endl;
00580                 int sm_num=ecid_xt.getID1();
00581                 int xt_num=ecid_xt.getID2();
00582                       
00583                 EBDetId ebdetid(sm_num,xt_num,EBDetId::SMCRYSTALMODE);
00584                 EcalPedestals::Item item;
00585                 item.mean_x1  =(unsigned int)rd_ped.getPedMeanG1() ;
00586                 item.mean_x6  =(unsigned int)rd_ped.getPedMeanG6();
00587                 item.mean_x12 =(unsigned int)rd_ped.getPedMeanG12();
00588                 item.rms_x1  =0.5 ;
00589                 item.rms_x6  =1. ;
00590                 item.rms_x12  =1.2 ;
00591                       
00592                 peds.insert(std::make_pair(ebdetid.rawId(),item));
00593                 ++icells;
00594               }else if (ecid_name=="EE_crystal_number"){
00595                       
00596                 // EE data
00597                 int z=ecid_xt.getID1();
00598                 int x=ecid_xt.getID2();
00599                 int y=ecid_xt.getID3();
00600                 EEDetId eedetid(x,y,z,EEDetId::XYMODE);
00601                 EcalPedestals::Item item;
00602                 item.mean_x1  =(unsigned int)rd_ped.getPedMeanG1();
00603                 item.mean_x6  =(unsigned int)rd_ped.getPedMeanG6();
00604                 item.mean_x12 =(unsigned int)rd_ped.getPedMeanG12();
00605                 item.rms_x1  =0.5 ;
00606                 item.rms_x6  =1. ;
00607                 item.rms_x12  =1.2 ;
00608                       
00609                 peds.insert(std::make_pair(eedetid.rawId(),item));
00610                 ++icells;
00611               }
00612             }
00613 
00614           pedMap = peds.getMap() ;
00615 
00616           EcalPedestalsMapIterator pedIter ; 
00617           int nPed = 0 ;
00618           for (pedIter = pedMap.begin() ; pedIter != pedMap.end() && nPed<10 ; ++pedIter, nPed++) {
00619             EcalPedestals::Item aped = (*pedIter);
00620             std::cout<<aped.mean_x12<<", "<<aped.mean_x6<<", "<<aped.mean_x1<<std::endl ;
00621           }
00622           
00623 
00624       } catch (exception &e) {
00625         cout << " error reading previous pedestals " << endl;
00626       } catch (...) {
00627         cout << "Unknown error reading previous pedestals " << endl;
00628       }
00629 
00630 
00631   }
00632 
00633 
00634 
00635   std::cout<<"...\n"<<std::endl ; 
00636 
00637 
00638   // Intercalib constants
00639   std::cout <<"Getting intercalib from offline DB..."<<endl;
00640   ESHandle<EcalIntercalibConstants> pIntercalib ;
00641   evtSetup.get<EcalIntercalibConstantsRcd>().get(pIntercalib) ;
00642   const EcalIntercalibConstants * intercalib = pIntercalib.product() ;
00643   const EcalIntercalibConstantMap & calibMap = intercalib->getMap() ;
00644   EcalIntercalibConstantMap::const_iterator calIter ;
00645   int nCal = 0 ;
00646   for (calIter = calibMap.begin() ; calIter != calibMap.end() && nCal<10 ; ++calIter, nCal++) {
00647     std::cout<<(*calIter)<<std::endl ;
00648   }  
00649   std::cout<<"...\n"<<std::endl ;
00650   float calibvec[1700] ;
00651   if (H2_) {
00652     std::cout<<"H2: overwriting IC coef with file"<<std::endl ;
00653     std::ifstream calibfile("calib_sm36.txt", std::ios::out) ;  
00654     int idata, icry ; 
00655     float fdata, fcali ;
00656     std::string strdata ;
00657     if (calibfile.is_open()) {
00658       calibfile >> strdata >>  strdata >>  strdata >>  strdata >>  strdata ;
00659       while(!calibfile.eof()) {
00660         calibfile >> idata >> icry >> fcali >> fdata >> fdata ;
00661         calibvec[icry-1] = fcali ;
00662         if(calibfile.eof()){
00663           break ; // avoid last line duplication
00664         }
00665       }
00666     }
00667   }    
00668 
00669   // Gain Ratios
00670   std::cout <<"Getting the gain ratios from offline DB..."<<endl;
00671   ESHandle<EcalGainRatios> pRatio;
00672   evtSetup.get<EcalGainRatiosRcd>().get(pRatio);
00673   const EcalGainRatioMap & gainMap = pRatio.product()->getMap();
00674   EcalGainRatioMap::const_iterator gainIter ;
00675   int nGain = 0 ;
00676   for (gainIter = gainMap.begin() ; gainIter != gainMap.end() && nGain<10 ; ++gainIter, nGain++) {
00677     const EcalMGPAGainRatio & aGain = (*gainIter) ;
00678     std::cout<<aGain.gain12Over6()<<", "<<aGain.gain6Over1() * aGain.gain12Over6()<<std::endl ;
00679   }  
00680   std::cout<<"...\n"<<std::endl ;    
00681 
00682 
00683   // ADCtoGeV
00684   std::cout <<"Getting the ADC to GEV from offline DB..."<<endl;
00685   ESHandle<EcalADCToGeVConstant> pADCToGeV ;
00686   evtSetup.get<EcalADCToGeVConstantRcd>().get(pADCToGeV) ;
00687   const EcalADCToGeVConstant * ADCToGeV = pADCToGeV.product() ;
00688   xtal_LSB_EB_ = ADCToGeV->getEBValue() ;
00689   xtal_LSB_EE_ = ADCToGeV->getEEValue() ;
00690   std::cout<<"xtal_LSB_EB_ = "<<xtal_LSB_EB_<<std::endl ;
00691   std::cout<<"xtal_LSB_EE_ = "<<xtal_LSB_EE_<<std::endl ;
00692   std::cout<<std::endl ;  
00693 
00694   
00695   vector<EcalLogicID> my_EcalLogicId;
00696   vector<EcalLogicID> my_TTEcalLogicId;
00697   vector<EcalLogicID> my_StripEcalLogicId;
00698   EcalLogicID my_EcalLogicId_EB;
00699     // Endcap identifiers
00700   EcalLogicID my_EcalLogicId_EE;
00701   vector<EcalLogicID> my_TTEcalLogicId_EE;
00702   vector<EcalLogicID> my_RTEcalLogicId_EE;
00703   vector<EcalLogicID> my_StripEcalLogicId1_EE;
00704   vector<EcalLogicID> my_StripEcalLogicId2_EE;
00705   vector<EcalLogicID> my_CrystalEcalLogicId_EE;
00706   vector<EcalLogicID> my_TTEcalLogicId_EB_by_TCC;
00707   vector<EcalLogicID> my_StripEcalLogicId_EE_strips_by_TCC;
00708 
00709 
00710   std::cout<<"going to get the ecal logic id set"<< endl;
00711   my_EcalLogicId_EB = db_->getEcalLogicID( "EB",EcalLogicID::NULLID,EcalLogicID::NULLID,EcalLogicID::NULLID,"EB");
00712   my_EcalLogicId_EE = db_->getEcalLogicID( "EE",EcalLogicID::NULLID,EcalLogicID::NULLID,EcalLogicID::NULLID,"EE");
00713   my_EcalLogicId = db_->getEcalLogicIDSetOrdered( "EB_crystal_number",
00714                                                   1, 36,
00715                                                   1, 1700,
00716                                                   EcalLogicID::NULLID,EcalLogicID::NULLID,
00717                                                   "EB_crystal_number",12 );
00718   my_TTEcalLogicId = db_->getEcalLogicIDSetOrdered( "EB_trigger_tower",
00719                                                     1, 36,
00720                                                     1, 68,
00721                                                     EcalLogicID::NULLID,EcalLogicID::NULLID,
00722                                                     "EB_trigger_tower",12 );
00723   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 
00724   std::cout<<"got the 3 ecal barrel logic id set"<< endl;
00725 
00726   // EE crystals identifiers
00727   my_CrystalEcalLogicId_EE = db_->getEcalLogicIDSetOrdered("EE_crystal_number",
00728                                                            -1, 1,
00729                                                            0, 200,
00730                                                            0, 200,
00731                                                            "EE_crystal_number",123 );
00732   
00733   // EE Strip identifiers
00734   // DCC=601-609 TT = ~40 EEstrip = 5
00735   my_StripEcalLogicId1_EE = db_->getEcalLogicIDSetOrdered( "ECAL_readout_strip",   
00736                                                            601, 609,   
00737                                                            1, 100,   
00738                                                            0,5 ,  
00739                                                            "ECAL_readout_strip",123 );
00740   // EE Strip identifiers
00741   // DCC=646-654 TT = ~40 EEstrip = 5
00742   my_StripEcalLogicId2_EE = db_->getEcalLogicIDSetOrdered( "ECAL_readout_strip",   
00743                                                            646, 654,   
00744                                                            1, 100,   
00745                                                            0,5 ,  
00746                                                            "ECAL_readout_strip",123 );
00747   // ----> modif here 31/1/2011   
00748   // EE Strip identifiers by TCC tower strip
00749   // TCC=1-108 TT = ~40 EEstrip = 1-5
00750   my_StripEcalLogicId_EE_strips_by_TCC = db_->getEcalLogicIDSetOrdered( "EE_trigger_strip",   
00751                                                            1, 108,   
00752                                                            1, 100,   
00753                                                            1,5 ,  
00754                                                            "EE_trigger_strip",123 );
00755 
00756   
00757   // ----> modif here 31/1/2011
00758   // ECAL Barrel trigger towers by TCC/tower identifiers
00759   // TTC=38-72 TT = 1-68
00760   my_TTEcalLogicId_EB_by_TCC = db_->getEcalLogicIDSetOrdered( "ECAL_trigger_tower",
00761                                                        1, 108,
00762                                                        1, 68,
00763                                                        EcalLogicID::NULLID,EcalLogicID::NULLID,
00764                                                        "ECAL_trigger_tower",12 );
00765   
00766   // EE TT identifiers
00767   // TTC=72 TT = 1440
00768   my_TTEcalLogicId_EE = db_->getEcalLogicIDSetOrdered( "EE_trigger_tower",
00769                                                        1, 108,
00770                                                        1, 40,
00771                                                        EcalLogicID::NULLID,EcalLogicID::NULLID,
00772                                                        "EE_trigger_tower",12 );
00773   
00774   // EE TT identifiers    
00775   // TTC=72 TT = 1440
00776   my_RTEcalLogicId_EE = db_->getEcalLogicIDSetOrdered( "EE_readout_tower",
00777                                                        1, 1000,
00778                                                        1, 100,
00779                                                        EcalLogicID::NULLID,EcalLogicID::NULLID,
00780                                                        "EE_readout_tower",12 );
00781   
00782   std::cout<<"got the end cap logic id set"<< endl;
00783 
00784   if (writeToDB_) {
00785     std::cout<<"Getting the latest ids for this tag (latest version) "<< endl;
00786 
00787     FEConfigMainInfo fe_main_info;
00788     fe_main_info.setConfigTag(tag_);
00789     try {
00790       std::cout << "trying to read previous tag if it exists tag="<< tag_<< ".version"<<version_<< endl;
00791 
00792       db_-> fetchConfigSet(&fe_main_info);
00793       if(fe_main_info.getPedId()>0 && ped_conf_id_ ==0 ) ped_conf_id_=fe_main_info.getPedId();
00794       lin_conf_id_=fe_main_info.getLinId();
00795       lut_conf_id_=fe_main_info.getLUTId();
00796       wei_conf_id_=fe_main_info.getWeiId();
00797       fgr_conf_id_=fe_main_info.getFgrId();
00798       sli_conf_id_=fe_main_info.getSliId();
00799       spi_conf_id_=fe_main_info.getSpiId(); //modif-alex 21/01/11
00800       del_conf_id_=fe_main_info.getTimId(); //modif-alex 21/01/11
00801       if(fe_main_info.getBxtId()>0 && bxt_conf_id_==0 ) bxt_conf_id_=fe_main_info.getBxtId();
00802       if(fe_main_info.getBttId()>0 && btt_conf_id_==0 ) btt_conf_id_=fe_main_info.getBttId();
00803       if(fe_main_info.getBstId()>0 && bst_conf_id_==0 ) bst_conf_id_=fe_main_info.getBstId();
00804       // those that are not written specifically in this program are propagated
00805       // from the previous record with the same tag and the highest version
00806 
00807       std::cout<<"got it "<< endl;
00808 
00809     } catch (exception &e) {
00810       cout << " tag did not exist a new tag will be created " << endl;
00811     } catch (...) {
00812       cout << "Unknown error caught" << endl;
00813     }
00814 
00815   }
00816 
00818   // Compute linearization coeff section //
00820 
00821   map<EcalLogicID, FEConfigPedDat> pedset ;
00822   map<EcalLogicID, FEConfigLinDat> linset ;
00823   map<EcalLogicID, FEConfigLinParamDat> linparamset ;
00824   map<EcalLogicID, FEConfigLUTParamDat> lutparamset ;
00825   map<EcalLogicID, FEConfigFgrParamDat> fgrparamset ;
00826 
00827   map<int, linStruc> linEtaSlice ;
00828   map< vector<int>, linStruc > linMap ;
00829 
00830   // count number of strip per tower
00831   int NbOfStripPerTCC[108][68] ; 
00832   for (int i=0 ; i<108 ; i++)
00833     for (int j=0 ; j<68 ; j++) 
00834       NbOfStripPerTCC[i][j] = 0 ;
00835   const std::vector<DetId> & ebCells = theBarrelGeometry_->getValidDetIds(DetId::Ecal, EcalBarrel);
00836   const std::vector<DetId> & eeCells = theEndcapGeometry_->getValidDetIds(DetId::Ecal, EcalEndcap);
00837   for (vector<DetId>::const_iterator it = ebCells.begin(); it != ebCells.end(); ++it) {
00838     EBDetId id(*it) ;
00839     const EcalTrigTowerDetId towid= id.tower();
00840     const EcalTriggerElectronicsId elId = theMapping_->getTriggerElectronicsId(id) ;
00841     int tccNb = theMapping_->TCCid(towid) ;
00842     int towerInTCC = theMapping_->iTT(towid) ;
00843     int stripInTower = elId.pseudoStripId() ;  
00844     if (stripInTower>NbOfStripPerTCC[tccNb-1][towerInTCC-1]) NbOfStripPerTCC[tccNb-1][towerInTCC-1] = stripInTower ;
00845   }
00846   for (vector<DetId>::const_iterator it = eeCells.begin(); it != eeCells.end(); ++it) {
00847     EEDetId id(*it) ;
00848     const EcalTrigTowerDetId towid= (*eTTmap_).towerOf(id) ;
00849     const EcalTriggerElectronicsId elId = theMapping_->getTriggerElectronicsId(id) ;
00850     int tccNb = theMapping_->TCCid(towid) ;
00851     int towerInTCC = theMapping_->iTT(towid) ; 
00852     int stripInTower = elId.pseudoStripId() ;
00853     if (stripInTower>NbOfStripPerTCC[tccNb-1][towerInTCC-1]) NbOfStripPerTCC[tccNb-1][towerInTCC-1] = stripInTower ;
00854   }
00855 
00856 
00857 
00858   // loop on EB xtals
00859   if (writeToFiles_) (*out_file_)<<"COMMENT ====== barrel crystals ====== "<<std::endl ;
00860 
00861   // special case of eta slices
00862   for (vector<DetId>::const_iterator it = ebCells.begin(); it != ebCells.end(); ++it) {
00863     EBDetId id(*it) ;
00864     double theta = theBarrelGeometry_->getGeometry(id)->getPosition().theta() ;
00865     if (!useTransverseEnergy_) theta = acos(0.) ;
00866     const EcalTrigTowerDetId towid= id.tower();
00867     const EcalTriggerElectronicsId elId = theMapping_->getTriggerElectronicsId(id) ;
00868     int dccNb = theMapping_->DCCid(towid) ;
00869     int tccNb = theMapping_->TCCid(towid) ;
00870     int towerInTCC = theMapping_->iTT(towid) ; // from 1 to 68 (EB)
00871     int stripInTower = elId.pseudoStripId() ;  // from 1 to 5
00872     int xtalInStrip = elId.channelId() ;       // from 1 to 5
00873     const EcalElectronicsId Id = theMapping_->getElectronicsId(id) ;
00874     int CCUid = Id.towerId() ;
00875     int VFEid = Id.stripId() ;
00876     int xtalInVFE = Id.xtalId() ;
00877     int xtalWithinCCUid = 5*(VFEid-1) + xtalInVFE -1 ; // Evgueni expects [0,24]
00878          
00879     (*geomFile_)<<"dccNb = "<<dccNb<<" tccNb = "<<tccNb<<" towerInTCC = "<<towerInTCC
00880                 <<" stripInTower = "<<stripInTower<<" xtalInStrip = "<<xtalInStrip
00881                 <<" CCUid = "<<CCUid<<" VFEid = "<<VFEid<<" xtalInVFE = "<<xtalInVFE
00882                 <<" xtalWithinCCUid = "<<xtalWithinCCUid<<" ieta = "<<id.ieta()<<" iphi = "<<id.iphi()
00883                 <<" xtalhashedId = "<<id.hashedIndex()<<" xtalNb = "<<id.ic()
00884                 <<" ietaTT = "<<towid.ieta()<<" iphiTT = "<<towid.iphi()<<endl ;
00885 
00886     int TCCch = towerInTCC ;
00887     int SLBslot = int((towerInTCC-1)/8.) + 1 ;
00888     int SLBch = (towerInTCC-1)%8 + 1 ;
00889     int cmsswId = id.rawId() ;
00890     int ixtal=(id.ism()-1)*1700+(id.ic()-1);
00891     EcalLogicID logicId =my_EcalLogicId[ixtal];
00892     int dbId = logicId.getLogicID() ;
00893     int val[] = {dccNb+600,tccNb,towerInTCC,stripInTower,
00894                  xtalInStrip,CCUid,VFEid,xtalInVFE,xtalWithinCCUid,id.ieta(),id.iphi(),
00895                  -999,-999,towid.ieta()/abs(towid.ieta()),id.hashedIndex(),
00896                  id.ic(), cmsswId, dbId, 
00897                  towid.ieta(),towid.iphi(), TCCch, getCrate(tccNb).second, SLBch, SLBslot, 
00898                  getGCTRegionEta(towid.ieta()),getGCTRegionPhi(towid.iphi())} ;
00899     for (int i=0 ; i<26 ; i++) ntupleInts_[i] = val[i] ;
00900     
00901     sprintf(ntupleDet_,getDet(tccNb).c_str()) ;
00902     sprintf(ntupleCrate_,getCrate(tccNb).first.c_str()) ;
00903     ntuple->Fill() ;
00904     
00905     
00906     if (tccNb == 37 && stripInTower == 3 && xtalInStrip == 3 && (towerInTCC-1)%4==0) {
00907       int etaSlice = towid.ietaAbs() ;
00908       coeffStruc coeff ;
00909       getCoeff(coeff, calibMap, id.rawId()) ;
00910       getCoeff(coeff, gainMap, id.rawId()) ;
00911       getCoeff(coeff, pedMap, id.rawId()) ;
00912       linStruc lin ;
00913       for (int i=0 ; i<3 ; i++) {
00914         int mult, shift ;
00915         bool ok = computeLinearizerParam(theta, coeff.gainRatio_[i], coeff.calibCoeff_, "EB", mult , shift) ;
00916         if (!ok) std::cout << "unable to compute the parameters for SM="<< id.ism()<<" xt="<< id.ic()<<" " <<dec<<id.rawId()<<std::endl ;  
00917         else {
00918           lin.pedestal_[i] = coeff.pedestals_[i] ;
00919           lin.mult_[i] = mult ;
00920           lin.shift_[i] = shift ;
00921         }
00922       }
00923 
00924       bool ok(true) ;
00925       if (forcedPedestalValue_ == -2) ok = realignBaseline(lin, 0) ;
00926       if (!ok) std::cout<<"SM="<< id.ism()<<" xt="<< id.ic()<<" " <<dec<<id.rawId()<<std::endl ; 
00927       linEtaSlice[etaSlice] = lin ;     
00928     }
00929   }    
00930 
00931   // general case
00932   for (vector<DetId>::const_iterator it = ebCells.begin(); it != ebCells.end(); ++it) {
00933     EBDetId id(*it) ;
00934     double theta = theBarrelGeometry_->getGeometry(id)->getPosition().theta() ;
00935     if (!useTransverseEnergy_) theta = acos(0.) ;
00936     const EcalTrigTowerDetId towid= id.tower();
00937     towerListEB.push_back(towid.rawId()) ;
00938     const EcalTriggerElectronicsId elId = theMapping_->getTriggerElectronicsId(id) ;
00939     stripListEB.push_back(elId.rawId() & 0xfffffff8) ;
00940     int dccNb = theMapping_->DCCid(towid) ;
00941     //int tccNb = theMapping_->TCCid(towid) ;
00942     int towerInTCC = theMapping_->iTT(towid) ; // from 1 to 68 (EB)
00943     //int stripInTower = elId.pseudoStripId() ;  // from 1 to 5
00944     //int xtalInStrip = elId.channelId() ;       // from 1 to 5
00945     const EcalElectronicsId Id = theMapping_->getElectronicsId(id) ;
00946     int CCUid = Id.towerId() ;
00947     int VFEid = Id.stripId() ;
00948     int xtalInVFE = Id.xtalId() ;
00949     int xtalWithinCCUid = 5*(VFEid-1) + xtalInVFE -1 ; // Evgueni expects [0,24]
00950     int etaSlice = towid.ietaAbs() ;
00951 
00952     // hashed index of strip EcalLogicID:
00953     int hashedStripLogicID = 68*5*(id.ism()-1) + 5*(towerInTCC-1) +  (VFEid-1) ;
00954     stripMapEB[hashedStripLogicID] = elId.rawId() & 0xfffffff8 ;
00955     stripMapEBsintheta[elId.rawId() & 0xfffffff8] = SFGVB_Threshold_ + abs(int(sin(theta)*pedestal_offset_)) ;
00956 
00957     //modif-debug
00958     /*
00959     FEConfigFgrEEStripDat stripdebug;
00960     EcalLogicID thestrip_debug = my_StripEcalLogicId[hashedStripLogicID] ;
00961     if(towid.ieta() == 6 && towid.iphi() == 7){
00962       std::cout << "xtal info=" << id << " VFE=" << VFEid << std::endl;
00963       std::cout << "TOWER DEBUG ieta=" << towid.ieta() << " iphi=" << towid.iphi() << " SFGVB=" << SFGVB_Threshold_ + abs(int(sin(theta)*pedestal_offset_)) 
00964                 << " dbId=" << (elId.rawId() & 0xfffffff8) << " " << hashedStripLogicID << " " << thestrip_debug.getLogicID() << std::endl;  //modif-debug
00965     }//EB+3 TT24
00966     */
00967     //std::cout<<std::dec<<SFGVB_Threshold_ + abs(int(sin(theta)*pedestal_offset_))<<" "<<SFGVB_Threshold_<<" "<<abs(int(sin(theta)*pedestal_offset_))<<std::endl ;
00968 
00969     FEConfigPedDat pedDB ;
00970     FEConfigLinDat linDB ;
00971     if (writeToFiles_) (*out_file_)<<"CRYSTAL "<<dec<<id.rawId()<<std::endl ;
00972     //  if (writeToDB_) logicId = db_->getEcalLogicID ("EB_crystal_number", id.ism(), id.ic()) ;
00973 
00974     coeffStruc coeff ;
00975     getCoeff(coeff, calibMap, id.rawId()) ;
00976     if (H2_) coeff.calibCoeff_ = calibvec[id.ic()-1] ;
00977     getCoeff(coeff, gainMap, id.rawId()) ;
00978     getCoeff(coeff, pedMap, id.rawId()) ;
00979     ICEB->Fill(id.iphi(), id.ieta(), coeff.calibCoeff_) ;  
00980     IC->Fill(theBarrelGeometry_->getGeometry(id)->getPosition().phi(), theBarrelGeometry_->getGeometry(id)->getPosition().eta(), coeff.calibCoeff_) ;  
00981 
00982     vector<int> xtalCCU ;
00983     xtalCCU.push_back(dccNb+600) ; 
00984     xtalCCU.push_back(CCUid) ; 
00985     xtalCCU.push_back(xtalWithinCCUid) ;
00986     xtalCCU.push_back(id.rawId()) ;
00987     
00988     // compute and fill linearization parameters
00989     // case of eta slice
00990     if (forceEtaSlice_) {
00991       map<int, linStruc>::const_iterator itLin = linEtaSlice.find(etaSlice);
00992       if (itLin != linEtaSlice.end()) {
00993         linMap[xtalCCU] = itLin->second ;
00994         if (writeToFiles_) {
00995           for (int i=0 ; i<3 ; i++) 
00996             (*out_file_) << hex <<" 0x"<<itLin->second.pedestal_[i]<<" 0x"<<itLin->second.mult_[i]<<" 0x"<<itLin->second.shift_[i]<<std::endl;
00997         }
00998         if (writeToDB_) {
00999           for (int i=0 ; i<3 ; i++) {
01000             if (i==0)  {pedDB.setPedMeanG12(itLin->second.pedestal_[i]) ; linDB.setMultX12(itLin->second.mult_[i]) ; linDB.setShift12(itLin->second.shift_[i]) ; } 
01001             if (i==1)  {pedDB.setPedMeanG6(itLin->second.pedestal_[i]) ; linDB.setMultX6(itLin->second.mult_[i]) ; linDB.setShift6(itLin->second.shift_[i]) ; } 
01002             if (i==2)  {pedDB.setPedMeanG1(itLin->second.pedestal_[i]) ; linDB.setMultX1(itLin->second.mult_[i]) ; linDB.setShift1(itLin->second.shift_[i]) ; } 
01003           }
01004         }
01005         float factor = float(itLin->second.mult_[0])*pow(2.,-itLin->second.shift_[0])/xtal_LSB_EB_ ;
01006         tpgFactorEB->Fill(id.iphi(), id.ieta(), factor) ;
01007         tpgFactor->Fill(theBarrelGeometry_->getGeometry(id)->getPosition().phi(), 
01008                         theBarrelGeometry_->getGeometry(id)->getPosition().eta(), factor) ;
01009       }
01010       else std::cout<<"current EtaSlice = "<<etaSlice<<" not found in the EtaSlice map"<<std::endl ;
01011     }
01012     else {
01013       // general case
01014       linStruc lin ;
01015       int forceBase12 = 0 ;
01016       for (int i=0 ; i<3 ; i++) {
01017         int mult, shift ;
01018         bool ok = computeLinearizerParam(theta, coeff.gainRatio_[i], coeff.calibCoeff_, "EB", mult , shift) ;
01019         if (!ok) std::cout << "unable to compute the parameters for SM="<< id.ism()<<" xt="<< id.ic()<<" " <<dec<<id.rawId()<<std::endl ;  
01020         else {
01021           //PP begin
01022           //  mult = 0 ; shift = 0 ;
01023           //  if (CCUid==1 && xtalWithinCCUid==21) {
01024           //    if (i==0) {mult = 0x80 ; shift = 0x3 ;}
01025           //    if (i==1) {mult = 0x80 ; shift = 0x2 ;}
01026           //    if (i==2) {mult = 0xc0 ; shift = 0x0 ;}
01027           //  } 
01028           //PP end
01029           double base = coeff.pedestals_[i] ;
01030           if (forcedPedestalValue_ == -3 && i==0) {
01031             double G = mult*pow(2.0,-(shift+2)) ;
01032             double g = G/sin(theta) ;
01033             // int pedestal = coeff.pedestals_[i] ;
01034             base = double(coeff.pedestals_[i]) - pedestal_offset_/g ;
01035             if (base<0.) base = 0 ;
01036             forceBase12 = int(base) ;
01037           }
01038           lin.pedestal_[i] = coeff.pedestals_[i] ;
01039           lin.mult_[i] = mult ;
01040           lin.shift_[i] = shift ;
01041 
01042 //        if (xtalWithinCCUid != 14) {
01043 //          forceBase12 = 0 ;
01044 //          lin.pedestal_[i] = 0 ; 
01045 //          lin.mult_[i] = 0 ;
01046 //          lin.shift_[i] = 0 ;
01047 //        }
01048 
01049         }
01050       }
01051 
01052       bool ok(true) ;
01053       if (forcedPedestalValue_ == -2) ok = realignBaseline(lin, 0) ;
01054       if (forcedPedestalValue_ == -3) ok = realignBaseline(lin, forceBase12) ;
01055       if (!ok) std::cout<<"SM="<< id.ism()<<" xt="<< id.ic()<<" " <<dec<<id.rawId()<<std::endl ; 
01056       
01057       for (int i=0 ; i<3 ; i++) {      
01058         if (writeToFiles_) (*out_file_) << hex <<" 0x"<<lin.pedestal_[i]<<" 0x"<<lin.mult_[i]<<" 0x"<<lin.shift_[i]<<std::endl; 
01059         if (writeToDB_) {
01060           if (i==0)  {pedDB.setPedMeanG12(lin.pedestal_[i]) ; linDB.setMultX12(lin.mult_[i]) ; linDB.setShift12(lin.shift_[i]) ; } 
01061           if (i==1)  {pedDB.setPedMeanG6(lin.pedestal_[i]) ; linDB.setMultX6(lin.mult_[i]) ; linDB.setShift6(lin.shift_[i]) ; } 
01062           if (i==2)  {pedDB.setPedMeanG1(lin.pedestal_[i]) ; linDB.setMultX1(lin.mult_[i]) ; linDB.setShift1(lin.shift_[i]) ; }
01063         }
01064         if (i==0) {
01065           float factor = float(lin.mult_[i])*pow(2.,-lin.shift_[i])/xtal_LSB_EB_ ;
01066           tpgFactorEB->Fill(id.iphi(), id.ieta(), factor) ;
01067           tpgFactor->Fill(theBarrelGeometry_->getGeometry(id)->getPosition().phi(), 
01068                           theBarrelGeometry_->getGeometry(id)->getPosition().eta(), factor) ;                       
01069         }
01070         double G = lin.mult_[i]*pow(2.0,-(lin.shift_[i]+2)) ;
01071         double g = G/sin(theta) ;
01072         float val[] = { float(i),
01073                         float(theta),
01074                         float(G),
01075                         float(g),
01076                         float(coeff.pedestals_[i]),
01077                         float(lin.pedestal_[i])};// first arg = gainId (0 means gain12)
01078         ntupleSpike->Fill(val) ;
01079       }
01080       linMap[xtalCCU] = lin ;
01081     }
01082     if (writeToDB_) {
01083     // 1700 crystals/SM in the ECAL barrel
01084       int ixtal=(id.ism()-1)*1700+(id.ic()-1);
01085       EcalLogicID logicId =my_EcalLogicId[ixtal];
01086       pedset[logicId] = pedDB ;
01087       linset[logicId] = linDB ; 
01088     }
01089   } //ebCells
01090 
01091   if (writeToDB_) {
01092     // EcalLogicID  of the whole barrel is: my_EcalLogicId_EB
01093     FEConfigLinParamDat linparam ;
01094     linparam.setETSat(Et_sat_EB_); 
01095     linparamset[my_EcalLogicId_EB] = linparam ;
01096 
01097     FEConfigFgrParamDat fgrparam ;
01098     fgrparam.setFGlowthresh(FG_lowThreshold_EB_); 
01099     fgrparam.setFGhighthresh(FG_highThreshold_EB_); 
01100     fgrparam.setFGlowratio(FG_lowRatio_EB_); 
01101     fgrparam.setFGhighratio(FG_highRatio_EB_);
01102     fgrparamset[my_EcalLogicId_EB] = fgrparam ;
01103 
01104 
01105     FEConfigLUTParamDat lutparam ;
01106     lutparam.setETSat(Et_sat_EB_); 
01107     lutparam.setTTThreshlow(TTF_lowThreshold_EB_); 
01108     lutparam.setTTThreshhigh(TTF_highThreshold_EB_); 
01109     lutparamset[my_EcalLogicId_EB] = lutparam ;
01110 
01111   }
01112 
01113 
01114   // loop on EE xtals
01115   if (writeToFiles_) (*out_file_)<<"COMMENT ====== endcap crystals ====== "<<std::endl ;
01116   
01117   // special case of eta slices
01118   for (vector<DetId>::const_iterator it = eeCells.begin(); it != eeCells.end(); ++it) {
01119     EEDetId id(*it) ;
01120     double theta = theEndcapGeometry_->getGeometry(id)->getPosition().theta() ;
01121     if (!useTransverseEnergy_) theta = acos(0.) ;
01122     const EcalTrigTowerDetId towid= (*eTTmap_).towerOf(id) ;
01123     const EcalTriggerElectronicsId elId = theMapping_->getTriggerElectronicsId(id) ;
01124     const EcalElectronicsId Id = theMapping_->getElectronicsId(id) ;
01125     int dccNb = Id.dccId() ;
01126     int tccNb = theMapping_->TCCid(towid) ;
01127     int towerInTCC = theMapping_->iTT(towid) ; 
01128     int stripInTower = elId.pseudoStripId() ;
01129     int xtalInStrip = elId.channelId() ;
01130     int CCUid = Id.towerId() ;
01131     int VFEid = Id.stripId() ;
01132     int xtalInVFE = Id.xtalId() ;
01133     int xtalWithinCCUid = 5*(VFEid-1) + xtalInVFE -1 ; // Evgueni expects [0,24]
01134 
01135     (*geomFile_)<<"dccNb = "<<dccNb<<" tccNb = "<<tccNb<<" towerInTCC = "<<towerInTCC
01136                 <<" stripInTower = "<<stripInTower<<" xtalInStrip = "<<xtalInStrip
01137                 <<" CCUid = "<<CCUid<<" VFEid = "<<VFEid<<" xtalInVFE = "<<xtalInVFE
01138                 <<" xtalWithinCCUid = "<<xtalWithinCCUid<<" ix = "<<id.ix()<<" iy = "<<id.iy()
01139                 <<" xtalhashedId = "<<id.hashedIndex()<<" xtalNb = "<<id.isc()
01140                 <<" ietaTT = "<<towid.ieta()<<" iphiTT = "<<towid.iphi()<<endl ;
01141 
01142     int TCCch = stripInTower ;
01143     int SLBslot, SLBch ;
01144     if (towerInTCC<5) {
01145       SLBslot = 1 ;
01146       SLBch = 4 + towerInTCC ;
01147     }
01148     else {
01149       SLBslot = int((towerInTCC-5)/8.) + 2 ;
01150       SLBch = (towerInTCC-5)%8 + 1 ;
01151     } 
01152     for (int j=0 ; j<towerInTCC-1 ; j++) TCCch += NbOfStripPerTCC[tccNb-1][j] ;
01153 
01154     int cmsswId = id.rawId() ;
01155     EcalLogicID logicId ;
01156     int iz = id.positiveZ() ;
01157     if (iz ==0) iz = -1 ;
01158     for (int k=0; k<(int)my_CrystalEcalLogicId_EE.size(); k++) {
01159       int z= my_CrystalEcalLogicId_EE[k].getID1() ;
01160       int x= my_CrystalEcalLogicId_EE[k].getID2() ;
01161       int y= my_CrystalEcalLogicId_EE[k].getID3() ;
01162       if (id.ix()==x && id.iy()==y && iz==z) logicId=my_CrystalEcalLogicId_EE[k];
01163     }
01164     int dbId = logicId.getLogicID() ;
01165     
01166     int val[] = {dccNb+600,tccNb,towerInTCC,stripInTower,
01167                  xtalInStrip,CCUid,VFEid,xtalInVFE,xtalWithinCCUid,-999,-999,
01168                  id.ix(),id.iy(),towid.ieta()/abs(towid.ieta()),id.hashedIndex(),
01169                  id.ic(),cmsswId, dbId, 
01170                  towid.ieta(),towid.iphi(),TCCch, getCrate(tccNb).second, SLBch, SLBslot, 
01171                  getGCTRegionEta(towid.ieta()),getGCTRegionPhi(towid.iphi())} ;
01172     for (int i=0 ; i<26 ; i++) ntupleInts_[i] = val[i] ;
01173     sprintf(ntupleDet_,getDet(tccNb).c_str()) ;
01174     sprintf(ntupleCrate_,getCrate(tccNb).first.c_str()) ;
01175     ntuple->Fill() ;
01176      
01177     if ((tccNb == 76 || tccNb == 94) && stripInTower == 1 && xtalInStrip == 3 && (towerInTCC-1)%4==0) {
01178       int etaSlice = towid.ietaAbs() ;
01179       coeffStruc coeff ;
01180       getCoeff(coeff, calibMap, id.rawId()) ;
01181       getCoeff(coeff, gainMap, id.rawId()) ;
01182       getCoeff(coeff, pedMap, id.rawId()) ;
01183       linStruc lin ;
01184       for (int i=0 ; i<3 ; i++) {
01185         int mult, shift ;
01186         bool ok = computeLinearizerParam(theta, coeff.gainRatio_[i], coeff.calibCoeff_, "EE", mult , shift) ;
01187         if (!ok) std::cout << "unable to compute the parameters for Quadrant="<< id.iquadrant()<<" xt="<< id.ic()<<" " <<dec<<id.rawId()<<std::endl ;  
01188         else {
01189           lin.pedestal_[i] = coeff.pedestals_[i] ;
01190           lin.mult_[i] = mult ;
01191           lin.shift_[i] = shift ;
01192         }
01193       }
01194 
01195       bool ok(true) ;
01196       if (forcedPedestalValue_ == -2 || forcedPedestalValue_ == -3) ok = realignBaseline(lin, 0) ;
01197       if (!ok) std::cout<<"Quadrant="<< id.iquadrant()<<" xt="<< id.ic()<<" " <<dec<<id.rawId()<<std::endl ; 
01198 
01199       linEtaSlice[etaSlice] = lin ;
01200     }
01201   }
01202 
01203   // general case
01204   for (vector<DetId>::const_iterator it = eeCells.begin(); it != eeCells.end(); ++it) {
01205     EEDetId id(*it);
01206     double theta = theEndcapGeometry_->getGeometry(id)->getPosition().theta() ;
01207     if (!useTransverseEnergy_) theta = acos(0.) ;
01208     const EcalTrigTowerDetId towid= (*eTTmap_).towerOf(id) ;
01209     const EcalTriggerElectronicsId elId = theMapping_->getTriggerElectronicsId(id) ;
01210     const EcalElectronicsId Id = theMapping_->getElectronicsId(id) ;
01211     towerListEE.push_back(towid.rawId()) ;
01212     // special case of towers in inner rings of EE
01213     if (towid.ietaAbs() == 27 || towid.ietaAbs() == 28) {
01214       EcalTrigTowerDetId additionalTower(towid.zside(), towid.subDet(), towid.ietaAbs(), towid.iphi()+1) ;
01215       towerListEE.push_back(additionalTower.rawId()) ;
01216     }
01217     stripListEE.push_back(elId.rawId() & 0xfffffff8) ;
01218     int dccNb = Id.dccId() ;
01219     //int tccNb = theMapping_->TCCid(towid) ;
01220     //int towerInTCC = theMapping_->iTT(towid) ;
01221     //int stripInTower = elId.pseudoStripId() ;
01222     //int xtalInStrip = elId.channelId() ;
01223     int CCUid = Id.towerId() ;
01224     int VFEid = Id.stripId() ;
01225     int xtalInVFE = Id.xtalId() ;
01226     int xtalWithinCCUid = 5*(VFEid-1) + xtalInVFE -1 ; // Evgueni expects [0,24]
01227     int etaSlice = towid.ietaAbs() ;
01228 
01229     EcalLogicID logicId ;
01230     FEConfigPedDat pedDB ;
01231     FEConfigLinDat linDB ;
01232     if (writeToFiles_) (*out_file_)<<"CRYSTAL "<<dec<<id.rawId()<<std::endl ;
01233     if (writeToDB_ && DBEE_) {
01234       int iz = id.positiveZ() ;
01235       if (iz ==0) iz = -1 ;      
01236       for (int k=0; k<(int)my_CrystalEcalLogicId_EE.size(); k++) {
01237         int z= my_CrystalEcalLogicId_EE[k].getID1() ;
01238         int x= my_CrystalEcalLogicId_EE[k].getID2() ;
01239         int y= my_CrystalEcalLogicId_EE[k].getID3() ;   
01240         if (id.ix()==x && id.iy()==y && iz==z) logicId=my_CrystalEcalLogicId_EE[k];
01241       }
01242     }
01243     
01244     coeffStruc coeff ;
01245     getCoeff(coeff, calibMap, id.rawId()) ;
01246     getCoeff(coeff, gainMap, id.rawId()) ;
01247     getCoeff(coeff, pedMap, id.rawId()) ;
01248     if (id.zside()>0) ICEEPlus->Fill(id.ix(), id.iy(), coeff.calibCoeff_) ;  
01249     else ICEEMinus->Fill(id.ix(), id.iy(), coeff.calibCoeff_) ;  
01250     IC->Fill(theEndcapGeometry_->getGeometry(id)->getPosition().phi(), theEndcapGeometry_->getGeometry(id)->getPosition().eta(), coeff.calibCoeff_) ;  
01251   
01252     vector<int> xtalCCU ;
01253     xtalCCU.push_back(dccNb+600) ; 
01254     xtalCCU.push_back(CCUid) ; 
01255     xtalCCU.push_back(xtalWithinCCUid) ;
01256     xtalCCU.push_back(id.rawId()) ;
01257 
01258     // compute and fill linearization parameters
01259     // case of eta slice
01260     if (forceEtaSlice_) {
01261       map<int, linStruc>::const_iterator itLin = linEtaSlice.find(etaSlice);
01262       if (itLin != linEtaSlice.end()) {
01263         linMap[xtalCCU] = itLin->second ;
01264         if (writeToFiles_) {
01265           for (int i=0 ; i<3 ; i++) 
01266             (*out_file_) << hex <<" 0x"<<itLin->second.pedestal_[i]<<" 0x"<<itLin->second.mult_[i]<<" 0x"<<itLin->second.shift_[i]<<std::endl;
01267         }
01268         if (writeToDB_ && DBEE_) {
01269           for (int i=0 ; i<3 ; i++) {
01270             if (i==0)  {pedDB.setPedMeanG12(itLin->second.pedestal_[i]) ; linDB.setMultX12(itLin->second.mult_[i]) ; linDB.setShift12(itLin->second.shift_[i]) ; } 
01271             if (i==1)  {pedDB.setPedMeanG6(itLin->second.pedestal_[i]) ; linDB.setMultX6(itLin->second.mult_[i]) ; linDB.setShift6(itLin->second.shift_[i]) ; } 
01272             if (i==2)  {pedDB.setPedMeanG1(itLin->second.pedestal_[i]) ; linDB.setMultX1(itLin->second.mult_[i]) ; linDB.setShift1(itLin->second.shift_[i]) ; } 
01273           }
01274         }
01275         float factor = float(itLin->second.mult_[0])*pow(2.,-itLin->second.shift_[0])/xtal_LSB_EE_ ;
01276         if (id.zside()>0) tpgFactorEEPlus->Fill(id.ix(), id.iy(), factor) ;
01277         else tpgFactorEEMinus->Fill(id.ix(), id.iy(), factor) ;
01278         tpgFactor->Fill(theEndcapGeometry_->getGeometry(id)->getPosition().phi(), 
01279                         theEndcapGeometry_->getGeometry(id)->getPosition().eta(), factor) ;
01280       }
01281       else std::cout<<"current EtaSlice = "<<etaSlice<<" not found in the EtaSlice map"<<std::endl ;      
01282     }
01283     else {
01284       // general case
01285       linStruc lin ;
01286       for (int i=0 ; i<3 ; i++) {
01287         int mult, shift ;
01288         bool ok = computeLinearizerParam(theta, coeff.gainRatio_[i], coeff.calibCoeff_, "EE", mult , shift) ;
01289         if (!ok) std::cout << "unable to compute the parameters for "<<dec<<id.rawId()<<std::endl ;  
01290         else {
01291           lin.pedestal_[i] = coeff.pedestals_[i] ;
01292           lin.mult_[i] = mult ;
01293           lin.shift_[i] = shift ;
01294         }
01295       }
01296 
01297       bool ok(true) ;
01298       if (forcedPedestalValue_ == -2 || forcedPedestalValue_ == -3) ok = realignBaseline(lin, 0) ;
01299       if (!ok) std::cout<<"Quadrant="<< id.iquadrant()<<" xt="<< id.ic()<<" " <<dec<<id.rawId()<<std::endl ; 
01300 
01301       for (int i=0 ; i<3 ; i++) {
01302         if (writeToFiles_) (*out_file_) << hex <<" 0x"<<lin.pedestal_[i]<<" 0x"<<lin.mult_[i]<<" 0x"<<lin.shift_[i]<<std::endl; 
01303         if (writeToDB_ && DBEE_) {
01304           if (i==0)  {pedDB.setPedMeanG12(lin.pedestal_[i]) ; linDB.setMultX12(lin.mult_[i]) ; linDB.setShift12(lin.shift_[i]) ; } 
01305           if (i==1)  {pedDB.setPedMeanG6(lin.pedestal_[i]) ; linDB.setMultX6(lin.mult_[i]) ; linDB.setShift6(lin.shift_[i]) ; } 
01306           if (i==2)  {pedDB.setPedMeanG1(lin.pedestal_[i]) ; linDB.setMultX1(lin.mult_[i]) ; linDB.setShift1(lin.shift_[i]) ; }
01307         }
01308         if (i==0) {
01309           float factor = float(lin.mult_[i])*pow(2.,-lin.shift_[i])/xtal_LSB_EE_ ;
01310           if (id.zside()>0) tpgFactorEEPlus->Fill(id.ix(), id.iy(), factor) ;
01311           else tpgFactorEEMinus->Fill(id.ix(), id.iy(), factor) ;           
01312           tpgFactor->Fill(theEndcapGeometry_->getGeometry(id)->getPosition().phi(), 
01313                           theEndcapGeometry_->getGeometry(id)->getPosition().eta(), factor) ;                               
01314         }
01315       }
01316       linMap[xtalCCU] = lin ;
01317     }
01318     if (writeToDB_ && DBEE_) {
01319       pedset[logicId] = pedDB ;
01320       linset[logicId] = linDB ;
01321     }
01322   } //eeCells
01323 
01324   if (writeToDB_ ) {
01325     // EcalLogicID  of the whole barrel is: my_EcalLogicId_EB
01326     FEConfigLinParamDat linparam ;
01327     linparam.setETSat(Et_sat_EE_);
01328     linparamset[my_EcalLogicId_EE] = linparam ;
01329 
01330     FEConfigLUTParamDat lutparam ;
01331     lutparam.setETSat(Et_sat_EE_);
01332     lutparam.setTTThreshlow(TTF_lowThreshold_EE_); 
01333     lutparam.setTTThreshhigh(TTF_highThreshold_EE_);
01334     lutparamset[my_EcalLogicId_EE] = lutparam ;
01335 
01336  
01337     FEConfigFgrParamDat fgrparam ;
01338     fgrparam.setFGlowthresh(FG_Threshold_EE_); 
01339     fgrparam.setFGhighthresh(FG_Threshold_EE_); 
01340     fgrparamset[my_EcalLogicId_EE] = fgrparam ;
01341 
01342   }
01343 
01344   if (writeToDB_) {
01345 
01346     ostringstream ltag;
01347     ltag.str("EB_"); ltag<<Et_sat_EB_<<"_EE_"<<Et_sat_EE_;
01348     std::string lin_tag=ltag.str();
01349     std::cout<< " LIN tag "<<lin_tag<<endl;
01350 
01351     if(m_write_ped==1) {
01352       ped_conf_id_=db_->writeToConfDB_TPGPedestals(pedset, 1, "from_OfflineDB") ;
01353     } else {
01354       std::cout<< "the ped id ="<<ped_conf_id_<<" will be used for the pedestals "<<std::endl; 
01355     }
01356 
01357     if(m_write_lin==1) lin_conf_id_=db_->writeToConfDB_TPGLinearCoef(linset,linparamset, 1, lin_tag) ;
01358   }
01359 
01361   // Evgueni interface
01363   std::ofstream evgueni("TPG_hardcoded.hh", std::ios::out) ;  
01364   evgueni<<"void getLinParamTPG_hardcoded(int fed, int ccu, int xtal,"<<endl ;
01365   evgueni<<"                        int & mult12, int & shift12, int & base12,"<<endl ;
01366   evgueni<<"                        int & mult6, int & shift6, int & base6,"<<endl ;
01367   evgueni<<"                        int & mult1, int & shift1, int & base1)"<<endl ;
01368   evgueni<<"{"<<endl;
01369   evgueni<<"  mult12 = 0 ; shift12 = 0 ; base12 = 0 ; mult6 = 0 ; shift6 = 0 ; base6 = 0 ; mult1 = 0 ; shift1 = 0 ; base1 = 0 ;"<<endl ;
01370   map< vector<int>, linStruc>::const_iterator itLinMap ;
01371   for (itLinMap = linMap.begin() ; itLinMap != linMap.end() ; itLinMap++) {
01372     vector<int> xtalInCCU = itLinMap->first ;
01373     evgueni<<"  if (fed=="<<xtalInCCU[0]<<" && ccu=="<<xtalInCCU[1]<<" && xtal=="<<xtalInCCU[2]<<") {" ;
01374     evgueni<<"  mult12 = "<<itLinMap->second.mult_[0]<<" ; shift12 = "<<itLinMap->second.shift_[0]<<" ; base12 = "<<itLinMap->second.pedestal_[0]<<" ; " ;
01375     evgueni<<"  mult6 = "<<itLinMap->second.mult_[1]<<" ; shift6 = "<<itLinMap->second.shift_[1]<<" ; base6 = "<<itLinMap->second.pedestal_[1]<<" ; " ;
01376     evgueni<<"  mult1 = "<<itLinMap->second.mult_[2]<<" ; shift1 = "<<itLinMap->second.shift_[2]<<" ; base1 = "<<itLinMap->second.pedestal_[2]<<" ; " ;
01377     evgueni<<"  return ;}" <<endl ;
01378   }
01379   evgueni<<"}" <<endl ;
01380   evgueni.close() ;
01381 
01382 
01383 
01385   // Compute weights section //
01387 
01388   const int NWEIGROUPS = 2 ; 
01389   std::vector<unsigned int> weights[NWEIGROUPS] ;
01390 
01391 #if (CMSSW_VERSION>=340)
01392   EBShape shapeEB ;
01393   EEShape shapeEE ;
01394   weights[0] = computeWeights(shapeEB, hshapeEB) ;
01395   weights[1] = computeWeights(shapeEE, hshapeEE) ;
01396 #else
01397   // loading reference signal representation
01398   EcalSimParameterMap parameterMap;  
01399   EBDetId   barrel(1,1);
01400   double    phase = parameterMap.simParameters(barrel).timePhase();
01401   EcalShape shape(phase); 
01402   weights[0] = computeWeights(shape, hshapeEB) ;
01403   weights[1] = weights[0] ;
01404 #endif
01405 
01406   map<EcalLogicID, FEConfigWeightGroupDat> dataset;
01407 
01408   for (int igrp=0 ; igrp<NWEIGROUPS ; igrp++) {
01409 
01410     if (weights[igrp].size() == 5) {
01411 
01412       if (writeToFiles_) {
01413         (*out_file_) <<std::endl ;
01414         (*out_file_) <<"WEIGHT "<<igrp<<endl ;
01415         for (unsigned int sample=0 ; sample<5 ; sample++) (*out_file_) << "0x" <<hex<<weights[igrp][sample]<<" " ;
01416         (*out_file_)<<std::endl ; 
01417         (*out_file_) <<std::endl ;
01418       }
01419       if (writeToDB_) {
01420         std::cout<<"going to write the weights for groupe:"<<igrp<<endl;
01421         FEConfigWeightGroupDat gut;
01422         gut.setWeightGroupId(igrp);
01423         //PP WARNING: weights order is reverted when stored in the DB 
01424         gut.setWeight0(weights[igrp][4] );
01425         gut.setWeight1(weights[igrp][3]+ 0x80); //0x80 to identify the max of the pulse in the FENIX (doesn't exist in emulator)
01426         gut.setWeight2(weights[igrp][2] );
01427         gut.setWeight3(weights[igrp][1] );
01428         gut.setWeight4(weights[igrp][0] );
01429         EcalLogicID ecid = EcalLogicID( "DUMMY", igrp,igrp); //1 dummy ID per group
01430         // Fill the dataset
01431         dataset[ecid] = gut;
01432       }
01433     }
01434   }
01435 
01436   if (writeToDB_) {
01437 
01438     // now we store in the DB the correspondence btw channels and groups 
01439     map<EcalLogicID, FEConfigWeightDat> dataset2;
01440 
01441     // EB loop
01442     for (int ich=0; ich<(int)my_StripEcalLogicId.size() ; ich++){
01443       FEConfigWeightDat wut;
01444       int igroup = 0; // this group is for EB
01445       wut.setWeightGroupId(igroup);
01446       dataset2[my_StripEcalLogicId[ich]] = wut;
01447     }
01448         
01449     // EE loop
01450     for (int ich=0; ich<(int)my_StripEcalLogicId1_EE.size() ; ich++){
01451       FEConfigWeightDat wut;
01452       int igroup = 1; // this group is for EE
01453       wut.setWeightGroupId(igroup);
01454       // Fill the dataset
01455       dataset2[my_StripEcalLogicId1_EE[ich]] = wut;
01456     }
01457     // EE loop 2 (we had to split the ids of EE in 2 vectors to avoid crash!)
01458     for (int ich=0; ich<(int)my_StripEcalLogicId2_EE.size() ; ich++){
01459       FEConfigWeightDat wut;
01460       int igroup = 1; // this group is for EE
01461       wut.setWeightGroupId(igroup);
01462       // Fill the dataset
01463       dataset2[my_StripEcalLogicId2_EE[ich]] = wut;
01464     }
01465 
01466     // Insert the datasets
01467     ostringstream wtag;
01468     wtag.str(""); wtag<<"Shape_NGroups_"<<NWEIGROUPS;
01469     std::string weight_tag=wtag.str();
01470     std::cout<< " weight tag "<<weight_tag<<endl; 
01471     if (m_write_wei==1) wei_conf_id_=db_->writeToConfDB_TPGWeight(dataset, dataset2, NWEIGROUPS, weight_tag) ;
01472   }      
01473 
01474 
01476   // Compute FG section //
01478 
01479   // barrel
01480   unsigned int lowRatio, highRatio, lowThreshold, highThreshold, lutFG ;
01481   computeFineGrainEBParameters(lowRatio, highRatio, lowThreshold, highThreshold, lutFG) ;
01482   if (writeToFiles_) {
01483     (*out_file_) <<std::endl ;
01484     (*out_file_) <<"FG 0"<<std::endl ;
01485     (*out_file_)<<hex<<"0x"<<lowThreshold<<" 0x"<<highThreshold
01486                   <<" 0x"<<lowRatio<<" 0x"<<highRatio<<" 0x"<<lutFG
01487                   <<std::endl ;
01488   }
01489 
01490   // endcap
01491   unsigned int threshold, lut_tower ;
01492   unsigned int lut_strip;
01493   computeFineGrainEEParameters(threshold, lut_strip, lut_tower) ; 
01494 
01495   // and here we store the fgr part
01496 
01497   
01498   if (writeToDB_) {
01499     std::cout<<"going to write the fgr "<< endl;
01500       map<EcalLogicID, FEConfigFgrGroupDat> dataset;
01501       // we create 1 group
01502       int NFGRGROUPS =1; 
01503       for (int ich=0; ich<NFGRGROUPS; ich++){
01504         FEConfigFgrGroupDat gut;
01505         gut.setFgrGroupId(ich);
01506         gut.setThreshLow(lowRatio );
01507         gut.setThreshHigh(highRatio);
01508         gut.setRatioLow(lowThreshold);
01509         gut.setRatioHigh(highThreshold);
01510         gut.setLUTValue(lutFG);
01511         EcalLogicID ecid = EcalLogicID( "DUMMY", ich,ich);
01512         // Fill the dataset
01513         dataset[ecid] = gut; // we use any logic id but different, because it is in any case ignored... 
01514       }
01515       
01516       // now we store in the DB the correspondence btw channels and groups 
01517       map<EcalLogicID, FEConfigFgrDat> dataset2;
01518       // in this case I decide in a stupid way which channel belongs to which group 
01519       for (int ich=0; ich<(int)my_TTEcalLogicId.size() ; ich++){
01520         FEConfigFgrDat wut;
01521         int igroup=0;
01522         wut.setFgrGroupId(igroup);
01523         // Fill the dataset
01524         // the logic ids are ordered by SM (1,...36) and TT (1,...68)  
01525         // you have to calculate the right index here 
01526         dataset2[my_TTEcalLogicId[ich]] = wut;
01527       }
01528 
01529       // endcap loop
01530       for (int ich=0; ich<(int)my_RTEcalLogicId_EE.size() ; ich++){
01531         //      std::cout << " endcap FGR " << std::endl;
01532         FEConfigFgrDat wut;
01533         int igroup=0;
01534         wut.setFgrGroupId(igroup);
01535         // Fill the dataset
01536         // the logic ids are ordered by .... ?  
01537         // you have to calculate the right index here 
01538         dataset2[my_RTEcalLogicId_EE[ich]] = wut;
01539       }
01540 
01541       // endcap TT loop for the FEfgr EE Tower
01542       map<EcalLogicID, FEConfigFgrEETowerDat> dataset3;
01543       for (int ich=0; ich<(int)my_TTEcalLogicId_EE.size() ; ich++){
01544         FEConfigFgrEETowerDat fgreett;
01545         fgreett.setLutValue(lut_tower);
01546         dataset3[my_TTEcalLogicId_EE[ich]]=fgreett;
01547       }
01548 
01549       // endcap strip loop for the FEfgr EE strip
01550       // and barrel strip loop for the spike parameters (same structure than EE FGr)
01551       map<EcalLogicID, FEConfigFgrEEStripDat> dataset4;
01552       for (int ich=0; ich<(int)my_StripEcalLogicId1_EE.size() ; ich++){
01553         FEConfigFgrEEStripDat zut;
01554         zut.setThreshold(threshold);
01555         zut.setLutFgr(lut_strip);
01556         dataset4[my_StripEcalLogicId1_EE[ich]] = zut;
01557       }
01558       for (int ich=0; ich<(int)my_StripEcalLogicId2_EE.size() ; ich++){
01559         FEConfigFgrEEStripDat zut;
01560         zut.setThreshold(threshold);
01561         zut.setLutFgr(lut_strip);
01562         // Fill the dataset
01563         dataset4[my_StripEcalLogicId2_EE[ich]] = zut;
01564       }
01565       for (int ich=0; ich<(int)my_StripEcalLogicId.size() ; ich++){
01566         // EB
01567         FEConfigFgrEEStripDat zut;
01568         EcalLogicID thestrip = my_StripEcalLogicId[ich] ;
01569         uint32_t elStripId = stripMapEB[ich] ;
01570         map<uint32_t, uint32_t>::const_iterator it = stripMapEBsintheta.find(elStripId) ;
01571         if (it != stripMapEBsintheta.end()) zut.setThreshold(it->second);
01572         else {
01573           cout<<"ERROR: strip SFGVB threshold parameter not found for that strip:"<<thestrip.getID1()<<" "<<thestrip.getID3()<<" "<<thestrip.getID3()<<endl ; 
01574           cout<<" using value = "<<SFGVB_Threshold_+pedestal_offset_<<endl ;
01575           zut.setThreshold(SFGVB_Threshold_+pedestal_offset_) ;
01576         }
01577         zut.setLutFgr(SFGVB_lut_);
01578         // Fill the dataset
01579         dataset4[thestrip] = zut;
01580       }
01581 
01582       // Insert the dataset
01583       ostringstream wtag;
01584       wtag.str(""); wtag<<"FGR_"<<lutFG<<"_N_"<<NFGRGROUPS<<"_eb_"<<FG_lowThreshold_EB_<<"_EB_"<<FG_highThreshold_EB_;
01585       std::string weight_tag=wtag.str();
01586       std::cout<< " weight tag "<<weight_tag<<endl; 
01587       if(m_write_fgr==1) fgr_conf_id_=db_->writeToConfDB_TPGFgr(dataset, dataset2,  fgrparamset,dataset3, dataset4, NFGRGROUPS, weight_tag) ;
01588 
01589       //modif-alex 21/01/11
01590       map<EcalLogicID, FEConfigSpikeDat> datasetspike; //loob EB TT
01591       for (int ich=0; ich<(int)my_TTEcalLogicId.size() ; ich++){
01592         FEConfigSpikeDat spiketh;
01593         spiketh.setSpikeThreshold(SFGVB_SpikeKillingThreshold_);
01594         datasetspike[my_TTEcalLogicId[ich]] = spiketh;
01595       }//loop EB TT towers    
01596 
01597       //modif-alex 21/01/11
01598       ostringstream stag;
01599       stag.str(""); stag<<"SpikeTh"<<SFGVB_SpikeKillingThreshold_;
01600       std::string spike_tag=stag.str();
01601       std::cout<< " spike tag "<<spike_tag<<endl;
01602       if(m_write_spi==1) spi_conf_id_=db_->writeToConfDB_Spike(datasetspike, spike_tag) ; //modif-alex 21/01/11
01603 
01604       //modif-alex 31/01/11
01605       //DELAYS EB
01606       map<EcalLogicID, FEConfigTimingDat> datasetdelay; // the loop goes from TCC 38 to 72 and throught the towers from 1 to 68 
01607       for (int ich=0; ich<(int)my_TTEcalLogicId_EB_by_TCC.size() ; ich++){
01608         FEConfigTimingDat delay;
01609 
01610         EcalLogicID logiciddelay = my_TTEcalLogicId_EB_by_TCC[ich] ;
01611         int id1_tcc=my_TTEcalLogicId_EB_by_TCC[ich].getID1(); // the TCC
01612         int id2_tt =my_TTEcalLogicId_EB_by_TCC[ich].getID2(); // the tower
01613         std::map<int, vector<int> >::const_iterator ittEB  = delays_EB_.find(id1_tcc); 
01614         std::vector<int> TimingDelaysEB = ittEB->second;
01615 
01616         if (ittEB != delays_EB_.end()){
01617           if(TimingDelaysEB[id2_tt-1] == -1){
01618               std::cout << "ERROR: Barrel timing delay not specified, check file, putting default value 1" << std::endl;
01619               delay.setTimingPar1(1); 
01620             }
01621           else delay.setTimingPar1(TimingDelaysEB[id2_tt-1]);
01622         } else {
01623           std::cout << "ERROR:Barrel Could not find delay parameter for that trigger tower " << std::endl;
01624           std::cout << "Using default value = 1" << std::endl;
01625           delay.setTimingPar1(1);
01626         }
01627 
01628         std::map<int, vector<int> >::const_iterator ittpEB  = phases_EB_.find(id1_tcc);
01629         std::vector<int> TimingPhasesEB = ittpEB->second;
01630 
01631         if (ittpEB != phases_EB_.end()){
01632           if(TimingPhasesEB[id2_tt-1] == -1){
01633             std::cout << "ERROR: Barrel timing phase not specified, check file, putting default value 0" << std::endl;
01634             delay.setTimingPar2(0);
01635           }
01636           else delay.setTimingPar2(TimingPhasesEB[id2_tt-1]);
01637         } else {
01638           std::cout << "ERROR:Barrel Could not find phase parameter for that trigger tower " << std::endl;
01639           std::cout << "Using default value = 0" << std::endl;
01640           delay.setTimingPar2(0);
01641         }
01642 
01643         std::cout << ich << " tcc=" << id1_tcc << " TT=" << id2_tt << " logicId=" << logiciddelay.getLogicID() 
01644                   << " delay=" << TimingDelaysEB[id2_tt-1] << " phase=" << TimingPhasesEB[id2_tt-1] <<  std::endl; 
01645         
01646         //delay.setTimingPar1(1);
01647         //delay.setTimingPar2(2);
01648         datasetdelay[my_TTEcalLogicId_EB_by_TCC[ich]] = delay;
01649       }//loop EB TT towers
01650 
01651       //DELAYS EE
01652       int stripindex = 0;
01653       int tccin = 1;
01654       for (int ich=0; ich<(int)my_StripEcalLogicId_EE_strips_by_TCC.size() ; ich++){
01655         FEConfigTimingDat delay;
01656         //int id1_strip=my_StripEcalLogicId_EE_strips_by_TCC[ich].getID1(); // the TCC
01657         //int id2_strip=my_StripEcalLogicId_EE_strips_by_TCC[ich].getID2(); // the Tower
01658         //int id3_strip=my_StripEcalLogicId_EE_strips_by_TCC[ich].getID3(); // the strip
01659 
01660         EcalLogicID logiciddelay = my_StripEcalLogicId_EE_strips_by_TCC[ich] ;
01661         int id1_tcc=my_StripEcalLogicId_EE_strips_by_TCC[ich].getID1(); // the TCC
01662         int id2_tt =my_StripEcalLogicId_EE_strips_by_TCC[ich].getID2(); // the tower
01663         int id3_st =my_StripEcalLogicId_EE_strips_by_TCC[ich].getID3(); // the strip
01664 
01665         //reset strip counter
01666         if(id1_tcc != tccin) {tccin = id1_tcc; stripindex=0;}
01667 
01668         std::map<int, vector<int> >::const_iterator ittEE  = delays_EE_.find(id1_tcc);
01669         std::vector<int> TimingDelaysEE = ittEE->second;
01670 
01671         if (ittEE != delays_EE_.end()){
01672           if(TimingDelaysEE[stripindex] == -1){
01673             std::cout << "ERROR: Endcap timing delay not specified, check file, putting default value 1" << std::endl;
01674             delay.setTimingPar1(1);
01675           }
01676           else delay.setTimingPar1(TimingDelaysEE[stripindex]);
01677         } else {
01678           std::cout << "ERROR:Endcap Could not find delay parameter for that trigger tower " << std::endl;
01679           std::cout << "Using default value = 1" << std::endl;
01680           delay.setTimingPar1(1);
01681         }
01682 
01683         std::map<int, vector<int> >::const_iterator ittpEE  = phases_EE_.find(id1_tcc);
01684         std::vector<int> TimingPhasesEE = ittpEE->second;
01685 
01686         if (ittpEE != phases_EE_.end()){
01687           if(TimingPhasesEE[stripindex] == -1){
01688             std::cout << "ERROR: Endcap timing phase not specified, check file, putting default value 0" << std::endl;
01689             delay.setTimingPar2(0);
01690           }
01691           else delay.setTimingPar2(TimingPhasesEE[stripindex]);
01692         } else {
01693           std::cout << "ERROR:Endcap Could not find phase parameter for that trigger tower " << std::endl;
01694           std::cout << "Using default value = 0" << std::endl;
01695           delay.setTimingPar2(0);
01696         }
01697 
01698         std::cout << ich << " stripindex=" << stripindex << " tcc=" << id1_tcc << " TT=" << id2_tt << " id3_st=" << id3_st 
01699                   << " logicId=" << logiciddelay.getLogicID()
01700                   << " delay=" << TimingDelaysEE[stripindex] << " phase=" << TimingPhasesEE[stripindex] <<  std::endl;
01701 
01702 
01703         //delay.setTimingPar1(1);
01704         //delay.setTimingPar2(2);
01705         datasetdelay[my_StripEcalLogicId_EE_strips_by_TCC[ich]] = delay;
01706         stripindex++;
01707       }//loop EE strip towers
01708 
01709       ostringstream de_tag;
01710       de_tag.str(""); de_tag<<"DelaysFromFile";
01711       std::string delay_tag=de_tag.str();
01712       std::cout<< " delay tag "<<delay_tag<<endl;
01713       if(m_write_del==1) del_conf_id_=db_->writeToConfDB_Delay(datasetdelay, delay_tag) ; //modif-alex 31/01/11   
01714 
01715 
01716   } //write to DB
01717 
01718   if (writeToDB_) {
01719     std::cout<<"going to write the sliding "<< endl;
01720       map<EcalLogicID, FEConfigSlidingDat> dataset;
01721       // in this case I decide in a stupid way which channel belongs to which group 
01722       for (int ich=0; ich<(int)my_StripEcalLogicId.size() ; ich++){
01723         FEConfigSlidingDat wut;
01724         wut.setSliding(sliding_);
01725         // Fill the dataset
01726         // the logic ids are ordered by SM (1,...36) , TT (1,...68) and strip (1..5) 
01727         // you have to calculate the right index here 
01728         dataset[my_StripEcalLogicId[ich]] = wut;
01729       }
01730 
01731       // endcap loop
01732       for (int ich=0; ich<(int)my_StripEcalLogicId1_EE.size() ; ich++){
01733         FEConfigSlidingDat wut;
01734         wut.setSliding(sliding_);
01735         // Fill the dataset
01736         // the logic ids are ordered by fed tower strip
01737         // you have to calculate the right index here 
01738         dataset[my_StripEcalLogicId1_EE[ich]] = wut;
01739       }
01740       for (int ich=0; ich<(int)my_StripEcalLogicId2_EE.size() ; ich++){
01741         FEConfigSlidingDat wut;
01742         wut.setSliding(sliding_);
01743         // Fill the dataset
01744         // the logic ids are ordered by ... ? 
01745         // you have to calculate the right index here 
01746         dataset[my_StripEcalLogicId2_EE[ich]] = wut;
01747       }
01748       
01749 
01750       // Insert the dataset
01751       ostringstream wtag;
01752       wtag.str(""); wtag<<"Sliding_"<<sliding_;
01753       std::string justatag=wtag.str();
01754       std::cout<< " sliding tag "<<justatag<<endl;
01755       int iov_id=0; // just a parameter ... 
01756       if(m_write_sli==1) sli_conf_id_=db_->writeToConfDB_TPGSliding(dataset,iov_id, justatag) ;
01757   }
01758 
01759   
01760 
01761 
01762 
01764   // Compute LUT section //
01766 
01767   int lut_EB[1024], lut_EE[1024] ;
01768 
01769   // barrel
01770   computeLUT(lut_EB, "EB") ; 
01771   if (writeToFiles_) {
01772     (*out_file_) <<std::endl ;
01773     (*out_file_) <<"LUT 0"<<std::endl ;
01774     for (int i=0 ; i<1024 ; i++) (*out_file_)<<"0x"<<hex<<lut_EB[i]<<endl ;
01775     (*out_file_)<<endl ;
01776   }
01777   
01778   // endcap
01779   computeLUT(lut_EE, "EE") ;
01780   // check first if lut_EB and lut_EE are the same
01781   bool newLUT(false) ;
01782   for (int i=0 ; i<1024 ; i++) if (lut_EE[i] != lut_EB[i]) newLUT = true ;
01783   if (newLUT && writeToFiles_) { 
01784     (*out_file_) <<std::endl ;
01785     (*out_file_) <<"LUT 1"<<std::endl ;
01786     for (int i=0 ; i<1024 ; i++) (*out_file_)<<"0x"<<hex<<lut_EE[i]<<endl ;
01787     (*out_file_)<<endl ;
01788   }
01789 
01790 
01791 
01792   if (writeToDB_) {
01793 
01794 
01795 
01796     map<EcalLogicID, FEConfigLUTGroupDat> dataset;
01797     // we create 1 LUT group
01798     int NLUTGROUPS =0; 
01799     int ich=0;
01800       FEConfigLUTGroupDat lut;
01801       lut.setLUTGroupId(ich);
01802       for (int i=0; i<1024; i++){
01803         lut.setLUTValue(i, lut_EB[i] );
01804       }
01805       EcalLogicID ecid = EcalLogicID( "DUMMY", ich,ich);
01806       // Fill the dataset
01807       dataset[ecid] = lut; // we use any logic id but different, because it is in any case ignored... 
01808       
01809       ich++;
01810       
01811       FEConfigLUTGroupDat lute;
01812       lute.setLUTGroupId(ich);
01813       for (int i=0; i<1024; i++){
01814         lute.setLUTValue(i, lut_EE[i] );
01815       }
01816       EcalLogicID ecide = EcalLogicID( "DUMMY", ich,ich);
01817       // Fill the dataset
01818       dataset[ecide] = lute; // we use any logic id but different, because it is in any case ignored...
01819 
01820       ich++;
01821 
01822       NLUTGROUPS=ich;
01823 
01824     // now we store in the DB the correspondence btw channels and LUT groups 
01825     map<EcalLogicID, FEConfigLUTDat> dataset2;
01826     // in this case I decide in a stupid way which channel belongs to which group 
01827     for (int ich=0; ich<(int)my_TTEcalLogicId.size() ; ich++){
01828       FEConfigLUTDat lut;
01829       int igroup=0;
01830       lut.setLUTGroupId(igroup);
01831       // calculate the right TT - in the vector they are ordered by SM and by TT 
01832       // Fill the dataset
01833       dataset2[my_TTEcalLogicId[ich]] = lut;
01834     }
01835 
01836     // endcap loop 
01837     for (int ich=0; ich<(int)my_TTEcalLogicId_EE.size() ; ich++){
01838       FEConfigLUTDat lut;
01839       int igroup=1;
01840       lut.setLUTGroupId(igroup);
01841       // calculate the right TT  
01842       // Fill the dataset
01843       dataset2[my_TTEcalLogicId_EE[ich]] = lut;
01844     }    
01845 
01846     // Insert the dataset
01847     ostringstream ltag;
01848     ltag.str(""); ltag<<LUT_option_<<"_NGroups_"<<NLUTGROUPS;
01849     std::string lut_tag=ltag.str();
01850     std::cout<< " LUT tag "<<lut_tag<<endl; 
01851     if(m_write_lut==1) lut_conf_id_=db_->writeToConfDB_TPGLUT(dataset, dataset2,lutparamset, NLUTGROUPS, lut_tag) ;
01852 
01853   }
01854 
01855   // last we insert the FE_CONFIG_MAIN table 
01856  if (writeToDB_) {
01857    
01858    //int conf_id_=db_->writeToConfDB_TPGMain(ped_conf_id_,lin_conf_id_, lut_conf_id_, fgr_conf_id_, 
01859    //                           sli_conf_id_, wei_conf_id_, bxt_conf_id_, btt_conf_id_, tag_, version_) ;
01860    int conf_id_=db_->writeToConfDB_TPGMain(ped_conf_id_,lin_conf_id_, lut_conf_id_, fgr_conf_id_, 
01861                                            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
01862    
01863    std::cout << "\n Conf ID = " << conf_id_ << std::endl;
01864 
01865  }
01866 
01867 
01869   // loop on strips and associate them with  values //
01871 
01872   // Barrel
01873   stripListEB.sort() ;
01874   stripListEB.unique() ;
01875   cout<<"Number of EB strips="<<dec<<stripListEB.size()<<endl ;
01876   if (writeToFiles_) {
01877     (*out_file_) <<std::endl ;
01878     for (itList = stripListEB.begin(); itList != stripListEB.end(); itList++ ) {
01879       (*out_file_) <<"STRIP_EB "<<dec<<(*itList)<<endl ;
01880       (*out_file_) << hex << "0x" <<sliding_<<std::endl ;
01881       (*out_file_) <<"0" <<std::endl ;
01882       (*out_file_) <<"0x"<<stripMapEBsintheta[(*itList)] << " 0x" << SFGVB_lut_ <<std::endl ;
01883     }
01884   }
01885 
01886   // Endcap
01887   stripListEE.sort() ;
01888   stripListEE.unique() ;
01889   cout<<"Number of EE strips="<<dec<<stripListEE.size()<<endl ;
01890   if (writeToFiles_) {
01891     (*out_file_) <<std::endl ;
01892     for (itList = stripListEE.begin(); itList != stripListEE.end(); itList++ ) {
01893       (*out_file_) <<"STRIP_EE "<<dec<<(*itList)<<endl ;
01894       (*out_file_) << hex << "0x" <<sliding_<<std::endl ;
01895       (*out_file_) <<" 0" << std::endl ;
01896       (*out_file_)<<hex<<"0x"<<threshold<<" 0x"<<lut_strip<<std::endl ;  
01897     }
01898   }
01899 
01900 
01902   // loop on towers and associate them with default values //
01904 
01905   // Barrel
01906   towerListEB.sort() ;
01907   towerListEB.unique() ;
01908   cout<<"Number of EB towers="<<dec<<towerListEB.size()<<endl ;
01909   if (writeToFiles_) {
01910     (*out_file_) <<std::endl ;
01911     for (itList = towerListEB.begin(); itList != towerListEB.end(); itList++ ) {
01912       (*out_file_) <<"TOWER_EB "<<dec<<(*itList)<<endl ;
01913       (*out_file_) <<" 0\n 0\n" ;
01914       (*out_file_) <<" " <<  SFGVB_SpikeKillingThreshold_ << std::endl; //modif-alex
01915     }
01916   }
01917 
01918   // Endcap
01919   towerListEE.sort() ;
01920   towerListEE.unique() ;
01921   cout<<"Number of EE towers="<<dec<<towerListEE.size()<<endl ;
01922   if (writeToFiles_) {
01923     (*out_file_) <<std::endl ;
01924     for (itList = towerListEE.begin(); itList != towerListEE.end(); itList++ ) {
01925       (*out_file_) <<"TOWER_EE "<<dec<<(*itList)<<endl ;
01926       if (newLUT) (*out_file_) <<" 1\n" ;
01927       else  (*out_file_) <<" 0\n" ;
01928       (*out_file_)<<hex<<"0x"<<lut_tower<<std::endl ;
01929     }
01930   }
01931 
01932 
01933 
01934 
01936   // store control histos //
01938   ICEB->Write() ;
01939   tpgFactorEB->Write() ;
01940   ICEEPlus->Write() ;
01941   tpgFactorEEPlus->Write() ;  
01942   ICEEMinus->Write() ;
01943   tpgFactorEEMinus->Write() ;
01944   IC->Write() ;
01945   tpgFactor->Write() ;
01946   hshapeEB->Write() ;
01947   hshapeEE->Write() ;
01948   ntuple->Write() ;
01949   ntupleSpike->Write() ;
01950   saving.Close () ;
01951 }
01952 
01953 void EcalTPGParamBuilder::beginJob()
01954 {
01955   using namespace edm;
01956   using namespace std;
01957 
01958   std::cout<<"we are in beginJob"<<endl;
01959 
01960   create_header() ; 
01961 
01962   DetId eb(DetId::Ecal,EcalBarrel) ;
01963   DetId ee(DetId::Ecal,EcalEndcap) ;
01964 
01965   if (writeToFiles_) {
01966     (*out_file_)<<"PHYSICS_EB "<<dec<<eb.rawId()<<std::endl ;
01967     (*out_file_)<<Et_sat_EB_<<" "<<TTF_lowThreshold_EB_<<" "<<TTF_highThreshold_EB_<<std::endl ;
01968     (*out_file_)<<FG_lowThreshold_EB_<<" "<<FG_highThreshold_EB_<<" "
01969                   <<FG_lowRatio_EB_<<" "<<FG_highRatio_EB_<<std::endl ;
01970     //(*out_file_) << SFGVB_SpikeKillingThreshold_ << std::endl; //modif-alex02/02/2011
01971     (*out_file_) <<std::endl ;
01972 
01973     (*out_file_)<<"PHYSICS_EE "<<dec<<ee.rawId()<<std::endl ;
01974     (*out_file_)<<Et_sat_EE_<<" "<<TTF_lowThreshold_EE_<<" "<<TTF_highThreshold_EE_<<std::endl ;
01975     (*out_file_)<<FG_Threshold_EE_<<" "<<-1<<" "
01976                   <<-1<<" "<<-1<<std::endl ;
01977     (*out_file_) <<std::endl ;
01978   }
01979 
01980 }
01981 
01982 
01983 
01984 bool EcalTPGParamBuilder::computeLinearizerParam(double theta, double gainRatio, double calibCoeff, std::string subdet, int & mult , int & shift) 
01985 {
01986   /*
01987     Linearization coefficient are determined in order to satisfy:
01988     tpg(ADC_sat) = 1024
01989     where: 
01990     tpg() is a model of the linearized tpg response on 10b 
01991     ADC_sat is the number of ADC count corresponding the Et_sat, the maximum scale of the transverse energy
01992     
01993     Since we have:
01994     Et_sat = xtal_LSB * ADC_sat * gainRatio * calibCoeff * sin(theta)
01995     and a simple model of tpg() being given by:
01996     tpg(X) = [ (X*mult) >> (shift+2) ] >> (sliding+shiftDet) 
01997     we must satisfy:
01998     [ (Et_sat/(xtal_LSB * gainRatio * calibCoeff * sin(theta)) * mult) >> (shift+2) ] >> (sliding+shiftDet) = 1024 
01999     that is:
02000     mult = 1024/Et_sat * xtal_LSB * gainRatio * calibCoeff * sin(theta) * 2^-(sliding+shiftDet+2) * 2^-shift
02001     mult = factor * 2^-shift
02002   */
02003 
02004   // case barrel:
02005   int shiftDet = 2 ; //fixed, due to FE FENIX TCP format
02006   double ratio = xtal_LSB_EB_/Et_sat_EB_ ;
02007   // case endcap:
02008   if (subdet=="EE") {
02009     shiftDet = 2 ; //applied in TCC-EE and not in FE FENIX TCP... This parameters is setable in the TCC-EE
02010     //shiftDet = 0 ; //was like this before with FE bug
02011     ratio = xtal_LSB_EE_/Et_sat_EE_ ;
02012   }
02013 
02014 
02015 
02016   double factor = 1024 * ratio * gainRatio * calibCoeff * sin(theta) * (1 << (sliding_ + shiftDet + 2)) ;
02017   // Let's try first with shift = 0 (trivial solution)
02018   mult = (int)(factor+0.5) ; 
02019   for (shift = 0 ; shift<15 ; shift++) {
02020     if (mult>=128  && mult<256) return true ;
02021     factor *= 2 ; 
02022     mult = (int)(factor+0.5) ;
02023   }
02024   std::cout << "too bad we did not manage to calculate the factor for calib="<<calibCoeff<<endl;
02025   return false ;
02026 }
02027 
02028 void EcalTPGParamBuilder::create_header() 
02029 {
02030   if (!writeToFiles_) return ;
02031   (*out_file_) <<"COMMENT put your comments here"<<std::endl ; 
02032 
02033   (*out_file_) <<"COMMENT ================================="<<std::endl ;
02034   (*out_file_) <<"COMMENT           physics EB structure"<<std::endl ;
02035   (*out_file_) <<"COMMENT"<<std::endl ;
02036   (*out_file_) <<"COMMENT  EtSaturation (GeV), ttf_threshold_Low (GeV), ttf_threshold_High (GeV)"<<std::endl ;
02037   (*out_file_) <<"COMMENT  FG_lowThreshold (GeV), FG_highThreshold (GeV), FG_lowRatio, FG_highRatio"<<std::endl ;
02038   //(*out_file_) <<"COMMENT  SFGVB_SpikeKillingThreshold (GeV)"<<std::endl ; //modif-alex-02/02/2011
02039   (*out_file_) <<"COMMENT ================================="<<std::endl ;
02040   (*out_file_) <<"COMMENT"<<std::endl ;
02041 
02042   (*out_file_) <<"COMMENT ================================="<<std::endl ;
02043   (*out_file_) <<"COMMENT           physics EE structure"<<std::endl ;
02044   (*out_file_) <<"COMMENT"<<std::endl ;
02045   (*out_file_) <<"COMMENT  EtSaturation (GeV), ttf_threshold_Low (GeV), ttf_threshold_High (GeV)"<<std::endl ;
02046   (*out_file_) <<"COMMENT  FG_Threshold (GeV), dummy, dummy, dummy"<<std::endl ;
02047   (*out_file_) <<"COMMENT ================================="<<std::endl ;
02048   (*out_file_) <<"COMMENT"<<std::endl ;
02049 
02050   (*out_file_) <<"COMMENT ================================="<<std::endl ;
02051   (*out_file_) <<"COMMENT           crystal structure (same for EB and EE)"<<std::endl ;
02052   (*out_file_) <<"COMMENT"<<std::endl ;
02053   (*out_file_) <<"COMMENT  ped, mult, shift [gain12]"<<std::endl ;
02054   (*out_file_) <<"COMMENT  ped, mult, shift [gain6]"<<std::endl ;
02055   (*out_file_) <<"COMMENT  ped, mult, shift [gain1]"<<std::endl ;
02056   (*out_file_) <<"COMMENT ================================="<<std::endl ;
02057   (*out_file_) <<"COMMENT"<<std::endl ;
02058 
02059   (*out_file_) <<"COMMENT ================================="<<std::endl ;
02060   (*out_file_) <<"COMMENT           strip EB structure"<<std::endl ;
02061   (*out_file_) <<"COMMENT"<<std::endl ;
02062   (*out_file_) <<"COMMENT  sliding_window"<<std::endl ;
02063   (*out_file_) <<"COMMENT  weightGroupId"<<std::endl ;
02064   (*out_file_) <<"COMMENT  threshold_sfg lut_sfg"<<std::endl ;
02065   (*out_file_) <<"COMMENT ================================="<<std::endl ;
02066   (*out_file_) <<"COMMENT"<<std::endl ;
02067 
02068   (*out_file_) <<"COMMENT ================================="<<std::endl ;
02069   (*out_file_) <<"COMMENT           strip EE structure"<<std::endl ;
02070   (*out_file_) <<"COMMENT"<<std::endl ;
02071   (*out_file_) <<"COMMENT  sliding_window"<<std::endl ;
02072   (*out_file_) <<"COMMENT  weightGroupId"<<std::endl ;
02073   (*out_file_) <<"COMMENT  threshold_fg lut_fg"<<std::endl ;
02074   (*out_file_) <<"COMMENT ================================="<<std::endl ;
02075   (*out_file_) <<"COMMENT"<<std::endl ;
02076 
02077   (*out_file_) <<"COMMENT ================================="<<std::endl ;
02078   (*out_file_) <<"COMMENT           tower EB structure"<<std::endl ;
02079   (*out_file_) <<"COMMENT"<<std::endl ;
02080   (*out_file_) <<"COMMENT  LUTGroupId"<<std::endl ;
02081   (*out_file_) <<"COMMENT  FgGroupId"<<std::endl ;
02082   (*out_file_) <<"COMMENT  spike_killing_threshold"<<std::endl ;//modif alex
02083   (*out_file_) <<"COMMENT ================================="<<std::endl ;
02084   (*out_file_) <<"COMMENT"<<std::endl ;
02085 
02086   (*out_file_) <<"COMMENT ================================="<<std::endl ;
02087   (*out_file_) <<"COMMENT           tower EE structure"<<std::endl ;
02088   (*out_file_) <<"COMMENT"<<std::endl ;
02089   (*out_file_) <<"COMMENT  LUTGroupId"<<std::endl ;
02090   (*out_file_) <<"COMMENT  tower_lut_fg"<<std::endl ;
02091   (*out_file_) <<"COMMENT ================================="<<std::endl ;
02092   (*out_file_) <<"COMMENT"<<std::endl ;
02093   
02094   (*out_file_) <<"COMMENT ================================="<<std::endl ;
02095   (*out_file_) <<"COMMENT           Weight structure"<<std::endl ;
02096   (*out_file_) <<"COMMENT"<<std::endl ;
02097   (*out_file_) <<"COMMENT  weightGroupId"<<std::endl ;
02098   (*out_file_) <<"COMMENT  w0, w1, w2, w3, w4"<<std::endl ;
02099   (*out_file_) <<"COMMENT ================================="<<std::endl ;
02100   (*out_file_) <<"COMMENT"<<std::endl ;
02101 
02102   (*out_file_) <<"COMMENT ================================="<<std::endl ;
02103   (*out_file_) <<"COMMENT           lut structure"<<std::endl ;
02104   (*out_file_) <<"COMMENT"<<std::endl ;
02105   (*out_file_) <<"COMMENT  LUTGroupId"<<std::endl ;
02106   (*out_file_) <<"COMMENT  LUT[1-1024]"<<std::endl ;
02107   (*out_file_) <<"COMMENT ================================="<<std::endl ;
02108   (*out_file_) <<"COMMENT"<<std::endl ;
02109 
02110   (*out_file_) <<"COMMENT ================================="<<std::endl ;
02111   (*out_file_) <<"COMMENT           fg EB structure"<<std::endl ;
02112   (*out_file_) <<"COMMENT"<<std::endl ;
02113   (*out_file_) <<"COMMENT  FgGroupId"<<std::endl ;
02114   (*out_file_) <<"COMMENT  el, eh, tl, th, lut_fg"<<std::endl ;
02115   (*out_file_) <<"COMMENT ================================="<<std::endl ;
02116   (*out_file_) <<"COMMENT"<<std::endl ;
02117 
02118   (*out_file_) <<std::endl ;
02119 }
02120 
02121 
02122 int EcalTPGParamBuilder::uncodeWeight(double weight, int complement2)
02123 {
02124   int iweight ;
02125   unsigned int max = (unsigned int)(pow(2.,complement2)-1) ;
02126   if (weight>0) iweight=int((1<<6)*weight+0.5) ; // +0.5 for rounding pb
02127   else iweight= max - int(-weight*(1<<6)+0.5) +1 ;
02128   iweight = iweight & max ;
02129   return iweight ;
02130 }
02131 
02132 double EcalTPGParamBuilder::uncodeWeight(int iweight, int complement2)
02133 {
02134   double weight = double(iweight)/pow(2., 6.) ;
02135   // test if negative weight:
02136   if ( (iweight & (1<<(complement2-1))) != 0) weight = (double(iweight)-pow(2., complement2))/pow(2., 6.) ;
02137   return weight ;
02138 }
02139 
02140 #if (CMSSW_VERSION>=340)
02141 std::vector<unsigned int> EcalTPGParamBuilder::computeWeights(EcalShapeBase & shape, TH1F * histo)
02142 #else
02143 std::vector<unsigned int> EcalTPGParamBuilder::computeWeights(EcalShape & shape, TH1F * histo)
02144 #endif
02145 {
02146   std::cout<<"Computing Weights..."<<std::endl ;
02147 #if (CMSSW_VERSION>=340)
02148   double timeMax = shape.timeOfMax() - shape.timeOfThr() ; // timeMax w.r.t begining of pulse
02149 #else
02150   double timeMax = shape.computeTimeOfMaximum() - shape.computeT0() ; // timeMax w.r.t begining of pulse
02151 #endif
02152   double max = shape(timeMax) ;
02153 
02154   double sumf = 0. ;
02155   double sumf2 = 0. ;
02156   for (unsigned int sample = 0 ; sample<nSample_ ; sample++) {
02157     double time = timeMax - ((double)sampleMax_-(double)sample)*25. ;
02158     time -= weight_timeShift_ ;
02159     sumf += shape(time)/max ;
02160     sumf2 += shape(time)/max * shape(time)/max ;
02161     for (int subtime = 0 ; subtime<25 ; subtime++) histo->Fill(float(sample*25. + subtime)/25., shape(time+subtime)) ;
02162   }
02163   double lambda = 1./(sumf2-sumf*sumf/nSample_) ;
02164   double gamma = -lambda*sumf/nSample_ ;
02165   double * weight = new double[nSample_] ;
02166   for (unsigned int sample = 0 ; sample<nSample_ ; sample++) {
02167     double time = timeMax - ((double)sampleMax_-(double)sample)*25. ;
02168     time -= weight_timeShift_ ;
02169     weight[sample] = lambda*shape(time)/max + gamma ;
02170   }
02171 
02172 
02173   int * iweight = new int[nSample_] ;
02174   for (unsigned int sample = 0 ; sample<nSample_ ; sample++)   iweight[sample] = uncodeWeight(weight[sample], complement2_) ;
02175 
02176   // Let's check:  
02177   int isumw  = 0 ;  
02178   for (unsigned int sample = 0 ; sample<nSample_ ; sample++) isumw  += iweight[sample] ;
02179   unsigned int imax = (unsigned int)(pow(2.,int(complement2_))-1) ;
02180   isumw = (isumw & imax ) ;
02181 
02182   double ampl = 0. ;
02183   double sumw = 0. ;
02184   for (unsigned int sample = 0 ; sample<nSample_ ; sample++) {
02185      double time = timeMax - ((double)sampleMax_-(double)sample)*25. ;
02186      time -= weight_timeShift_ ;
02187      ampl += weight[sample]*shape(time) ;
02188      sumw += weight[sample] ;
02189      std::cout<<"weight="<<weight[sample]<<" shape="<<shape(time)<<std::endl ;
02190   }
02191   std::cout<<"Weights: sum="<<isumw<<" in float ="<<uncodeWeight(isumw, complement2_)<<" sum of floats ="<<sumw<<std::endl ;
02192   std::cout<<"Weights: sum (weight*shape) = "<<ampl<<std::endl ;
02193 
02194 
02195 
02196   // Let's correct for bias if any
02197   if (weight_unbias_recovery_) {
02198     int count = 0 ;
02199     while (isumw != 0 && count<10) {
02200       double min = 99. ;
02201       unsigned int index = 0 ;
02202       if ( (isumw & (1<<(complement2_-1))) != 0) {
02203         // add 1:
02204         std::cout<<"Correcting for bias: adding 1"<<std::endl ;
02205         for (unsigned int sample = 0 ; sample<nSample_ ; sample++) {
02206           int new_iweight = iweight[sample]+1 ; 
02207           double new_weight = uncodeWeight(new_iweight, complement2_) ;
02208           if (fabs(new_weight-weight[sample])<min) {
02209             min = fabs(new_weight-weight[sample]) ;
02210             index = sample ;
02211           }
02212         }
02213         iweight[index] ++ ; 
02214       } else {
02215         // Sub 1:
02216         std::cout<<"Correcting for bias: subtracting 1"<<std::endl ;
02217         for (unsigned int sample = 0 ; sample<nSample_ ; sample++) {
02218           int new_iweight = iweight[sample]-1 ;    
02219           double new_weight = uncodeWeight(new_iweight, complement2_) ;
02220           if (fabs(new_weight-weight[sample])<min) {
02221             min = fabs(new_weight-weight[sample]) ;
02222             index = sample ;
02223           }
02224         }
02225         iweight[index] -- ; 
02226       } 
02227       isumw  = 0 ;  
02228       for (unsigned int sample = 0 ; sample<nSample_ ; sample++) isumw  += iweight[sample] ; 
02229       imax = (unsigned int)(pow(2.,int(complement2_))-1) ;
02230       isumw = (isumw & imax ) ;
02231       std::cout<<"Correcting weight number: "<<index<<" sum weights = "<<isumw<<std::endl ;
02232       count ++ ;
02233     }
02234   }
02235   
02236   // let's check again
02237   isumw  = 0 ;  
02238   for (unsigned int sample = 0 ; sample<nSample_ ; sample++) isumw  += iweight[sample] ;
02239   imax = (unsigned int)(pow(2.,int(complement2_))-1) ;
02240   isumw = (isumw & imax ) ;
02241   ampl = 0. ;
02242   for (unsigned int sample = 0 ; sample<nSample_ ; sample++) {
02243      double time = timeMax - ((double)sampleMax_-(double)sample)*25. ;
02244      time -= weight_timeShift_ ;
02245      double new_weight = uncodeWeight(iweight[sample], complement2_) ;
02246      sumw += uncodeWeight(iweight[sample], complement2_) ;
02247      ampl += new_weight*shape(time) ;
02248      std::cout<<"weight unbiased after integer conversion="<<new_weight<<" shape="<<shape(time)<<std::endl ;
02249   }
02250   std::cout<<"Weights: sum="<<isumw<<" in float ="<<uncodeWeight(isumw, complement2_)<<" sum of floats ="<<sumw<<std::endl ;
02251   std::cout<<"Weights: sum (weight*shape) = "<<ampl<<std::endl ;
02252 
02253 
02254 
02255   std::vector<unsigned int> theWeights ;
02256   for (unsigned int sample = 0 ; sample<nSample_ ; sample++) theWeights.push_back(iweight[sample]) ;
02257   std::cout<<std::endl ;
02258 
02259   delete weight ;
02260   delete iweight ;
02261   return theWeights ;
02262 }
02263 
02264 void EcalTPGParamBuilder::computeLUT(int * lut, std::string det) 
02265 {
02266   double Et_sat = Et_sat_EB_ ;
02267   double LUT_threshold = LUT_threshold_EB_ ;
02268   double LUT_stochastic = LUT_stochastic_EB_ ;
02269   double LUT_noise = LUT_noise_EB_ ;
02270   double LUT_constant = LUT_constant_EB_ ;
02271   double TTF_lowThreshold = TTF_lowThreshold_EB_ ;
02272   double TTF_highThreshold = TTF_highThreshold_EB_ ;
02273   if (det == "EE") {
02274     Et_sat = Et_sat_EE_ ;
02275     LUT_threshold = LUT_threshold_EE_ ;
02276     LUT_stochastic = LUT_stochastic_EE_ ;
02277     LUT_noise = LUT_noise_EE_ ;
02278     LUT_constant = LUT_constant_EE_ ;
02279     TTF_lowThreshold = TTF_lowThreshold_EE_ ;
02280     TTF_highThreshold = TTF_highThreshold_EE_ ;
02281   }
02282 
02283   // initialisation with identity
02284   for (int i=0 ; i<1024 ; i++) {
02285     lut[i] = i ;
02286     if (lut[i]>0xff) lut[i] = 0xff ;
02287   }
02288 
02289   // case linear LUT
02290   if (LUT_option_ == "Linear") {
02291     int mylut = 0 ;
02292     for (int i=0 ; i<1024 ; i++) {
02293       lut[i] = mylut ;
02294       if ((i+1)%4 == 0 ) mylut++ ; 
02295       //if ((i+1)%8 == 0 ) mylut++ ;//modif-alex 16/12/2010 LSB==500MeV ONLY USED FOR BEAMV4 key 
02296     }
02297   }
02298 
02299   // case LUT following Ecal resolution
02300   if (LUT_option_ == "EcalResolution") {
02301     TF1 * func = new TF1("func",oneOverEtResolEt, 0., Et_sat,3) ;
02302     func->SetParameters(LUT_stochastic, LUT_noise, LUT_constant) ;
02303     double norm = func->Integral(0., Et_sat) ;
02304     for (int i=0 ; i<1024 ; i++) {   
02305       double Et = i*Et_sat/1024. ;
02306       lut[i] =  int(0xff*func->Integral(0., Et)/norm + 0.5) ;
02307     }
02308   }
02309 
02310   // Now, add TTF thresholds to LUT and apply LUT threshold if needed
02311   for (int j=0 ; j<1024 ; j++) {
02312     double Et_GeV = Et_sat/1024*(j+0.5) ;
02313     if (Et_GeV <= LUT_threshold) lut[j] = 0 ; // LUT threshold
02314     int ttf = 0x0 ;    
02315     if (Et_GeV >= TTF_highThreshold) ttf = 3 ;
02316     if (Et_GeV >= TTF_lowThreshold && Et_GeV < TTF_highThreshold) ttf = 1 ;
02317     ttf = ttf << 8 ;
02318     lut[j] += ttf ;
02319   }
02320 
02321 }
02322 
02323 void EcalTPGParamBuilder::getCoeff(coeffStruc & coeff, const EcalIntercalibConstantMap & calibMap, unsigned int rawId)
02324 {
02325   // get current intercalibration coeff
02326   coeff.calibCoeff_ = 1. ;
02327   if (!useInterCalibration_) return ;
02328   EcalIntercalibConstantMap::const_iterator icalit = calibMap.find(rawId);
02329   if( icalit != calibMap.end() ) coeff.calibCoeff_ = (*icalit) ;
02330   else std::cout<<"getCoeff: "<<rawId<<" not found in EcalIntercalibConstantMap"<<std::endl ;
02331 }
02332 
02333 void EcalTPGParamBuilder::getCoeff(coeffStruc & coeff, const EcalGainRatioMap & gainMap, unsigned int rawId)
02334 {
02335   // get current gain ratio
02336   coeff.gainRatio_[0]  = 1. ;
02337   coeff.gainRatio_[1]  = 2. ;
02338   coeff.gainRatio_[2]  = 12. ;
02339   EcalGainRatioMap::const_iterator gainIter = gainMap.find(rawId);
02340   if (gainIter != gainMap.end()) {
02341     const EcalMGPAGainRatio & aGain = (*gainIter) ;
02342     coeff.gainRatio_[1] = aGain.gain12Over6() ;
02343     coeff.gainRatio_[2] = aGain.gain6Over1() * aGain.gain12Over6() ;
02344   }
02345   else std::cout<<"getCoeff: "<<rawId<<" not found in EcalGainRatioMap"<<std::endl ;
02346 }
02347 
02348 void EcalTPGParamBuilder::getCoeff(coeffStruc & coeff, const EcalPedestalsMap & pedMap, unsigned int rawId)
02349 {
02350   coeff.pedestals_[0] = 0 ;
02351   coeff.pedestals_[1] = 0 ;
02352   coeff.pedestals_[2] = 0 ;
02353 
02354   if (forcedPedestalValue_ >= 0) {
02355     coeff.pedestals_[0] = forcedPedestalValue_ ;
02356     coeff.pedestals_[1] = forcedPedestalValue_ ;
02357     coeff.pedestals_[2] = forcedPedestalValue_ ;  
02358     return ;
02359   }
02360 
02361   // get current pedestal
02362   EcalPedestalsMapIterator pedIter = pedMap.find(rawId);
02363   if (pedIter != pedMap.end()) {
02364     EcalPedestals::Item aped = (*pedIter);
02365     coeff.pedestals_[0] = int(aped.mean_x12 + 0.5) ; 
02366     coeff.pedestals_[1] = int(aped.mean_x6 + 0.5) ;
02367     coeff.pedestals_[2] = int(aped.mean_x1 + 0.5) ;
02368   }
02369   else std::cout<<"getCoeff: "<<rawId<<" not found in EcalPedestalsMap"<<std::endl ;
02370 }
02371 
02372 void EcalTPGParamBuilder::getCoeff(coeffStruc & coeff, const map<EcalLogicID, MonPedestalsDat> & pedMap, const EcalLogicID & logicId)
02373 {
02374   // get current pedestal
02375   coeff.pedestals_[0] = 0 ;
02376   coeff.pedestals_[1] = 0 ;
02377   coeff.pedestals_[2] = 0 ;
02378 
02379   map<EcalLogicID, MonPedestalsDat>::const_iterator it =  pedMap.find(logicId);
02380   if (it != pedMap.end()) {
02381     MonPedestalsDat ped = it->second ;
02382     coeff.pedestals_[0] = int(ped.getPedMeanG12() + 0.5) ; 
02383     coeff.pedestals_[1] = int(ped.getPedMeanG6() + 0.5) ; 
02384     coeff.pedestals_[2] = int(ped.getPedMeanG1() + 0.5) ; 
02385   } 
02386   else std::cout<<"getCoeff: "<<logicId.getID1()<<", "<<logicId.getID2()<<", "<<logicId.getID3()
02387                 <<" not found in map<EcalLogicID, MonPedestalsDat>"<<std::endl ;
02388 }
02389 
02390 void EcalTPGParamBuilder::computeFineGrainEBParameters(unsigned int & lowRatio, unsigned int & highRatio,
02391                                                        unsigned int & lowThreshold, unsigned int & highThreshold, unsigned int & lut)
02392 {
02393   lowRatio = int(0x80*FG_lowRatio_EB_ + 0.5) ;
02394   if (lowRatio>0x7f) lowRatio = 0x7f ;
02395   highRatio = int(0x80*FG_highRatio_EB_ + 0.5) ;
02396   if (highRatio>0x7f) highRatio = 0x7f ;
02397   
02398   // lsb at the stage of the FG calculation is:
02399   double lsb_FG = Et_sat_EB_/1024./4 ;
02400   lowThreshold = int(FG_lowThreshold_EB_/lsb_FG+0.5) ;
02401   if (lowThreshold>0xff) lowThreshold = 0xff ;
02402   highThreshold = int(FG_highThreshold_EB_/lsb_FG+0.5) ;
02403   if (highThreshold>0xff) highThreshold = 0xff ;
02404 
02405   // FG lut: FGVB response is LUT(adress) where adress is: 
02406   // bit3: maxof2/ET >= lowRatio, bit2: maxof2/ET >= highRatio, bit1: ET >= lowThreshold, bit0: ET >= highThreshold
02407   // FGVB =1 if jet-like (veto active), =0 if E.M.-like
02408   // the condition for jet-like is: ET>Threshold and  maxof2/ET < Ratio (only TT with enough energy are vetoed)
02409 
02410   // With the following lut, what matters is only max(TLow, Thigh) and max(Elow, Ehigh)
02411   // So, jet-like if maxof2/ettot<max(TLow, Thigh) && ettot >= max(Elow, Ehigh)
02412   if (FG_lut_EB_ == 0) lut = 0x0888 ; 
02413   else lut = FG_lut_EB_ ; // let's use the users value (hope he/she knows what he/she does!)
02414 }
02415 
02416 void EcalTPGParamBuilder::computeFineGrainEEParameters(unsigned int & threshold, unsigned int & lut_strip, unsigned int & lut_tower) 
02417 {
02418   // lsb for EE:
02419   double lsb_FG = Et_sat_EE_/1024. ; // FIXME is it true????
02420   threshold = int(FG_Threshold_EE_/lsb_FG+0.5) ;
02421   lut_strip = FG_lut_strip_EE_  ;
02422   lut_tower = FG_lut_tower_EE_  ;
02423 }
02424 
02425 bool EcalTPGParamBuilder::realignBaseline(linStruc & lin, float forceBase12)
02426 {
02427   bool ok(true) ;
02428   float base[3] = {forceBase12, float(lin.pedestal_[1]), float(lin.pedestal_[2])} ;
02429   for (int i=1 ; i<3 ; i++) {
02430     if (lin.mult_[i]>0) {
02431       base[i] = float(lin.pedestal_[i]) - 
02432         float(lin.mult_[0])/float(lin.mult_[i])*pow(2., -(lin.shift_[0]-lin.shift_[i]))*(lin.pedestal_[0]-base[0]) ;
02433     } else base[i] = 0 ;
02434   }
02435 
02436   for (int i=0 ; i<3 ; i++) {
02437     lin.pedestal_[i] = base[i] ;
02438     //cout<<lin.pedestal_[i]<<" "<<base[i]<<endl ;
02439     if (base[i]<0 || lin.pedestal_[i]>1000) {
02440       cout<<"WARNING: base= "<<base[i]<<", "<<lin.pedestal_[i]<<" for gainId[0-2]="<<i<<" ==> forcing at 0"<<endl ;
02441       lin.pedestal_[i] = 0 ; 
02442       ok = false ;
02443     }
02444   }
02445   return ok ;
02446 
02447 }
02448 
02449 int EcalTPGParamBuilder::getGCTRegionPhi(int ttphi){
02450   int gctphi=0;
02451   gctphi = (ttphi+1)/4;
02452   if(ttphi<=2) gctphi=0;
02453   if(ttphi>=71) gctphi=0;
02454   return gctphi;
02455 }
02456 
02457 int EcalTPGParamBuilder::getGCTRegionEta(int tteta){
02458   int gcteta = 0;
02459   if(tteta>0) gcteta = (tteta-1)/4 + 11;
02460   else if(tteta<0) gcteta = (tteta+1)/4 + 10;
02461   return gcteta;
02462 }
02463 
02464 std::string EcalTPGParamBuilder::getDet(int tcc){
02465   std::stringstream sdet ;
02466 
02467   if (tcc>36 && tcc<55) sdet<<"EB-"<<tcc-36 ;
02468   else if (tcc>=55 && tcc<73) sdet<<"EB+"<<tcc-54 ;
02469   else if (tcc<=36) sdet<<"EE-" ;
02470   else sdet<<"EE+" ;
02471   
02472   if (tcc<=36 || tcc>=73) {
02473     if (tcc>=73) tcc-=72 ;
02474     if (tcc==1 || tcc==18 || tcc==19 || tcc==36)  sdet<<7 ;
02475     else if (tcc==2 || tcc==3 || tcc==20 || tcc==21)  sdet<<8 ;
02476     else if (tcc==4 || tcc==5 || tcc==22 || tcc==23)  sdet<<9 ;
02477     else if (tcc==6 || tcc==7 || tcc==24 || tcc==25)  sdet<<1 ;
02478     else if (tcc==8 || tcc==9 || tcc==26 || tcc==27)  sdet<<2 ;
02479     else if (tcc==10 || tcc==11 || tcc==28 || tcc==29)  sdet<<3 ;
02480     else if (tcc==12 || tcc==13 || tcc==30 || tcc==31)  sdet<<4 ;
02481     else if (tcc==14 || tcc==15 || tcc==32 || tcc==33)  sdet<<5 ;
02482     else if (tcc==16 || tcc==17 || tcc==34 || tcc==35)  sdet<<6 ;
02483   }
02484 
02485   return sdet.str() ;
02486 }
02487 
02488 std::pair < std::string, int >  EcalTPGParamBuilder::getCrate(int tcc){
02489   std::stringstream crate ;
02490   std::string pos ;
02491   int slot = 0 ;
02492 
02493   crate<<"S2D" ;  
02494   if (tcc>=40 && tcc<=42) {crate<<"02d" ; slot = 5 + (tcc-40)*6 ;}
02495   if (tcc>=43 && tcc<=45) {crate<<"03d" ; slot = 5 + (tcc-43)*6 ;}
02496   if (tcc>=37 && tcc<=39) {crate<<"04d" ; slot = 5 + (tcc-37)*6 ;}
02497   if (tcc>=52 && tcc<=54) {crate<<"06d" ; slot = 5 + (tcc-52)*6 ;}
02498   if (tcc>=46 && tcc<=48) {crate<<"07d" ; slot = 5 + (tcc-46)*6 ;}
02499   if (tcc>=49 && tcc<=51) {crate<<"08d" ; slot = 5 + (tcc-49)*6 ;}
02500   if (tcc>=58 && tcc<=60) {crate<<"02h" ; slot = 5 + (tcc-58)*6 ;}
02501   if (tcc>=61 && tcc<=63) {crate<<"03h" ; slot = 5 + (tcc-61)*6 ;}
02502   if (tcc>=55 && tcc<=57) {crate<<"04h" ; slot = 5 + (tcc-55)*6 ;}
02503   if (tcc>=70 && tcc<=72) {crate<<"06h" ; slot = 5 + (tcc-70)*6 ;}
02504   if (tcc>=64 && tcc<=66) {crate<<"07h" ; slot = 5 + (tcc-64)*6 ;}
02505   if (tcc>=67 && tcc<=69) {crate<<"08h" ; slot = 5 + (tcc-67)*6 ;}
02506 
02507   if (tcc>=76 && tcc<=81) {
02508     crate<<"02l" ; 
02509     if (tcc%2==0) slot = 2 + (tcc-76)*3 ;
02510     else slot = 4 + (tcc-77)*3 ;
02511   }
02512   if (tcc>=94 && tcc<=99) {
02513     crate<<"02l" ; 
02514     if (tcc%2==0) slot = 3 + (tcc-94)*3 ;
02515     else slot = 5 + (tcc-95)*3 ;
02516   }
02517 
02518   if (tcc>=22 && tcc<=27) {
02519     crate<<"03l" ; 
02520     if (tcc%2==0) slot = 2 + (tcc-22)*3 ;
02521     else slot = 4 + (tcc-23)*3 ;
02522   }
02523   if (tcc>=4 && tcc<=9) {
02524     crate<<"03l" ; 
02525     if (tcc%2==0) slot = 3 + (tcc-4)*3 ;
02526     else slot = 5 + (tcc-5)*3 ;
02527   }
02528 
02529   if (tcc>=82 && tcc<=87) {
02530     crate<<"07l" ; 
02531     if (tcc%2==0) slot = 2 + (tcc-82)*3 ;
02532     else slot = 4 + (tcc-83)*3 ;
02533   }
02534   if (tcc>=100 && tcc<=105) {
02535     crate<<"07l" ; 
02536     if (tcc%2==0) slot = 3 + (tcc-100)*3 ;
02537     else slot = 5 + (tcc-101)*3 ;
02538   }
02539 
02540   if (tcc>=28 && tcc<=33) {
02541     crate<<"08l" ; 
02542     if (tcc%2==0) slot = 2 + (tcc-28)*3 ;
02543     else slot = 4 + (tcc-29)*3 ;
02544   }
02545   if (tcc>=10 && tcc<=15) {
02546     crate<<"08l" ; 
02547     if (tcc%2==0) slot = 3 + (tcc-10)*3 ;
02548     else slot = 5 + (tcc-11)*3 ;
02549   }
02550 
02551   if (tcc==34) {crate<<"04l" ; slot = 2 ;}
02552   if (tcc==16) {crate<<"04l" ; slot = 3 ;}
02553   if (tcc==35) {crate<<"04l" ; slot = 4 ;}
02554   if (tcc==17) {crate<<"04l" ; slot = 5 ;}
02555   if (tcc==36) {crate<<"04l" ; slot = 8 ;}
02556   if (tcc==18) {crate<<"04l" ; slot = 9 ;}
02557   if (tcc==19) {crate<<"04l" ; slot = 10 ;}
02558   if (tcc==1)  {crate<<"04l" ; slot = 11 ;}
02559   if (tcc==20) {crate<<"04l" ; slot = 14 ;}
02560   if (tcc==2)  {crate<<"04l" ; slot = 15 ;}
02561   if (tcc==21) {crate<<"04l" ; slot = 16 ;}
02562   if (tcc==3)  {crate<<"04l" ; slot = 17 ;}
02563 
02564   if (tcc==88)  {crate<<"06l" ; slot = 2 ;}
02565   if (tcc==106) {crate<<"06l" ; slot = 3 ;}
02566   if (tcc==89)  {crate<<"06l" ; slot = 4 ;}
02567   if (tcc==107) {crate<<"06l" ; slot = 5 ;}
02568   if (tcc==90)  {crate<<"06l" ; slot = 8 ;}
02569   if (tcc==108) {crate<<"06l" ; slot = 9 ;}
02570   if (tcc==73)  {crate<<"06l" ; slot = 10 ;}
02571   if (tcc==91)  {crate<<"06l" ; slot = 11 ;}
02572   if (tcc==74)  {crate<<"06l" ; slot = 14 ;}
02573   if (tcc==92)  {crate<<"06l" ; slot = 15 ;}
02574   if (tcc==75)  {crate<<"06l" ; slot = 16 ;}
02575   if (tcc==93)  {crate<<"06l" ; slot = 17 ;}
02576 
02577   return std::pair< std::string, int > (crate.str(),slot) ;
02578 }