Tools/bin/coeffCompareInvEE.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/EcalAlCaRecoProducers/bin/trivialParser.h"

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

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

0038 //

0039
0040 //-------------------------------------------------------------------------

0041
0042 int main(int argc, char* argv[]) {
0043   int EEradStart = 15;
0044   int EEradEnd = 50;
0045   int EEphiStart = 20;
0046   int EEphiEnd = 45;
0047
0048   std::string endcapfile = argv[1];
0049   std::string calibEndcapfile = argv[2];
0050
0051   //  std::string endcapfile = "/afs/cern.ch/user/g/govoni/scratch1/CMSSW/CALIB/CMSSW_1_6_0/src/CalibCalorimetry/CaloMiscalibTools/data/ecal_endcap_startup.xml" ;

0052   //  std::string calibEndcapfile = "/afs/cern.ch/user/g/govoni/scratch1/CMSSW/CALIB/CMSSW_1_6_0/src/CalibCalorimetry/CaloMiscalibTools/data/inv_ecal_endcap_startup.xml" ;

0053
0054   //PG get the miscalibration files for EB and EE

0055   //PG ------------------------------------------

0056
0057   CaloMiscalibMapEcal EEscalibMap;
0058   EEscalibMap.prefillMap();
0059   MiscalibReaderFromXMLEcalEndcap endcapreader(EEscalibMap);
0060   if (!endcapfile.empty())
0061     endcapreader.parseXMLMiscalibFile(endcapfile);
0062   EcalIntercalibConstants* EEconstants = new EcalIntercalibConstants(EEscalibMap.get());
0063   const EcalIntercalibConstantMap& iEEscalibMap = EEconstants->getMap();  //MF prende i vecchi coeff

0064
0065   //PG get the recalibration files for EB and EE

0066   //PG -----------------------------------------

0067
0068   CaloMiscalibMapEcal EEcalibMap;
0069   EEcalibMap.prefillMap();
0070   MiscalibReaderFromXMLEcalEndcap calibEndcapreader(EEcalibMap);
0071   if (!calibEndcapfile.empty())
0072     calibEndcapreader.parseXMLMiscalibFile(calibEndcapfile);
0073   EcalIntercalibConstants* EECconstants = new EcalIntercalibConstants(EEcalibMap.get());
0074   const EcalIntercalibConstantMap& iEEcalibMap = EECconstants->getMap();  //MF prende i vecchi coeff

0075
0076   //PG fill the histograms

0077   //PG -------------------

0078
0079   TH1F EEPCompareCoeffDistr("EEPCompareCoeffDistr", "EEPCompareCoeffDistr", 5000, 0, 2);
0080   TH2F EEPCompareCoeffMap("EEPCompareCoeffMap", "EEPCompareCoeffMap", 101, 0, 101, 101, 0, 101);
0081   TH2F EEPCompareCoeffEtaTrend("EEPCompareCoeffEtaTrend", "EEPCompareCoeffEtaTrend", 51, 0, 50, 500, 0, 2);
0082   TProfile EEPCompareCoeffEtaProfile("EEPCompareCoeffEtaProfile", "EEPCompareCoeffEtaProfile", 51, 0, 50, 0, 2);
0083
0084   // ECAL endcap +

0085   for (int ix = 1; ix <= 100; ++ix)
0086     for (int iy = 1; iy <= 100; ++iy) {
0087       int rad = static_cast<int>(sqrt((ix - 50) * (ix - 50) + (iy - 50) * (iy - 50)));
0088       if (rad < EEradStart || rad > EEradEnd)
0089         continue;
0090       double phiTemp = atan2(iy - 50, ix - 50);
0091       if (phiTemp < 0)
0092         phiTemp += 2 * PI_GRECO;
0093       int phi = static_cast<int>(phiTemp * 180 / PI_GRECO);
0094       if (phi < EEphiStart || phi > EEphiEnd)
0095         continue;
0096       if (!EEDetId::validDetId(ix, iy, 1))
0097         continue;
0098       EEDetId det = EEDetId(ix, iy, 1, EEDetId::XYMODE);
0099       if (*(iEEscalibMap.find(det.rawId())) == 0)
0100         continue;
0101       double factor = *(iEEcalibMap.find(det.rawId())) / *(iEEscalibMap.find(det.rawId()));
0102       EEPCompareCoeffDistr.Fill(factor);
0103       EEPCompareCoeffMap.Fill(ix, iy, factor);
0104       EEPCompareCoeffEtaTrend.Fill(rad, factor);
0105       EEPCompareCoeffEtaProfile.Fill(rad, factor);
0106     }  // ECAL endcap +

0107
0108   // ECAL endcap-

0109   TH1F EEMCompareCoeffDistr("EEMCompareCoeffDistr", "EEMCompareCoeffDistr", 200, 0, 2);
0110   TH2F EEMCompareCoeffMap("EEMCompareCoeffMap", "EEMCompareCoeffMap", 100, 0, 100, 100, 0, 100);
0111   TH2F EEMCompareCoeffEtaTrend("EEMCompareCoeffEtaTrend", "EEMCompareCoeffEtaTrend", 51, 0, 50, 500, 0, 2);
0112   TProfile EEMCompareCoeffEtaProfile("EEMCompareCoeffEtaProfile", "EEMCompareCoeffEtaProfile", 51, 0, 50, 0, 2);
0113
0114   // ECAL endcap -

0115   for (int ix = 1; ix <= 100; ++ix)
0116     for (int iy = 1; iy <= 100; ++iy) {
0117       int rad = static_cast<int>(sqrt((ix - 50) * (ix - 50) + (iy - 50) * (iy - 50)));
0118       if (rad < EEradStart || rad > EEradEnd)
0119         continue;
0120       if (!EEDetId::validDetId(ix, iy, -1))
0121         continue;
0122       EEDetId det = EEDetId(ix, iy, -1, EEDetId::XYMODE);
0123       if (*(iEEscalibMap.find(det.rawId())) == 0)
0124         continue;
0125       double factor = *(iEEcalibMap.find(det.rawId())) / *(iEEscalibMap.find(det.rawId()));
0126       EEMCompareCoeffDistr.Fill(factor);
0127       EEMCompareCoeffMap.Fill(ix, iy, factor);
0128       EEMCompareCoeffEtaTrend.Fill(rad, factor);
0129       EEMCompareCoeffEtaProfile.Fill(rad, factor);
0130     }  // ECAL endcap -

0131
0132   std::string filename = "coeffcompareEE.root";
0133   TFile out(filename.c_str(), "recreate");
0134   EEMCompareCoeffMap.Write();
0135   EEMCompareCoeffDistr.Write();
0136   EEMCompareCoeffEtaTrend.Write();
0137   EEMCompareCoeffEtaProfile.Write();
0138   EEPCompareCoeffMap.Write();
0139   EEPCompareCoeffDistr.Write();
0140   EEPCompareCoeffEtaTrend.Write();
0141   EEPCompareCoeffEtaProfile.Write();
0142   out.Close();
0143 }