Tools/bin/coeffCompare.cpp

0001 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0002 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0003 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0004 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0005 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0006 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0007 #include "CondFormats/DataRecord/interface/EcalIntercalibConstantsRcd.h"
0008 #include "Calibration/Tools/interface/calibXMLwriter.h"
0009 #include "CalibCalorimetry/CaloMiscalibTools/interface/CaloMiscalibTools.h"
0010 #include "CalibCalorimetry/CaloMiscalibTools/interface/CaloMiscalibMapEcal.h"
0011 #include "CalibCalorimetry/CaloMiscalibTools/interface/MiscalibReaderFromXMLEcalBarrel.h"
0012 #include "CalibCalorimetry/CaloMiscalibTools/interface/MiscalibReaderFromXMLEcalEndcap.h"
0013 #include "CondFormats/EcalObjects/interface/EcalIntercalibConstants.h"
0014 #include "FWCore/Framework/interface/EventSetup.h"
0015 #include "FWCore/Framework/interface/ESHandle.h"
0016
0017 //#include "Calibration/EcalAlCaRecoProducers/interface/trivialParser.h"

0018 //#include "Calibration/Tools/bin/trivialParser.h"

0019
0020 #include "TH2.h"
0021 #include "TH1.h"
0022 #include "TFile.h"
0023
0024 #define PI_GRECO 3.14159265
0025
0026 inline int etaShifter(const int etaOld) {
0027   if (etaOld < 0)
0028     return etaOld + 85;
0029   else if (etaOld > 0)
0030     return etaOld + 84;
0031   assert(false && "eta cannot be 0");
0032 }
0033
0034 // ------------------------------------------------------------------------

0035
0036 // MF questo file prende due set di coefficienti e li confronta

0037 //

0038
0039 //-------------------------------------------------------------------------

0040
0041 int main(int argc, char* argv[]) {
0042   //  std::cout << "parsing cfg file: " << argv[1] << std::endl ;

0043   //  trivialParser configParams (static_cast<std::string> (argv[1])) ;

0044   //

0045   //  int EBetaStart = static_cast<int> (configParams.getVal ("EBetaStart")) ;

0046   //  int EBetaEnd = static_cast<int> (configParams.getVal ("EBetaEnd")) ;

0047   //  int EBphiStart = static_cast<int> (configParams.getVal ("EBphiStart")) ;

0048   //  int EBphiEnd = static_cast<int> (configParams.getVal ("EBphiEnd")) ;

0049   //  int EEradStart = static_cast<int> (configParams.getVal ("EEradStart")) ;

0050   //  int EEradEnd = static_cast<int> (configParams.getVal ("EEradEnd")) ;

0051   //  int EEphiStart = static_cast<int> (configParams.getVal ("EEphiStart")) ;

0052   //  int EEphiEnd = static_cast<int> (configParams.getVal ("EEphiEnd")) ;

0053
0054   //  int EBetaStart = -85 ;

0055   //  int EBetaEnd = 86 ;

0056   //  int EBphiStart = 1 ;

0057   //  int EBphiEnd = 361 ;

0058   //  int EEradStart = 0 ;

0059   //  int EEradEnd = 50 ;

0060   //  int EEphiStart = 1 ;

0061   //  int EEphiEnd = 361 ;

0062
0063   int EBetaStart = 0;
0064   int EBetaEnd = 86;
0065   int EBphiStart = 20;
0066   int EBphiEnd = 60;
0067   int EEradStart = 15;
0068   int EEradEnd = 50;
0069   int EEphiStart = 20;
0070   int EEphiEnd = 60;
0071
0072   //  std::string barrelfile ="miscalib_barrel_0.05.xml" ;

0073   //  std::string endcapfile ="miscalib_endcap_0.05.xml" ;

0074   //  std::string calibBarrelfile ="recalib_calibBarrel_0.05.xml" ;

0075   //  std::string calibEndcapfile ="recalib_calibEndcap_0.05.xml" ;

0076
0077   std::string barrelfile =
0078       "/afs/cern.ch/user/g/govoni/scratch1/CMSSW/CALIB/CMSSW_1_6_0/src/CalibCalorimetry/CaloMiscalibTools/data/"
0079       "ecal_barrel_startup.xml";
0080   std::string endcapfile =
0081       "/afs/cern.ch/user/g/govoni/scratch1/CMSSW/CALIB/CMSSW_1_6_0/src/CalibCalorimetry/CaloMiscalibTools/data/"
0082       "ecal_endcap_startup.xml";
0083   std::string calibBarrelfile =
0084       "/afs/cern.ch/user/g/govoni/scratch1/CMSSW/CALIB/CMSSW_1_6_0/src/CalibCalorimetry/CaloMiscalibTools/data/"
0085       "inv_ecal_barrel_startup.xml";
0086   std::string calibEndcapfile =
0087       "/afs/cern.ch/user/g/govoni/scratch1/CMSSW/CALIB/CMSSW_1_6_0/src/CalibCalorimetry/CaloMiscalibTools/data/"
0088       "inv_ecal_endcap_startup.xml";
0089
0090   //PG get the miscalibration files for EB and EE

0091   //PG ------------------------------------------

0092
0093   CaloMiscalibMapEcal EBscalibMap;
0094   EBscalibMap.prefillMap();
0095   MiscalibReaderFromXMLEcalBarrel barrelreader(EBscalibMap);
0096   if (!barrelfile.empty())
0097     barrelreader.parseXMLMiscalibFile(barrelfile);
0098   EcalIntercalibConstants* EBconstants = new EcalIntercalibConstants(EBscalibMap.get());
0099   const EcalIntercalibConstantMap& iEBscalibMap = EBconstants->getMap();  //MF prende i vecchi coeff

0100
0101   CaloMiscalibMapEcal EEscalibMap;
0102   EEscalibMap.prefillMap();
0103   MiscalibReaderFromXMLEcalEndcap endcapreader(EEscalibMap);
0104   if (!endcapfile.empty())
0105     endcapreader.parseXMLMiscalibFile(endcapfile);
0106   EcalIntercalibConstants* EEconstants = new EcalIntercalibConstants(EEscalibMap.get());
0107   const EcalIntercalibConstantMap& iEEscalibMap = EEconstants->getMap();  //MF prende i vecchi coeff

0108
0109   //PG get the recalibration files for EB and EE

0110   //PG -----------------------------------------

0111
0112   CaloMiscalibMapEcal EBcalibMap;
0113   EBcalibMap.prefillMap();
0114   MiscalibReaderFromXMLEcalBarrel calibBarrelreader(EBcalibMap);
0115   if (!calibBarrelfile.empty())
0116     calibBarrelreader.parseXMLMiscalibFile(calibBarrelfile);
0117   EcalIntercalibConstants* EBCconstants = new EcalIntercalibConstants(EBcalibMap.get());
0118   const EcalIntercalibConstantMap& iEBcalibMap = EBCconstants->getMap();  //MF prende i vecchi coeff

0119
0120   CaloMiscalibMapEcal EEcalibMap;
0121   EEcalibMap.prefillMap();
0122   MiscalibReaderFromXMLEcalEndcap calibEndcapreader(EEcalibMap);
0123   if (!calibEndcapfile.empty())
0124     calibEndcapreader.parseXMLMiscalibFile(calibEndcapfile);
0125   EcalIntercalibConstants* EECconstants = new EcalIntercalibConstants(EEcalibMap.get());
0126   const EcalIntercalibConstantMap& iEEcalibMap = EECconstants->getMap();  //MF prende i vecchi coeff

0127
0128   //PG fill the histograms

0129   //PG -------------------

0130
0131   // ECAL barrel

0132   TH1F EBCompareCoeffDistr("EBCompareCoeff", "EBCompareCoeff", 300, 0, 2);
0133   TH2F EBCompareCoeffMap("EBCompareCoeffMap", "EBCompareCoeffMap", 171, -85, 86, 360, 1, 361);
0134
0135   // ECAL barrel

0136   for (int ieta = EBetaStart; ieta < EBetaEnd; ++ieta)
0137     for (int iphi = EBphiStart; iphi <= EBphiEnd; ++iphi) {
0138       if (!EBDetId::validDetId(ieta, iphi))
0139         continue;
0140       EBDetId det = EBDetId(ieta, iphi, EBDetId::ETAPHIMODE);
0141       double factor = *(iEBcalibMap.find(det.rawId())) * *(iEBscalibMap.find(det.rawId()));
0142       EBCompareCoeffDistr.Fill(factor);
0143       EBCompareCoeffMap.Fill(ieta, iphi, factor);
0144     }  // ECAL barrel

0145
0146   TH1F EEPCompareCoeffDistr("EEPCompareCoeffDistr", "EEPCompareCoeffDistr", 200, 0, 2);
0147   TH2F EEPCompareCoeffMap("EEPCompareCoeffMap", "EEPCompareCoeffMap", 101, 0, 101, 101, 0, 101);
0148
0149   // ECAL endcap +

0150   for (int ix = 1; ix <= 100; ++ix)
0151     for (int iy = 1; iy <= 100; ++iy) {
0152       int rad = static_cast<int>(sqrt((ix - 50) * (ix - 50) + (iy - 50) * (iy - 50)));
0153       if (rad < EEradStart || rad > EEradEnd)
0154         continue;
0155       double phiTemp = atan2(iy - 50, ix - 50);
0156       if (phiTemp < 0)
0157         phiTemp += 2 * PI_GRECO;
0158       int phi = static_cast<int>(phiTemp * 180 / PI_GRECO);
0159       if (phi < EEphiStart || phi > EEphiEnd)
0160         continue;
0161       if (!EEDetId::validDetId(ix, iy, 1))
0162         continue;
0163       EEDetId det = EEDetId(ix, iy, 1, EEDetId::XYMODE);
0164       double factor = *(iEEcalibMap.find(det.rawId())) * *(iEEscalibMap.find(det.rawId()));
0165       EEPCompareCoeffDistr.Fill(factor);
0166       EEPCompareCoeffMap.Fill(ix, iy, factor);
0167     }  // ECAL endcap +

0168
0169   // ECAL endcap-

0170   TH1F EEMCompareCoeffDistr("EEMCompareCoeffDistr", "EEMCompareCoeffDistr", 200, 0, 2);
0171   TH2F EEMCompareCoeffMap("EEMCompareCoeffMap", "EEMCompareCoeffMap", 100, 0, 100, 100, 0, 100);
0172
0173   // ECAL endcap -

0174   for (int ix = 1; ix <= 100; ++ix)
0175     for (int iy = 1; iy <= 100; ++iy) {
0176       int rad = static_cast<int>(sqrt((ix - 50) * (ix - 50) + (iy - 50) * (iy - 50)));
0177       if (rad < EEradStart || rad > EEradEnd)
0178         continue;
0179       if (!EEDetId::validDetId(ix, iy, -1))
0180         continue;
0181       EEDetId det = EEDetId(ix, iy, -1, EEDetId::XYMODE);
0182       double factor = *(iEEcalibMap.find(det.rawId())) * *(iEEscalibMap.find(det.rawId()));
0183       EEMCompareCoeffDistr.Fill(factor);
0184       EEMCompareCoeffMap.Fill(ix, iy, factor);
0185     }  // ECAL endcap -

0186
0187   std::string filename = "coeffcompare.root";
0188   TFile out(filename.c_str(), "recreate");
0189   EEMCompareCoeffMap.Write();
0190   EEMCompareCoeffDistr.Write();
0191   EEPCompareCoeffMap.Write();
0192   EEPCompareCoeffDistr.Write();
0193   EBCompareCoeffMap.Write();
0194   EBCompareCoeffDistr.Write();
0195   out.Close();
0196 }