Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:32:11

0001 #include "Validation/CaloTowers/interface/CaloTowersValidation.h"
0002 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0003 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0004 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0005 
0006 CaloTowersValidation::CaloTowersValidation(edm::ParameterSet const& conf) {
0007   tok_calo_ = consumes<CaloTowerCollection>(conf.getUntrackedParameter<edm::InputTag>("CaloTowerCollectionLabel"));
0008   tok_evt_ = consumes<edm::HepMCProduct>(edm::InputTag("generatorSmeared"));
0009 
0010   // DQM ROOT output
0011   outputFile_ = conf.getUntrackedParameter<std::string>("outputFile", "myfile.root");
0012 
0013   hcalselector_ = conf.getUntrackedParameter<std::string>("hcalselector", "all");
0014 
0015   mc_ = conf.getUntrackedParameter<bool>("mc", true);
0016   useAllHistos_ = conf.getUntrackedParameter<bool>("useAllHistos", false);
0017 
0018   etaMin[0] = 0.;
0019   etaMax[0] = 1.4;
0020   etaMin[1] = 1.4;
0021   etaMax[1] = 2.9;
0022   etaMin[2] = 2.9;
0023   etaMax[2] = 5.2;
0024 
0025   isub = 0;
0026   if (hcalselector_ == "HB")
0027     isub = 1;
0028   if (hcalselector_ == "HE")
0029     isub = 2;
0030   if (hcalselector_ == "HF")
0031     isub = 3;
0032 
0033   imc = 1;
0034   if (!mc_)
0035     imc = 0;
0036 
0037   if (!outputFile_.empty()) {
0038     edm::LogInfo("OutputInfo") << " Hcal RecHit Task histograms will be saved to '" << outputFile_.c_str() << "'";
0039   } else {
0040     edm::LogInfo("OutputInfo") << " Hcal RecHit Task histograms will NOT be saved";
0041   }
0042 
0043   nevent = 0;
0044   // const char * sub = hcalselector_.c_str();
0045 }
0046 
0047 CaloTowersValidation::~CaloTowersValidation() {}
0048 
0049 void CaloTowersValidation::bookHistograms(DQMStore::IBooker& ibooker,
0050                                           edm::Run const& irun,
0051                                           edm::EventSetup const& isetup) {
0052   Char_t histo[100];
0053 
0054   ibooker.setCurrentFolder("CaloTowersV/CaloTowersTask");
0055 
0056   //These two histos are not drawn by our macros, but they are used
0057   //in the EndJob for norms and such so I am leaving them alone for now
0058   //-------------------------------------------------------------------------------------------
0059   sprintf(histo, "Ntowers_per_event_vs_ieta");
0060   Ntowers_vs_ieta = ibooker.book1D(histo, histo, 83, -41.5, 41.5);
0061 
0062   sprintf(histo, "CaloTowersTask_map_Nentries");
0063   mapEnergy_N = ibooker.book2D(histo, histo, 83, -41.5, 41.5, 72, 0., 72.);
0064   //-------------------------------------------------------------------------------------------
0065 
0066   //These the single pion scan histos
0067   //-------------------------------------------------------------------------------------------
0068   sprintf(histo, "emean_vs_ieta_E");
0069   emean_vs_ieta_E = ibooker.bookProfile(histo, histo, 83, -41.5, 41.5, -100., 2000., " ");
0070   sprintf(histo, "emean_vs_ieta_H");
0071   emean_vs_ieta_H = ibooker.bookProfile(histo, histo, 83, -41.5, 41.5, -100., 2000., " ");
0072   sprintf(histo, "emean_vs_ieta_EH");
0073   emean_vs_ieta_EH = ibooker.bookProfile(histo, histo, 83, -41.5, 41.5, -100., 2000., " ");
0074   //These are drawn
0075   sprintf(histo, "emean_vs_ieta_E1");
0076   emean_vs_ieta_E1 = ibooker.bookProfile(histo, histo, 83, -41.5, 41.5, -100., 2000., " ");
0077   sprintf(histo, "emean_vs_ieta_H1");
0078   emean_vs_ieta_H1 = ibooker.bookProfile(histo, histo, 83, -41.5, 41.5, -100., 2000., " ");
0079   sprintf(histo, "emean_vs_ieta_EH1");
0080   emean_vs_ieta_EH1 = ibooker.bookProfile(histo, histo, 83, -41.5, 41.5, -100., 2000., " ");
0081   //-------------------------------------------------------------------------------------------
0082 
0083   //Map energy histos are not drawn
0084   if (useAllHistos_) {
0085     sprintf(histo, "CaloTowersTask_map_energy_E");
0086     mapEnergy_E = ibooker.book2D(histo, histo, 83, -41.5, 41.5, 72, 0., 72.);
0087     sprintf(histo, "CaloTowersTask_map_energy_H");
0088     mapEnergy_H = ibooker.book2D(histo, histo, 83, -41.5, 41.5, 72, 0., 72.);
0089     sprintf(histo, "CaloTowersTask_map_energy_EH");
0090     mapEnergy_EH = ibooker.book2D(histo, histo, 83, -41.5, 41.5, 72, 0., 72.);
0091   }
0092 
0093   //All ECAL cell histos are used
0094   // XXX: ECAL 0-25 [0-26, 26 bins]   HCAL 0-4 [0-5, 5 bins]
0095   sprintf(histo, "number_of_bad_cells_Ecal_EB");
0096   numBadCellsEcal_EB = ibooker.book1D(histo, histo, 26, 0, 26);
0097   sprintf(histo, "number_of_bad_cells_Ecal_EE");
0098   numBadCellsEcal_EE = ibooker.book1D(histo, histo, 26, 0, 26);
0099   sprintf(histo, "number_of_recovered_cells_Ecal_EB");
0100   numRcvCellsEcal_EB = ibooker.book1D(histo, histo, 26, 0, 26);
0101   sprintf(histo, "number_of_recovered_cells_Ecal_EE");
0102   numRcvCellsEcal_EE = ibooker.book1D(histo, histo, 26, 0, 26);
0103   sprintf(histo, "number_of_problematic_cells_Ecal_EB");
0104   numPrbCellsEcal_EB = ibooker.book1D(histo, histo, 26, 0, 26);
0105   sprintf(histo, "number_of_problematic_cells_Ecal_EE");
0106   numPrbCellsEcal_EE = ibooker.book1D(histo, histo, 26, 0, 26);
0107 
0108   //Occupancy vs. ieta is drawn, occupancy map is needed to draw it
0109   sprintf(histo, "CaloTowersTask_map_occupancy");
0110   occupancy_map = ibooker.book2D(histo, histo, 83, -41.5, 41.5, 72, 0., 72.);
0111 
0112   sprintf(histo, "CaloTowersTask_occupancy_vs_ieta");
0113   occupancy_vs_ieta = ibooker.book1D(histo, histo, 83, -41.5, 41.5);
0114 
0115   if (isub == 1 || isub == 0) {
0116     //All cell histos are used
0117     sprintf(histo, "number_of_bad_cells_Hcal_HB");
0118     numBadCellsHcal_HB = ibooker.book1D(histo, histo, 5, 0, 5);
0119     sprintf(histo, "number_of_recovered_cells_Hcal_HB");
0120     numRcvCellsHcal_HB = ibooker.book1D(histo, histo, 5, 0, 5);
0121     sprintf(histo, "number_of_problematic_cells_Hcal_HB");
0122     numPrbCellsHcal_HB = ibooker.book1D(histo, histo, 5, 0, 5);
0123 
0124     //These are the five oldest CaloTower histos used: NTowers, E in HCAL/ECAL, MET and SET
0125     //-------------------------------------------------------------------------------------------
0126     sprintf(histo, "CaloTowersTask_energy_HCAL_HB");
0127     meEnergyHcal_HB = ibooker.book1D(histo, histo, 4100, -200, 8000);
0128 
0129     sprintf(histo, "CaloTowersTask_energy_ECAL_HB");
0130     meEnergyEcal_HB = ibooker.book1D(histo, histo, 3100, -200, 6000);
0131 
0132     sprintf(histo, "CaloTowersTask_number_of_fired_towers_HB");
0133     meNumFiredTowers_HB = ibooker.book1D(histo, histo, 1000, 0, 2000);
0134 
0135     sprintf(histo, "CaloTowersTask_MET_HB");
0136     MET_HB = ibooker.book1D(histo, histo, 3000, 0., 3000.);
0137 
0138     sprintf(histo, "CaloTowersTask_SET_HB");
0139     SET_HB = ibooker.book1D(histo, histo, 8000, 0., 8000.);
0140     //-------------------------------------------------------------------------------------------
0141 
0142     //Timing histos and profiles -- all six are necessary
0143     //-------------------------------------------------------------------------------------------
0144     sprintf(histo, "CaloTowersTask_EM_Timing_HB");
0145     emTiming_HB = ibooker.book1D(histo, histo, 110, -120., 100.);
0146 
0147     sprintf(histo, "CaloTowersTask_HAD_Timing_HB");
0148     hadTiming_HB = ibooker.book1D(histo, histo, 70, -48., 92.);
0149 
0150     //Energy-Timing histos are divided into low, medium and high to reduce memory usage
0151     //EM
0152     sprintf(histo, "CaloTowersTask_EM_Energy_Timing_Low_HB");
0153     emEnergyTiming_Low_HB = ibooker.book2D(histo, histo, 40, 0., 40., 110, -120., 100.);
0154 
0155     sprintf(histo, "CaloTowersTask_EM_Energy_Timing_HB");
0156     emEnergyTiming_HB = ibooker.book2D(histo, histo, 200, 0., 400., 110, -120., 100.);
0157 
0158     sprintf(histo, "CaloTowersTask_EM_Energy_Timing_High_HB");
0159     emEnergyTiming_High_HB = ibooker.book2D(histo, histo, 200, 0., 3000., 110, -120., 100.);
0160 
0161     sprintf(histo, "CaloTowersTask_EM_Energy_Timing_profile_Low_HB");
0162     emEnergyTiming_profile_Low_HB = ibooker.bookProfile(histo, histo, 40, 0., 40., 110, -120., 100.);
0163 
0164     sprintf(histo, "CaloTowersTask_EM_Energy_Timing_profile_HB");
0165     emEnergyTiming_profile_HB = ibooker.bookProfile(histo, histo, 200, 0., 400., 110, -120., 100.);
0166 
0167     sprintf(histo, "CaloTowersTask_EM_Energy_Timing_profile_High_HB");
0168     emEnergyTiming_profile_High_HB = ibooker.bookProfile(histo, histo, 200, 0., 3000., 110, -120., 100.);
0169 
0170     //HAD
0171     sprintf(histo, "CaloTowersTask_HAD_Energy_Timing_Low_HB");
0172     hadEnergyTiming_Low_HB = ibooker.book2D(histo, histo, 40, 0., 40., 70, -48., 92.);
0173 
0174     sprintf(histo, "CaloTowersTask_HAD_Energy_Timing_HB");
0175     hadEnergyTiming_HB = ibooker.book2D(histo, histo, 100, 0., 200., 70, -48., 92.);
0176 
0177     sprintf(histo, "CaloTowersTask_HAD_Energy_Timing_High_HB");
0178     hadEnergyTiming_High_HB = ibooker.book2D(histo, histo, 300, 0., 3000., 70, -48., 92.);
0179 
0180     sprintf(histo, "CaloTowersTask_HAD_Energy_Timing_profile_Low_HB");
0181     hadEnergyTiming_profile_Low_HB = ibooker.bookProfile(histo, histo, 40, 0., 40., 70, -48., 92.);
0182 
0183     sprintf(histo, "CaloTowersTask_HAD_Energy_Timing_profile_HB");
0184     hadEnergyTiming_profile_HB = ibooker.bookProfile(histo, histo, 100, 0., 200., 70, -48., 92.);
0185 
0186     sprintf(histo, "CaloTowersTask_HAD_Energy_Timing_profile_High_HB");
0187     hadEnergyTiming_profile_High_HB = ibooker.bookProfile(histo, histo, 300, 0., 3000., 70, -48., 92.);
0188     //-------------------------------------------------------------------------------------------
0189 
0190     //Everything else is not drawn
0191     if (useAllHistos_) {
0192       sprintf(histo, "CaloTowersTask_sum_of_energy_HCAL_vs_ECAL_HB");
0193       meEnergyHcalvsEcal_HB = ibooker.book2D(histo, histo, 500, 0., 500., 500, 0., 500.);
0194 
0195       sprintf(histo, "CaloTowersTask_energy_OUTER_HB");
0196       meEnergyHO_HB = ibooker.book1D(histo, histo, 1640, -200, 8000);
0197 
0198       sprintf(histo, "CaloTowersTask_energy_of_ECAL_component_of_tower_HB");
0199       meEnergyEcalTower_HB = ibooker.book1D(histo, histo, 440, -200, 2000);
0200 
0201       sprintf(histo, "CaloTowersTask_energy_of_HCAL_component_of_tower_HB");
0202       meEnergyHcalTower_HB = ibooker.book1D(histo, histo, 440, -200, 2000);
0203 
0204       sprintf(histo, "CaloTowersTask_energy_HcalPlusEcalPlusHO_HB");
0205       meTotEnergy_HB = ibooker.book1D(histo, histo, 400, 0., 2000.);
0206 
0207       sprintf(histo, "CaloTowersTask_map_energy_HB");
0208       mapEnergy_HB = ibooker.book2D(histo, histo, 83, -41.5, 41.5, 72, 0., 72.);
0209       sprintf(histo, "CaloTowersTask_map_energy_HCAL_HB");
0210       mapEnergyHcal_HB = ibooker.book2D(histo, histo, 83, -41.5, 41.5, 72, 0., 72.);
0211       sprintf(histo, "CaloTowersTask_map_energy_ECAL_HB");
0212       mapEnergyEcal_HB = ibooker.book2D(histo, histo, 83, -41.5, 41.5, 72, 0., 72.);
0213 
0214       sprintf(histo, "CaloTowersTask_phi_MET_HB");
0215       phiMET_HB = ibooker.book1D(histo, histo, 72, -3.1415926535898, 3.1415926535898);
0216     }
0217   }
0218 
0219   if (isub == 2 || isub == 0) {
0220     //All cell histos are used
0221     sprintf(histo, "number_of_bad_cells_Hcal_HE");
0222     numBadCellsHcal_HE = ibooker.book1D(histo, histo, 5, 0, 5);
0223     sprintf(histo, "number_of_recovered_cells_Hcal_HE");
0224     numRcvCellsHcal_HE = ibooker.book1D(histo, histo, 5, 0, 5);
0225     sprintf(histo, "number_of_problematic_cells_Hcal_HE");
0226     numPrbCellsHcal_HE = ibooker.book1D(histo, histo, 5, 0, 5);
0227 
0228     //These are the five oldest CaloTower histos used: NTowers, E in HCAL/ECAL, MET and SET
0229     //-------------------------------------------------------------------------------------------
0230     sprintf(histo, "CaloTowersTask_energy_HCAL_HE");
0231     meEnergyHcal_HE = ibooker.book1D(histo, histo, 1240, -200, 6000);
0232 
0233     sprintf(histo, "CaloTowersTask_energy_ECAL_HE");
0234     meEnergyEcal_HE = ibooker.book1D(histo, histo, 1240, -200, 6000);
0235 
0236     sprintf(histo, "CaloTowersTask_number_of_fired_towers_HE");
0237     meNumFiredTowers_HE = ibooker.book1D(histo, histo, 1000, 0, 2000);
0238 
0239     sprintf(histo, "CaloTowersTask_MET_HE");
0240     MET_HE = ibooker.book1D(histo, histo, 1000, 0., 1000.);
0241 
0242     sprintf(histo, "CaloTowersTask_SET_HE");
0243     SET_HE = ibooker.book1D(histo, histo, 2000, 0., 2000.);
0244     //-------------------------------------------------------------------------------------------
0245 
0246     //Timing histos and profiles -- all six are necessary
0247     //-------------------------------------------------------------------------------------------
0248     sprintf(histo, "CaloTowersTask_EM_Timing_HE");
0249     emTiming_HE = ibooker.book1D(histo, histo, 110, -120., 100.);
0250 
0251     sprintf(histo, "CaloTowersTask_HAD_Timing_HE");
0252     hadTiming_HE = ibooker.book1D(histo, histo, 70, -48., 92.);
0253 
0254     //Energy-Timing histos are divided into low and normal to reduce memory usage
0255     //EM
0256     sprintf(histo, "CaloTowersTask_EM_Energy_Timing_Low_HE");
0257     emEnergyTiming_Low_HE = ibooker.book2D(histo, histo, 160, 0., 160., 110, -120., 100.);
0258 
0259     sprintf(histo, "CaloTowersTask_EM_Energy_Timing_HE");
0260     emEnergyTiming_HE = ibooker.book2D(histo, histo, 200, 0., 800., 110, -120., 100.);
0261 
0262     sprintf(histo, "CaloTowersTask_EM_Energy_Timing_profile_Low_HE");
0263     emEnergyTiming_profile_Low_HE = ibooker.bookProfile(histo, histo, 160, 0., 160., 110, -120., 100.);
0264 
0265     sprintf(histo, "CaloTowersTask_EM_Energy_Timing_profile_HE");
0266     emEnergyTiming_profile_HE = ibooker.bookProfile(histo, histo, 200, 0., 800., 110, -120., 100.);
0267 
0268     //HAD
0269     sprintf(histo, "CaloTowersTask_HAD_Energy_Timing_Low_HE");
0270     hadEnergyTiming_Low_HE = ibooker.book2D(histo, histo, 160, 0., 160., 70, -48., 92.);
0271 
0272     sprintf(histo, "CaloTowersTask_HAD_Energy_Timing_HE");
0273     hadEnergyTiming_HE = ibooker.book2D(histo, histo, 200, 0., 800., 70, -48., 92.);
0274 
0275     sprintf(histo, "CaloTowersTask_HAD_Energy_Timing_profile_Low_HE");
0276     hadEnergyTiming_profile_Low_HE = ibooker.bookProfile(histo, histo, 160, 0., 160., 70, -48., 92.);
0277 
0278     sprintf(histo, "CaloTowersTask_HAD_Energy_Timing_profile_HE");
0279     hadEnergyTiming_profile_HE = ibooker.bookProfile(histo, histo, 200, 0., 800., 70, -48., 92.);
0280     //-------------------------------------------------------------------------------------------
0281 
0282     //Everything else is not drawn
0283     if (useAllHistos_) {
0284       sprintf(histo, "CaloTowersTask_sum_of_energy_HCAL_vs_ECAL_HE");
0285       meEnergyHcalvsEcal_HE = ibooker.book2D(histo, histo, 500, 0., 500., 500, 0., 500.);
0286 
0287       sprintf(histo, "CaloTowersTask_energy_OUTER_HE");
0288       meEnergyHO_HE = ibooker.book1D(histo, histo, 440, -200, 2000);
0289 
0290       sprintf(histo, "CaloTowersTask_energy_of_ECAL_component_of_tower_HE");
0291       meEnergyEcalTower_HE = ibooker.book1D(histo, histo, 1100, -200, 2000);
0292 
0293       sprintf(histo, "CaloTowersTask_energy_of_HCAL_component_of_tower_HE");
0294       meEnergyHcalTower_HE = ibooker.book1D(histo, histo, 1100, -200, 2000);
0295 
0296       sprintf(histo, "CaloTowersTask_energy_HcalPlusEcalPlusHO_HE");
0297       meTotEnergy_HE = ibooker.book1D(histo, histo, 400, 0., 2000.);
0298 
0299       sprintf(histo, "CaloTowersTask_map_energy_HE");
0300       mapEnergy_HE = ibooker.book2D(histo, histo, 83, -41.5, 41.5, 72, 0., 72.);
0301       sprintf(histo, "CaloTowersTask_map_energy_HCAL_HE");
0302       mapEnergyHcal_HE = ibooker.book2D(histo, histo, 83, -41.5, 41.5, 72, 0., 72.);
0303       sprintf(histo, "CaloTowersTask_map_energy_ECAL_HE");
0304       mapEnergyEcal_HE = ibooker.book2D(histo, histo, 83, -41.5, 41.5, 72, 0., 72.);
0305 
0306       sprintf(histo, "CaloTowersTask_phi_MET_HE");
0307       phiMET_HE = ibooker.book1D(histo, histo, 72, -3.1415926535898, 3.1415926535898);
0308     }
0309   }
0310 
0311   if (isub == 3 || isub == 0) {
0312     //All cell histos are used
0313     sprintf(histo, "number_of_bad_cells_Hcal_HF");
0314     numBadCellsHcal_HF = ibooker.book1D(histo, histo, 5, 0, 5);
0315     sprintf(histo, "number_of_recovered_cells_Hcal_HF");
0316     numRcvCellsHcal_HF = ibooker.book1D(histo, histo, 5, 0, 5);
0317     sprintf(histo, "number_of_problematic_cells_Hcal_HF");
0318     numPrbCellsHcal_HF = ibooker.book1D(histo, histo, 5, 0, 5);
0319 
0320     //These are the five oldest CaloTower histos used: NTowers, E in HCAL/ECAL, MET and SET
0321     //-------------------------------------------------------------------------------------------
0322     sprintf(histo, "CaloTowersTask_energy_HCAL_HF");
0323     meEnergyHcal_HF = ibooker.book1D(histo, histo, 4040, -200, 20000);
0324 
0325     sprintf(histo, "CaloTowersTask_energy_ECAL_HF");
0326     meEnergyEcal_HF = ibooker.book1D(histo, histo, 2440, -200, 12000);
0327 
0328     sprintf(histo, "CaloTowersTask_number_of_fired_towers_HF");
0329     meNumFiredTowers_HF = ibooker.book1D(histo, histo, 1000, 0, 2000);
0330 
0331     sprintf(histo, "CaloTowersTask_MET_HF");
0332     MET_HF = ibooker.book1D(histo, histo, 500, 0., 500.);
0333 
0334     sprintf(histo, "CaloTowersTask_SET_HF");
0335     SET_HF = ibooker.book1D(histo, histo, 2000, 0., 2000.);
0336     //-------------------------------------------------------------------------------------------
0337 
0338     //Timing histos and profiles -- all six are necessary
0339     //-------------------------------------------------------------------------------------------
0340     sprintf(histo, "CaloTowersTask_EM_Timing_HF");
0341     emTiming_HF = ibooker.book1D(histo, histo, 110, -120., 100.);
0342 
0343     sprintf(histo, "CaloTowersTask_HAD_Timing_HF");
0344     hadTiming_HF = ibooker.book1D(histo, histo, 70, -48., 92.);
0345 
0346     //EM
0347     sprintf(histo, "CaloTowersTask_EM_Energy_Timing_HF");
0348     emEnergyTiming_HF = ibooker.book2D(histo, histo, 150, 0., 300., 110, -120., 100.);
0349 
0350     sprintf(histo, "CaloTowersTask_EM_Energy_Timing_profile_HF");
0351     emEnergyTiming_profile_HF = ibooker.bookProfile(histo, histo, 150, 0., 300., 110, -120., 100.);
0352 
0353     //HAD (requires two different sets of histograms to lower RAM usage)
0354     sprintf(histo, "CaloTowersTask_HAD_Energy_Timing_Low_HF");
0355     hadEnergyTiming_Low_HF = ibooker.book2D(histo, histo, 40, 0., 40., 70, -48., 92.);
0356 
0357     sprintf(histo, "CaloTowersTask_HAD_Energy_Timing_HF");
0358     hadEnergyTiming_HF = ibooker.book2D(histo, histo, 200, 0., 600., 70, -48., 92.);
0359 
0360     sprintf(histo, "CaloTowersTask_HAD_Energy_Timing_profile_Low_HF");
0361     hadEnergyTiming_profile_Low_HF = ibooker.bookProfile(histo, histo, 40, 0., 40., 70, -48., 92.);
0362 
0363     sprintf(histo, "CaloTowersTask_HAD_Energy_Timing_profile_HF");
0364     hadEnergyTiming_profile_HF = ibooker.bookProfile(histo, histo, 200, 0., 600., 70, -48., 92.);
0365     //-------------------------------------------------------------------------------------------
0366 
0367     //Everything else is not drawn
0368     if (useAllHistos_) {
0369       sprintf(histo, "CaloTowersTask_sum_of_energy_HCAL_vs_ECAL_HF");
0370       meEnergyHcalvsEcal_HF = ibooker.book2D(histo, histo, 500, 0., 500., 500, 0., 500.);
0371 
0372       sprintf(histo, "CaloTowersTask_energy_OUTER_HF");
0373       meEnergyHO_HF = ibooker.book1D(histo, histo, 440, -200, 2000);
0374 
0375       sprintf(histo, "CaloTowersTask_energy_of_ECAL_component_of_tower_HF");
0376       meEnergyEcalTower_HF = ibooker.book1D(histo, histo, 440, -200, 2000);
0377 
0378       sprintf(histo, "CaloTowersTask_energy_of_HCAL_component_of_tower_HF");
0379       meEnergyHcalTower_HF = ibooker.book1D(histo, histo, 440, -200, 2000);
0380 
0381       sprintf(histo, "CaloTowersTask_energy_HcalPlusEcalPlusHO_HF");
0382       meTotEnergy_HF = ibooker.book1D(histo, histo, 400, 0., 2000.);
0383 
0384       sprintf(histo, "CaloTowersTask_map_energy_HF");
0385       mapEnergy_HF = ibooker.book2D(histo, histo, 83, -41.5, 41.5, 72, 0., 72.);
0386       sprintf(histo, "CaloTowersTask_map_energy_HCAL_HF");
0387       mapEnergyHcal_HF = ibooker.book2D(histo, histo, 83, -41.5, 41.5, 72, 0., 72.);
0388       sprintf(histo, "CaloTowersTask_map_energy_ECAL_HF");
0389       mapEnergyEcal_HF = ibooker.book2D(histo, histo, 83, -41.5, 41.5, 72, 0., 72.);
0390 
0391       sprintf(histo, "CaloTowersTask_phi_MET_HF");
0392       phiMET_HF = ibooker.book1D(histo, histo, 72, -3.1415926535898, 3.1415926535898);
0393     }
0394   }
0395 }
0396 
0397 void CaloTowersValidation::analyze(edm::Event const& event, edm::EventSetup const& c) {
0398   nevent++;
0399 
0400   // bool     MC = false; // UNUSED
0401   double phi_MC = 9999.;
0402   double eta_MC = 9999.;
0403 
0404   if (imc != 0) {
0405     edm::Handle<edm::HepMCProduct> evtMC;
0406     event.getByToken(tok_evt_, evtMC);  // generator in late 310_preX
0407     if (!evtMC.isValid()) {
0408       std::cout << "no HepMCProduct found" << std::endl;
0409     } else {
0410       // MC=true; // UNUSED
0411       //    std::cout << "*** source HepMCProduct found"<< std::endl;
0412     }
0413 
0414     // MC particle with highest pt is taken as a direction reference
0415     double maxPt = -99999.;
0416     int npart = 0;
0417     const HepMC::GenEvent* myGenEvent = evtMC->GetEvent();
0418     for (HepMC::GenEvent::particle_const_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end();
0419          ++p) {
0420       double phip = (*p)->momentum().phi();
0421       double etap = (*p)->momentum().eta();
0422       //    phi_MC = phip;
0423       //    eta_MC = etap;
0424       double pt = (*p)->momentum().perp();
0425       if (pt > maxPt) {
0426         npart++;
0427         maxPt = pt;
0428         phi_MC = phip;
0429         eta_MC = etap;
0430       }
0431     }
0432     //  std::cout << "*** Max pT = " << maxPt <<  std::endl;
0433   }
0434 
0435   edm::Handle<CaloTowerCollection> towers;
0436   event.getByToken(tok_calo_, towers);
0437   CaloTowerCollection::const_iterator cal;
0438 
0439   double met;
0440   double phimet;
0441 
0442   // ieta scan
0443   double partR = 0.3;
0444   double Rmin = 9999.;
0445   double Econe = 0.;
0446   double Hcone = 0.;
0447   double Ee1 = 0.;
0448   double Eh1 = 0.;
0449   double ieta_MC = 9999;
0450   double iphi_MC = 9999;
0451   //  double  etaM   = 9999.;
0452 
0453   // HB
0454   double sumEnergyHcal_HB = 0.;
0455   double sumEnergyEcal_HB = 0.;
0456   double sumEnergyHO_HB = 0.;
0457   Int_t numFiredTowers_HB = 0;
0458   double metx_HB = 0.;
0459   double mety_HB = 0.;
0460   double metz_HB = 0.;
0461   double sEt_HB = 0.;
0462   // HE
0463   double sumEnergyHcal_HE = 0.;
0464   double sumEnergyEcal_HE = 0.;
0465   double sumEnergyHO_HE = 0.;
0466   Int_t numFiredTowers_HE = 0;
0467   double metx_HE = 0.;
0468   double mety_HE = 0.;
0469   double metz_HE = 0.;
0470   double sEt_HE = 0.;
0471   // HF
0472   double sumEnergyHcal_HF = 0.;
0473   double sumEnergyEcal_HF = 0.;
0474   double sumEnergyHO_HF = 0.;
0475   Int_t numFiredTowers_HF = 0;
0476   double metx_HF = 0.;
0477   double mety_HF = 0.;
0478   double metz_HF = 0.;
0479   double sEt_HF = 0.;
0480 
0481   for (cal = towers->begin(); cal != towers->end(); ++cal) {
0482     double eE = cal->emEnergy();
0483     double eH = cal->hadEnergy();
0484     double eHO = cal->outerEnergy();
0485     double etaT = cal->eta();
0486     double phiT = cal->phi();
0487     double en = cal->energy();
0488     double etT = cal->et();
0489     double had_tm = cal->hcalTime();
0490     double em_tm = cal->ecalTime();
0491 
0492     int numBadEcalCells = cal->numBadEcalCells();
0493     int numRcvEcalCells = cal->numRecoveredEcalCells();
0494     int numPrbEcalCells = cal->numProblematicEcalCells();
0495 
0496     int numBadHcalCells = cal->numBadHcalCells();
0497     int numRcvHcalCells = cal->numRecoveredHcalCells();
0498     int numPrbHcalCells = cal->numProblematicHcalCells();
0499 
0500     math::RhoEtaPhiVector mom(cal->et(), cal->eta(), cal->phi());
0501     //  Vector mom  = cal->momentum();
0502 
0503     // cell properties
0504     CaloTowerDetId idT = cal->id();
0505     int ieta = idT.ieta();
0506     int iphi = idT.iphi();
0507 
0508     // ecal:  0 EcalBarrel  1 EcalEndcap
0509     // hcal:  0 hcalBarrel  1 HcalEndcap  2 HcalForward
0510     std::vector<int> inEcals(2), inHcals(3);
0511     unsigned int constitSize = cal->constituentsSize();
0512     for (unsigned int ic = 0; ic < constitSize; ic++) {
0513       DetId detId = cal->constituent(ic);
0514       if (detId.det() == DetId::Ecal) {
0515         if (detId.subdetId() == EcalBarrel)
0516           inEcals[0] = 1;
0517         else if (detId.subdetId() == EcalEndcap)
0518           inEcals[1] = 1;
0519       }
0520       if (detId.det() == DetId::Hcal) {
0521         if (HcalDetId(detId).subdet() == HcalBarrel)
0522           inHcals[0] = 1;
0523         else if (HcalDetId(detId).subdet() == HcalEndcap)
0524           inHcals[1] = 1;
0525         else if (HcalDetId(detId).subdet() == HcalForward)
0526           inHcals[2] = 1;
0527       }
0528     }
0529     //All cell histos are used
0530     if (inEcals[0]) {
0531       numBadCellsEcal_EB->Fill(numBadEcalCells);
0532       numRcvCellsEcal_EB->Fill(numRcvEcalCells);
0533       numPrbCellsEcal_EB->Fill(numPrbEcalCells);
0534     }
0535     if (inEcals[1]) {
0536       numBadCellsEcal_EE->Fill(numBadEcalCells);
0537       numRcvCellsEcal_EE->Fill(numRcvEcalCells);
0538       numPrbCellsEcal_EE->Fill(numPrbEcalCells);
0539     }
0540 
0541     if (imc != 0) {
0542       double r = dR(eta_MC, phi_MC, etaT, phiT);
0543 
0544       if (r < partR) {
0545         Econe += eE;
0546         Hcone += eH;
0547 
0548         // closest to MC
0549         if (r < Rmin) {
0550           if (fabs(eta_MC) < 3.0 && (ieta > 29 || ieta < -29)) {
0551             ;
0552           } else {
0553             Rmin = r;
0554             ieta_MC = ieta;
0555             iphi_MC = iphi;
0556             Ee1 = eE;
0557             Eh1 = eH;
0558           }
0559         }
0560       }
0561     }
0562 
0563     //Ntowers is used in EndJob, occupancy_map is used for occupancy vs ieta
0564     Ntowers_vs_ieta->Fill(double(ieta), 1.);
0565     occupancy_map->Fill(double(ieta), double(iphi));
0566 
0567     if ((isub == 0 || isub == 1) && (fabs(etaT) < etaMax[0] && fabs(etaT) >= etaMin[0])) {
0568       //All cell histos are used
0569       numBadCellsHcal_HB->Fill(numBadHcalCells);
0570       numRcvCellsHcal_HB->Fill(numRcvHcalCells);
0571       numPrbCellsHcal_HB->Fill(numPrbHcalCells);
0572 
0573       //Map energy histos are not used
0574       if (useAllHistos_) {
0575         mapEnergy_HB->Fill(double(ieta), double(iphi), en);
0576         mapEnergyHcal_HB->Fill(double(ieta), double(iphi), eH);
0577         mapEnergyEcal_HB->Fill(double(ieta), double(iphi), eE);
0578       }
0579       //      std::cout << " e_ecal = " << eE << std::endl;
0580 
0581       //  simple sums
0582       sumEnergyHcal_HB += eH;
0583       sumEnergyEcal_HB += eE;
0584       sumEnergyHO_HB += eHO;
0585 
0586       numFiredTowers_HB++;
0587 
0588       //Not used
0589       if (useAllHistos_) {
0590         meEnergyEcalTower_HB->Fill(eE);
0591         meEnergyHcalTower_HB->Fill(eH);
0592       }
0593 
0594       // MET, SET & phimet
0595       //  double  etT = cal->et();
0596       metx_HB += mom.x();
0597       mety_HB += mom.y();  //etT * sin(phiT);
0598       sEt_HB += etT;
0599 
0600       //Timing (all histos are used)
0601       emTiming_HB->Fill(em_tm);
0602       hadTiming_HB->Fill(had_tm);
0603 
0604       emEnergyTiming_Low_HB->Fill(eE, em_tm);
0605       emEnergyTiming_HB->Fill(eE, em_tm);
0606       emEnergyTiming_High_HB->Fill(eE, em_tm);
0607       emEnergyTiming_profile_Low_HB->Fill(eE, em_tm);
0608       emEnergyTiming_profile_HB->Fill(eE, em_tm);
0609       emEnergyTiming_profile_High_HB->Fill(eE, em_tm);
0610 
0611       hadEnergyTiming_Low_HB->Fill(eH, had_tm);
0612       hadEnergyTiming_HB->Fill(eH, had_tm);
0613       hadEnergyTiming_High_HB->Fill(eH, had_tm);
0614       hadEnergyTiming_profile_Low_HB->Fill(eH, had_tm);
0615       hadEnergyTiming_profile_HB->Fill(eH, had_tm);
0616       hadEnergyTiming_profile_High_HB->Fill(eH, had_tm);
0617     }
0618 
0619     if ((isub == 0 || isub == 2) && (fabs(etaT) < etaMax[1] && fabs(etaT) >= etaMin[1])) {
0620       //All cell histos are used
0621       numBadCellsHcal_HE->Fill(numBadHcalCells);
0622       numRcvCellsHcal_HE->Fill(numRcvHcalCells);
0623       numPrbCellsHcal_HE->Fill(numPrbHcalCells);
0624 
0625       //Map energy histos are not used
0626       if (useAllHistos_) {
0627         mapEnergy_HE->Fill(double(ieta), double(iphi), en);
0628         mapEnergyHcal_HE->Fill(double(ieta), double(iphi), eH);
0629         mapEnergyEcal_HE->Fill(double(ieta), double(iphi), eE);
0630       }
0631       //      std::cout << " e_ecal = " << eE << std::endl;
0632 
0633       //  simple sums
0634       sumEnergyHcal_HE += eH;
0635       sumEnergyEcal_HE += eE;
0636       sumEnergyHO_HE += eHO;
0637 
0638       numFiredTowers_HE++;
0639 
0640       //Not used
0641       if (useAllHistos_) {
0642         meEnergyEcalTower_HE->Fill(eE);
0643         meEnergyHcalTower_HE->Fill(eH);
0644       }
0645       // MET, SET & phimet
0646       //  double  etT = cal->et();
0647       metx_HE += mom.x();
0648       mety_HE += mom.y();  //etT * sin(phiT);
0649       sEt_HE += etT;
0650 
0651       //Timing (all histos are used)
0652       emTiming_HE->Fill(em_tm);
0653       hadTiming_HE->Fill(had_tm);
0654 
0655       emEnergyTiming_Low_HE->Fill(eE, em_tm);
0656       emEnergyTiming_HE->Fill(eE, em_tm);
0657       emEnergyTiming_profile_Low_HE->Fill(eE, em_tm);
0658       emEnergyTiming_profile_HE->Fill(eE, em_tm);
0659 
0660       hadEnergyTiming_Low_HE->Fill(eH, had_tm);
0661       hadEnergyTiming_HE->Fill(eH, had_tm);
0662       hadEnergyTiming_profile_Low_HE->Fill(eH, had_tm);
0663       hadEnergyTiming_profile_HE->Fill(eH, had_tm);
0664     }
0665 
0666     if ((isub == 0 || isub == 3) && (fabs(etaT) < etaMax[2] && fabs(etaT) >= etaMin[2])) {
0667       //All cell histos are used
0668       numBadCellsHcal_HF->Fill(numBadHcalCells);
0669       numRcvCellsHcal_HF->Fill(numRcvHcalCells);
0670       numPrbCellsHcal_HF->Fill(numPrbHcalCells);
0671 
0672       //Map energy histos are not used
0673       if (useAllHistos_) {
0674         mapEnergy_HF->Fill(double(ieta), double(iphi), en);
0675         mapEnergyHcal_HF->Fill(double(ieta), double(iphi), eH);
0676         mapEnergyEcal_HF->Fill(double(ieta), double(iphi), eE);
0677       }
0678       //      std::cout << " e_ecal = " << eE << std::endl;
0679 
0680       //  simple sums
0681       sumEnergyHcal_HF += eH;
0682       sumEnergyEcal_HF += eE;
0683       sumEnergyHO_HF += eHO;
0684 
0685       numFiredTowers_HF++;
0686 
0687       //Not used
0688       if (useAllHistos_) {
0689         meEnergyEcalTower_HF->Fill(eE);
0690         meEnergyHcalTower_HF->Fill(eH);
0691       }
0692       // MET, SET & phimet
0693       //  double  etT = cal->et();
0694       metx_HF += mom.x();
0695       mety_HF += mom.y();  //etT * sin(phiT);
0696       sEt_HF += etT;
0697 
0698       //Timing (all histos are used)
0699       emTiming_HF->Fill(em_tm);
0700       hadTiming_HF->Fill(had_tm);
0701       emEnergyTiming_HF->Fill(eE, em_tm);
0702       emEnergyTiming_profile_HF->Fill(eE, em_tm);
0703 
0704       hadEnergyTiming_Low_HF->Fill(eH, had_tm);
0705       hadEnergyTiming_HF->Fill(eH, had_tm);
0706       hadEnergyTiming_profile_Low_HF->Fill(eH, had_tm);
0707       hadEnergyTiming_profile_HF->Fill(eH, had_tm);
0708     }
0709 
0710   }  // end of Towers cycle
0711 
0712   //These are the six single pion histos
0713   emean_vs_ieta_E->Fill(double(ieta_MC), Econe);
0714   emean_vs_ieta_H->Fill(double(ieta_MC), Hcone);
0715   emean_vs_ieta_EH->Fill(double(ieta_MC), Econe + Hcone);
0716   emean_vs_ieta_E1->Fill(double(ieta_MC), Ee1);
0717   emean_vs_ieta_H1->Fill(double(ieta_MC), Eh1);
0718   emean_vs_ieta_EH1->Fill(double(ieta_MC), Ee1 + Eh1);
0719 
0720   //Map histos are not used except the last one in EndJob
0721   if (useAllHistos_) {
0722     mapEnergy_E->Fill(double(ieta_MC), double(iphi_MC), Ee1);
0723     mapEnergy_H->Fill(double(ieta_MC), double(iphi_MC), Eh1);
0724     mapEnergy_EH->Fill(double(ieta_MC), double(iphi_MC), Ee1 + Eh1);
0725   }
0726   mapEnergy_N->Fill(double(ieta_MC), double(iphi_MC), 1.);
0727 
0728   if (isub == 0 || isub == 1) {
0729     met = sqrt(metx_HB * metx_HB + mety_HB * mety_HB);
0730     Vector metv(metx_HB, mety_HB, metz_HB);
0731     phimet = metv.phi();
0732 
0733     //Five oldest drawn histos first; the rest are not used
0734     meEnergyHcal_HB->Fill(sumEnergyHcal_HB);
0735     meEnergyEcal_HB->Fill(sumEnergyEcal_HB);
0736     meNumFiredTowers_HB->Fill(numFiredTowers_HB);
0737     MET_HB->Fill(met);
0738     SET_HB->Fill(sEt_HB);
0739 
0740     if (useAllHistos_) {
0741       meEnergyHcalvsEcal_HB->Fill(sumEnergyEcal_HB, sumEnergyHcal_HB);
0742       meEnergyHO_HB->Fill(sumEnergyHO_HB);
0743       meTotEnergy_HB->Fill(sumEnergyHcal_HB + sumEnergyEcal_HB + sumEnergyHO_HB);
0744       phiMET_HB->Fill(phimet);
0745     }
0746   }
0747 
0748   if (isub == 0 || isub == 2) {
0749     met = sqrt(metx_HE * metx_HE + mety_HE * mety_HE);
0750     Vector metv(metx_HE, mety_HE, metz_HE);
0751     phimet = metv.phi();
0752 
0753     //Five oldest drawn histos first; the rest are not used
0754     meEnergyHcal_HE->Fill(sumEnergyHcal_HE);
0755     meEnergyEcal_HE->Fill(sumEnergyEcal_HE);
0756     meNumFiredTowers_HE->Fill(numFiredTowers_HE);
0757     MET_HE->Fill(met);
0758     SET_HE->Fill(sEt_HE);
0759 
0760     if (useAllHistos_) {
0761       meEnergyHcalvsEcal_HE->Fill(sumEnergyEcal_HE, sumEnergyHcal_HE);
0762       meEnergyHO_HE->Fill(sumEnergyHO_HE);
0763       meTotEnergy_HE->Fill(sumEnergyHcal_HE + sumEnergyEcal_HE + sumEnergyHO_HE);
0764       phiMET_HE->Fill(phimet);
0765     }
0766   }
0767 
0768   if (isub == 0 || isub == 3) {
0769     met = sqrt(metx_HF * metx_HF + mety_HF * mety_HF);
0770     Vector metv(metx_HF, mety_HF, metz_HF);
0771     phimet = metv.phi();
0772 
0773     //Five oldest drawn histos first; the rest are not used
0774     meEnergyHcal_HF->Fill(sumEnergyHcal_HF);
0775     meEnergyEcal_HF->Fill(sumEnergyEcal_HF);
0776     meNumFiredTowers_HF->Fill(numFiredTowers_HF);
0777     MET_HF->Fill(met);
0778     SET_HF->Fill(sEt_HF);
0779 
0780     if (useAllHistos_) {
0781       meEnergyHcalvsEcal_HF->Fill(sumEnergyEcal_HF, sumEnergyHcal_HF);
0782       meEnergyHO_HF->Fill(sumEnergyHO_HF);
0783       meTotEnergy_HF->Fill(sumEnergyHcal_HF + sumEnergyEcal_HF + sumEnergyHO_HF);
0784       phiMET_HF->Fill(phimet);
0785     }
0786   }
0787 }
0788 
0789 double CaloTowersValidation::dR(double eta1, double phi1, double eta2, double phi2) {
0790   double PI = 3.1415926535898;
0791   double deltaphi = phi1 - phi2;
0792   if (phi2 > phi1) {
0793     deltaphi = phi2 - phi1;
0794   }
0795   if (deltaphi > PI) {
0796     deltaphi = 2. * PI - deltaphi;
0797   }
0798   double deltaeta = eta2 - eta1;
0799   double tmp = sqrt(deltaeta * deltaeta + deltaphi * deltaphi);
0800   return tmp;
0801 }
0802 
0803 DEFINE_FWK_MODULE(CaloTowersValidation);