Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-08-13 23:26:56 

0001 #include "FWCore/Framework/interface/Event.h"
0002 #include "FWCore/Framework/interface/EventSetup.h"
0003 #include "FWCore/Framework/interface/Frameworkfwd.h"
0004 #include "FWCore/Framework/interface/MakerMacros.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007 #include "FWCore/ServiceRegistry/interface/Service.h"
0008 #include "FWCore/Utilities/interface/InputTag.h"
0009 
0010 // DQM include files
0011 
0012 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0013 #include "DQMServices/Core/interface/DQMStore.h"
0014 
0015 // work on collections
0016 #include "DataFormats/DetId/interface/DetId.h"
0017 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0018 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0019 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0020 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0021 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0022 #include "DataFormats/Math/interface/Point3D.h"
0023 #include "DataFormats/Common/interface/Handle.h"
0024 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidate.h"
0025 #include "DataFormats/RecoCandidate/interface/RecoEcalCandidateFwd.h"
0026 
0027 // Geometry
0028 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0029 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0030 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0031 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0032 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
0033 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0034 #include "Geometry/CaloTopology/interface/EcalBarrelTopology.h"
0035 #include "Geometry/CaloTopology/interface/EcalEndcapTopology.h"
0036 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0037 #include "RecoEcal/EgammaCoreTools/interface/PositionCalc.h"
0038 
0039 #include "TVector3.h"
0040 
0041 // Less than operator for sorting EcalRecHits according to energy.
0042 inline bool ecalRecHitGreater(EcalRecHit x, EcalRecHit y) { return (x.energy() > y.energy()); }
0043 
0044 class DQMSourcePi0 : public DQMEDAnalyzer {
0045 public:
0046   DQMSourcePi0(const edm::ParameterSet &);
0047   ~DQMSourcePi0() override;
0048 
0049 protected:
0050   void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override;
0051   void analyze(const edm::Event &e, const edm::EventSetup &c) override;
0052 
0053   void convxtalid(int &, int &);
0054   int diff_neta_s(int, int);
0055   int diff_nphi_s(int, int);
0056 
0057 private:
0058   int eventCounter_;
0059   PositionCalc posCalculator_;
0060 
0061   /// Distribution of rechits in iPhi (pi0)
0062   MonitorElement *hiPhiDistrEBpi0_;
0063 
0064   /// Distribution of rechits in ix EE  (pi0)
0065   MonitorElement *hiXDistrEEpi0_;
0066 
0067   /// Distribution of rechits in iPhi (eta)
0068   MonitorElement *hiPhiDistrEBeta_;
0069 
0070   /// Distribution of rechits in ix EE  (eta)
0071   MonitorElement *hiXDistrEEeta_;
0072 
0073   /// Distribution of rechits in iEta (pi0)
0074   MonitorElement *hiEtaDistrEBpi0_;
0075 
0076   /// Distribution of rechits in iy EE (pi0)
0077   MonitorElement *hiYDistrEEpi0_;
0078 
0079   /// Distribution of rechits in iEta (eta)
0080   MonitorElement *hiEtaDistrEBeta_;
0081 
0082   /// Distribution of rechits in iy EE (eta)
0083   MonitorElement *hiYDistrEEeta_;
0084 
0085   /// Energy Distribution of rechits EB (pi0)
0086   MonitorElement *hRechitEnergyEBpi0_;
0087 
0088   /// Energy Distribution of rechits EE (pi0)
0089   MonitorElement *hRechitEnergyEEpi0_;
0090 
0091   /// Energy Distribution of rechits EB (eta)
0092   MonitorElement *hRechitEnergyEBeta_;
0093 
0094   /// Energy Distribution of rechits EE (eta)
0095   MonitorElement *hRechitEnergyEEeta_;
0096 
0097   /// Distribution of total event energy EB (pi0)
0098   MonitorElement *hEventEnergyEBpi0_;
0099 
0100   /// Distribution of total event energy EE (pi0)
0101   MonitorElement *hEventEnergyEEpi0_;
0102 
0103   /// Distribution of total event energy EB (eta)
0104   MonitorElement *hEventEnergyEBeta_;
0105 
0106   /// Distribution of total event energy EE (eta)
0107   MonitorElement *hEventEnergyEEeta_;
0108 
0109   /// Distribution of number of RecHits EB (pi0)
0110   MonitorElement *hNRecHitsEBpi0_;
0111 
0112   /// Distribution of number of RecHits EE (pi0)
0113   MonitorElement *hNRecHitsEEpi0_;
0114 
0115   /// Distribution of number of RecHits EB (eta)
0116   MonitorElement *hNRecHitsEBeta_;
0117 
0118   /// Distribution of number of RecHits EE (eta)
0119   MonitorElement *hNRecHitsEEeta_;
0120 
0121   /// Distribution of Mean energy per rechit EB (pi0)
0122   MonitorElement *hMeanRecHitEnergyEBpi0_;
0123 
0124   /// Distribution of Mean energy per rechit EE (pi0)
0125   MonitorElement *hMeanRecHitEnergyEEpi0_;
0126 
0127   /// Distribution of Mean energy per rechit EB (eta)
0128   MonitorElement *hMeanRecHitEnergyEBeta_;
0129 
0130   /// Distribution of Mean energy per rechit EE (eta)
0131   MonitorElement *hMeanRecHitEnergyEEeta_;
0132 
0133   /// Pi0 invariant mass in EB
0134   MonitorElement *hMinvPi0EB_;
0135 
0136   /// Pi0 invariant mass in EE
0137   MonitorElement *hMinvPi0EE_;
0138 
0139   /// Eta invariant mass in EB
0140   MonitorElement *hMinvEtaEB_;
0141 
0142   /// Eta invariant mass in EE
0143   MonitorElement *hMinvEtaEE_;
0144 
0145   /// Pt of the 1st most energetic Pi0 photon in EB
0146   MonitorElement *hPt1Pi0EB_;
0147 
0148   /// Pt of the 1st most energetic Pi0 photon in EE
0149   MonitorElement *hPt1Pi0EE_;
0150 
0151   /// Pt of the 1st most energetic Eta photon in EB
0152   MonitorElement *hPt1EtaEB_;
0153 
0154   /// Pt of the 1st most energetic Eta photon in EE
0155   MonitorElement *hPt1EtaEE_;
0156 
0157   /// Pt of the 2nd most energetic Pi0 photon in EB
0158   MonitorElement *hPt2Pi0EB_;
0159 
0160   /// Pt of the 2nd most energetic Pi0 photon in EE
0161   MonitorElement *hPt2Pi0EE_;
0162 
0163   /// Pt of the 2nd most energetic Eta photon in EB
0164   MonitorElement *hPt2EtaEB_;
0165 
0166   /// Pt of the 2nd most energetic Eta photon in EE
0167   MonitorElement *hPt2EtaEE_;
0168 
0169   /// Pi0 Pt in EB
0170   MonitorElement *hPtPi0EB_;
0171 
0172   /// Pi0 Pt in EE
0173   MonitorElement *hPtPi0EE_;
0174 
0175   /// Eta Pt in EB
0176   MonitorElement *hPtEtaEB_;
0177 
0178   /// Eta Pt in EE
0179   MonitorElement *hPtEtaEE_;
0180 
0181   /// Pi0 Iso EB
0182   MonitorElement *hIsoPi0EB_;
0183 
0184   /// Pi0 Iso EE
0185   MonitorElement *hIsoPi0EE_;
0186 
0187   /// Eta Iso EB
0188   MonitorElement *hIsoEtaEB_;
0189 
0190   /// Eta Iso EE
0191   MonitorElement *hIsoEtaEE_;
0192 
0193   /// S4S9 of the 1st most energetic pi0 photon
0194   MonitorElement *hS4S91Pi0EB_;
0195 
0196   /// S4S9 of the 1st most energetic pi0 photon EE
0197   MonitorElement *hS4S91Pi0EE_;
0198 
0199   /// S4S9 of the 1st most energetic eta photon
0200   MonitorElement *hS4S91EtaEB_;
0201 
0202   /// S4S9 of the 1st most energetic eta photon EE
0203   MonitorElement *hS4S91EtaEE_;
0204 
0205   /// S4S9 of the 2nd most energetic pi0 photon
0206   MonitorElement *hS4S92Pi0EB_;
0207 
0208   /// S4S9 of the 2nd most energetic pi0 photon EE
0209   MonitorElement *hS4S92Pi0EE_;
0210 
0211   /// S4S9 of the 2nd most energetic eta photon
0212   MonitorElement *hS4S92EtaEB_;
0213 
0214   /// S4S9 of the 2nd most energetic eta photon EE
0215   MonitorElement *hS4S92EtaEE_;
0216 
0217   /// object to monitor
0218   edm::EDGetTokenT<EcalRecHitCollection> productMonitoredEBpi0_;
0219   edm::EDGetTokenT<EcalRecHitCollection> productMonitoredEBeta_;
0220 
0221   /// object to monitor
0222   edm::EDGetTokenT<EcalRecHitCollection> productMonitoredEEpi0_;
0223   edm::EDGetTokenT<EcalRecHitCollection> productMonitoredEEeta_;
0224 
0225   edm::ESGetToken<CaloTopology, CaloTopologyRecord> caloTopoToken_;
0226   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeomToken_;
0227 
0228   int gammaCandEtaSize_;
0229   int gammaCandPhiSize_;
0230 
0231   double seleXtalMinEnergy_;
0232   double seleXtalMinEnergyEndCap_;
0233 
0234   double clusSeedThr_;
0235   int clusEtaSize_;
0236   int clusPhiSize_;
0237 
0238   double clusSeedThrEndCap_;
0239 
0240   //// for pi0->gg barrel
0241   double selePtGamma_;
0242   double selePtPi0_;
0243   double seleMinvMaxPi0_;
0244   double seleMinvMinPi0_;
0245   double seleS4S9Gamma_;
0246   double selePi0BeltDR_;
0247   double selePi0BeltDeta_;
0248   double selePi0Iso_;
0249   double ptMinForIsolation_;
0250 
0251   /// for pi0->gg endcap
0252   double selePtGammaEndCap_;
0253   double selePtPi0EndCap_;
0254   double seleMinvMaxPi0EndCap_;
0255   double seleMinvMinPi0EndCap_;
0256   double seleS4S9GammaEndCap_;
0257   double selePi0IsoEndCap_;
0258   double selePi0BeltDREndCap_;
0259   double selePi0BeltDetaEndCap_;
0260   double ptMinForIsolationEndCap_;
0261 
0262   /// for eta->gg barrel
0263   double selePtGammaEta_;
0264   double selePtEta_;
0265   double seleS4S9GammaEta_;
0266   double seleS9S25GammaEta_;
0267   double seleMinvMaxEta_;
0268   double seleMinvMinEta_;
0269   double ptMinForIsolationEta_;
0270   double seleEtaIso_;
0271   double seleEtaBeltDR_;
0272   double seleEtaBeltDeta_;
0273 
0274   /// for eta->gg endcap
0275   double selePtGammaEtaEndCap_;
0276   double seleS4S9GammaEtaEndCap_;
0277   double seleS9S25GammaEtaEndCap_;
0278   double selePtEtaEndCap_;
0279   double seleMinvMaxEtaEndCap_;
0280   double seleMinvMinEtaEndCap_;
0281   double ptMinForIsolationEtaEndCap_;
0282   double seleEtaIsoEndCap_;
0283   double seleEtaBeltDREndCap_;
0284   double seleEtaBeltDetaEndCap_;
0285 
0286   bool ParameterLogWeighted_;
0287   double ParameterX0_;
0288   double ParameterT0_barl_;
0289   double ParameterT0_endc_;
0290   double ParameterT0_endcPresh_;
0291   double ParameterW0_;
0292 
0293   std::vector<EBDetId> detIdEBRecHits;
0294   std::vector<EcalRecHit> EBRecHits;
0295 
0296   std::vector<EEDetId> detIdEERecHits;
0297   std::vector<EcalRecHit> EERecHits;
0298 
0299   /// Monitor every prescaleFactor_ events
0300   unsigned int prescaleFactor_;
0301 
0302   /// DQM folder name
0303   std::string folderName_;
0304 
0305   /// Write to file
0306   bool saveToFile_;
0307 
0308   /// which subdet will be monitored
0309   bool isMonEBpi0_;
0310   bool isMonEBeta_;
0311   bool isMonEEpi0_;
0312   bool isMonEEeta_;
0313 
0314   /// Output file name if required
0315   std::string fileName_;
0316 };
0317 
0318 #define TWOPI 6.283185308
0319 
0320 // ******************************************
0321 // constructors
0322 // *****************************************
0323 
0324 DQMSourcePi0::DQMSourcePi0(const edm::ParameterSet &ps) : eventCounter_(0) {
0325   folderName_ = ps.getUntrackedParameter<std::string>("FolderName", "HLT/AlCaEcalPi0");
0326   prescaleFactor_ = ps.getUntrackedParameter<int>("prescaleFactor", 1);
0327   productMonitoredEBpi0_ =
0328       consumes<EcalRecHitCollection>(ps.getUntrackedParameter<edm::InputTag>("AlCaStreamEBpi0Tag"));
0329   productMonitoredEBeta_ =
0330       consumes<EcalRecHitCollection>(ps.getUntrackedParameter<edm::InputTag>("AlCaStreamEBetaTag"));
0331   productMonitoredEEpi0_ =
0332       consumes<EcalRecHitCollection>(ps.getUntrackedParameter<edm::InputTag>("AlCaStreamEEpi0Tag"));
0333   productMonitoredEEeta_ =
0334       consumes<EcalRecHitCollection>(ps.getUntrackedParameter<edm::InputTag>("AlCaStreamEEetaTag"));
0335   caloTopoToken_ = esConsumes();
0336   caloGeomToken_ = esConsumes();
0337 
0338   isMonEBpi0_ = ps.getUntrackedParameter<bool>("isMonEBpi0", false);
0339   isMonEBeta_ = ps.getUntrackedParameter<bool>("isMonEBeta", false);
0340   isMonEEpi0_ = ps.getUntrackedParameter<bool>("isMonEEpi0", false);
0341   isMonEEeta_ = ps.getUntrackedParameter<bool>("isMonEEeta", false);
0342 
0343   saveToFile_ = ps.getUntrackedParameter<bool>("SaveToFile", false);
0344   fileName_ = ps.getUntrackedParameter<std::string>("FileName", "MonitorAlCaEcalPi0.root");
0345 
0346   clusSeedThr_ = ps.getParameter<double>("clusSeedThr");
0347   clusSeedThrEndCap_ = ps.getParameter<double>("clusSeedThrEndCap");
0348   clusEtaSize_ = ps.getParameter<int>("clusEtaSize");
0349   clusPhiSize_ = ps.getParameter<int>("clusPhiSize");
0350   if (clusPhiSize_ % 2 == 0 || clusEtaSize_ % 2 == 0)
0351     edm::LogError("AlCaPi0RecHitsProducerError") << "Size of eta/phi for simple clustering should be odd numbers";
0352 
0353   seleXtalMinEnergy_ = ps.getParameter<double>("seleXtalMinEnergy");
0354   seleXtalMinEnergyEndCap_ = ps.getParameter<double>("seleXtalMinEnergyEndCap");
0355 
0356   /// for Pi0 barrel selection
0357   selePtGamma_ = ps.getParameter<double>("selePtGamma");
0358   selePtPi0_ = ps.getParameter<double>("selePtPi0");
0359   seleMinvMaxPi0_ = ps.getParameter<double>("seleMinvMaxPi0");
0360   seleMinvMinPi0_ = ps.getParameter<double>("seleMinvMinPi0");
0361   seleS4S9Gamma_ = ps.getParameter<double>("seleS4S9Gamma");
0362   selePi0Iso_ = ps.getParameter<double>("selePi0Iso");
0363   ptMinForIsolation_ = ps.getParameter<double>("ptMinForIsolation");
0364   selePi0BeltDR_ = ps.getParameter<double>("selePi0BeltDR");
0365   selePi0BeltDeta_ = ps.getParameter<double>("selePi0BeltDeta");
0366 
0367   /// for Pi0 endcap selection
0368   selePtGammaEndCap_ = ps.getParameter<double>("selePtGammaEndCap");
0369   selePtPi0EndCap_ = ps.getParameter<double>("selePtPi0EndCap");
0370   seleS4S9GammaEndCap_ = ps.getParameter<double>("seleS4S9GammaEndCap");
0371   seleMinvMaxPi0EndCap_ = ps.getParameter<double>("seleMinvMaxPi0EndCap");
0372   seleMinvMinPi0EndCap_ = ps.getParameter<double>("seleMinvMinPi0EndCap");
0373   ptMinForIsolationEndCap_ = ps.getParameter<double>("ptMinForIsolationEndCap");
0374   selePi0BeltDREndCap_ = ps.getParameter<double>("selePi0BeltDREndCap");
0375   selePi0BeltDetaEndCap_ = ps.getParameter<double>("selePi0BeltDetaEndCap");
0376   selePi0IsoEndCap_ = ps.getParameter<double>("selePi0IsoEndCap");
0377 
0378   /// for Eta barrel selection
0379   selePtGammaEta_ = ps.getParameter<double>("selePtGammaEta");
0380   selePtEta_ = ps.getParameter<double>("selePtEta");
0381   seleS4S9GammaEta_ = ps.getParameter<double>("seleS4S9GammaEta");
0382   seleS9S25GammaEta_ = ps.getParameter<double>("seleS9S25GammaEta");
0383   seleMinvMaxEta_ = ps.getParameter<double>("seleMinvMaxEta");
0384   seleMinvMinEta_ = ps.getParameter<double>("seleMinvMinEta");
0385   ptMinForIsolationEta_ = ps.getParameter<double>("ptMinForIsolationEta");
0386   seleEtaIso_ = ps.getParameter<double>("seleEtaIso");
0387   seleEtaBeltDR_ = ps.getParameter<double>("seleEtaBeltDR");
0388   seleEtaBeltDeta_ = ps.getParameter<double>("seleEtaBeltDeta");
0389 
0390   /// for Eta endcap selection
0391   selePtGammaEtaEndCap_ = ps.getParameter<double>("selePtGammaEtaEndCap");
0392   selePtEtaEndCap_ = ps.getParameter<double>("selePtEtaEndCap");
0393   seleS4S9GammaEtaEndCap_ = ps.getParameter<double>("seleS4S9GammaEtaEndCap");
0394   seleS9S25GammaEtaEndCap_ = ps.getParameter<double>("seleS9S25GammaEtaEndCap");
0395   seleMinvMaxEtaEndCap_ = ps.getParameter<double>("seleMinvMaxEtaEndCap");
0396   seleMinvMinEtaEndCap_ = ps.getParameter<double>("seleMinvMinEtaEndCap");
0397   ptMinForIsolationEtaEndCap_ = ps.getParameter<double>("ptMinForIsolationEtaEndCap");
0398   seleEtaIsoEndCap_ = ps.getParameter<double>("seleEtaIsoEndCap");
0399   seleEtaBeltDREndCap_ = ps.getParameter<double>("seleEtaBeltDREndCap");
0400   seleEtaBeltDetaEndCap_ = ps.getParameter<double>("seleEtaBeltDetaEndCap");
0401 
0402   // Parameters for the position calculation:
0403   edm::ParameterSet posCalcParameters = ps.getParameter<edm::ParameterSet>("posCalcParameters");
0404   posCalculator_ = PositionCalc(posCalcParameters);
0405 }
0406 
0407 DQMSourcePi0::~DQMSourcePi0() {}
0408 
0409 //--------------------------------------------------------
0410 void DQMSourcePi0::bookHistograms(DQMStore::IBooker &ibooker, edm::Run const &irun, edm::EventSetup const &isetup) {
0411   // create and cd into new folder
0412   ibooker.setCurrentFolder(folderName_);
0413 
0414   // book some histograms 1D
0415 
0416   hiPhiDistrEBpi0_ = ibooker.book1D("iphiDistributionEBpi0", "RechitEB pi0 iphi", 361, 1, 361);
0417   hiPhiDistrEBpi0_->setAxisTitle("i#phi ", 1);
0418   hiPhiDistrEBpi0_->setAxisTitle("# rechits", 2);
0419 
0420   hiXDistrEEpi0_ = ibooker.book1D("iXDistributionEEpi0", "RechitEE pi0 ix", 100, 0, 100);
0421   hiXDistrEEpi0_->setAxisTitle("ix ", 1);
0422   hiXDistrEEpi0_->setAxisTitle("# rechits", 2);
0423 
0424   hiPhiDistrEBeta_ = ibooker.book1D("iphiDistributionEBeta", "RechitEB eta iphi", 361, 1, 361);
0425   hiPhiDistrEBeta_->setAxisTitle("i#phi ", 1);
0426   hiPhiDistrEBeta_->setAxisTitle("# rechits", 2);
0427 
0428   hiXDistrEEeta_ = ibooker.book1D("iXDistributionEEeta", "RechitEE eta ix", 100, 0, 100);
0429   hiXDistrEEeta_->setAxisTitle("ix ", 1);
0430   hiXDistrEEeta_->setAxisTitle("# rechits", 2);
0431 
0432   hiEtaDistrEBpi0_ = ibooker.book1D("iEtaDistributionEBpi0", "RechitEB pi0 ieta", 171, -85, 86);
0433   hiEtaDistrEBpi0_->setAxisTitle("i#eta", 1);
0434   hiEtaDistrEBpi0_->setAxisTitle("#rechits", 2);
0435 
0436   hiYDistrEEpi0_ = ibooker.book1D("iYDistributionEEpi0", "RechitEE pi0 iY", 100, 0, 100);
0437   hiYDistrEEpi0_->setAxisTitle("iy", 1);
0438   hiYDistrEEpi0_->setAxisTitle("#rechits", 2);
0439 
0440   hiEtaDistrEBeta_ = ibooker.book1D("iEtaDistributionEBeta", "RechitEB eta ieta", 171, -85, 86);
0441   hiEtaDistrEBeta_->setAxisTitle("i#eta", 1);
0442   hiEtaDistrEBeta_->setAxisTitle("#rechits", 2);
0443 
0444   hiYDistrEEeta_ = ibooker.book1D("iYDistributionEEeta", "RechitEE eta iY", 100, 0, 100);
0445   hiYDistrEEeta_->setAxisTitle("iy", 1);
0446   hiYDistrEEeta_->setAxisTitle("#rechits", 2);
0447 
0448   hRechitEnergyEBpi0_ = ibooker.book1D("rhEnergyEBpi0", "Pi0 rechits energy EB", 160, 0., 2.0);
0449   hRechitEnergyEBpi0_->setAxisTitle("energy (GeV) ", 1);
0450   hRechitEnergyEBpi0_->setAxisTitle("#rechits", 2);
0451 
0452   hRechitEnergyEEpi0_ = ibooker.book1D("rhEnergyEEpi0", "Pi0 rechits energy EE", 160, 0., 3.0);
0453   hRechitEnergyEEpi0_->setAxisTitle("energy (GeV) ", 1);
0454   hRechitEnergyEEpi0_->setAxisTitle("#rechits", 2);
0455 
0456   hRechitEnergyEBeta_ = ibooker.book1D("rhEnergyEBeta", "Eta rechits energy EB", 160, 0., 2.0);
0457   hRechitEnergyEBeta_->setAxisTitle("energy (GeV) ", 1);
0458   hRechitEnergyEBeta_->setAxisTitle("#rechits", 2);
0459 
0460   hRechitEnergyEEeta_ = ibooker.book1D("rhEnergyEEeta", "Eta rechits energy EE", 160, 0., 3.0);
0461   hRechitEnergyEEeta_->setAxisTitle("energy (GeV) ", 1);
0462   hRechitEnergyEEeta_->setAxisTitle("#rechits", 2);
0463 
0464   hEventEnergyEBpi0_ = ibooker.book1D("eventEnergyEBpi0", "Pi0 event energy EB", 100, 0., 20.0);
0465   hEventEnergyEBpi0_->setAxisTitle("energy (GeV) ", 1);
0466 
0467   hEventEnergyEEpi0_ = ibooker.book1D("eventEnergyEEpi0", "Pi0 event energy EE", 100, 0., 50.0);
0468   hEventEnergyEEpi0_->setAxisTitle("energy (GeV) ", 1);
0469 
0470   hEventEnergyEBeta_ = ibooker.book1D("eventEnergyEBeta", "Eta event energy EB", 100, 0., 20.0);
0471   hEventEnergyEBeta_->setAxisTitle("energy (GeV) ", 1);
0472 
0473   hEventEnergyEEeta_ = ibooker.book1D("eventEnergyEEeta", "Eta event energy EE", 100, 0., 50.0);
0474   hEventEnergyEEeta_->setAxisTitle("energy (GeV) ", 1);
0475 
0476   hNRecHitsEBpi0_ = ibooker.book1D("nRechitsEBpi0", "#rechits in pi0 collection EB", 100, 0., 250.);
0477   hNRecHitsEBpi0_->setAxisTitle("rechits ", 1);
0478 
0479   hNRecHitsEEpi0_ = ibooker.book1D("nRechitsEEpi0", "#rechits in pi0 collection EE", 100, 0., 250.);
0480   hNRecHitsEEpi0_->setAxisTitle("rechits ", 1);
0481 
0482   hNRecHitsEBeta_ = ibooker.book1D("nRechitsEBeta", "#rechits in eta collection EB", 100, 0., 250.);
0483   hNRecHitsEBeta_->setAxisTitle("rechits ", 1);
0484 
0485   hNRecHitsEEeta_ = ibooker.book1D("nRechitsEEeta", "#rechits in eta collection EE", 100, 0., 250.);
0486   hNRecHitsEEeta_->setAxisTitle("rechits ", 1);
0487 
0488   hMeanRecHitEnergyEBpi0_ = ibooker.book1D("meanEnergyEBpi0", "Mean rechit energy in pi0 collection EB", 50, 0., 2.);
0489   hMeanRecHitEnergyEBpi0_->setAxisTitle("Mean Energy [GeV] ", 1);
0490 
0491   hMeanRecHitEnergyEEpi0_ = ibooker.book1D("meanEnergyEEpi0", "Mean rechit energy in pi0 collection EE", 100, 0., 5.);
0492   hMeanRecHitEnergyEEpi0_->setAxisTitle("Mean Energy [GeV] ", 1);
0493 
0494   hMeanRecHitEnergyEBeta_ = ibooker.book1D("meanEnergyEBeta", "Mean rechit energy in eta collection EB", 50, 0., 2.);
0495   hMeanRecHitEnergyEBeta_->setAxisTitle("Mean Energy [GeV] ", 1);
0496 
0497   hMeanRecHitEnergyEEeta_ = ibooker.book1D("meanEnergyEEeta", "Mean rechit energy in eta collection EE", 100, 0., 5.);
0498   hMeanRecHitEnergyEEeta_->setAxisTitle("Mean Energy [GeV] ", 1);
0499 
0500   hMinvPi0EB_ = ibooker.book1D("Pi0InvmassEB", "Pi0 Invariant Mass in EB", 100, 0., 0.5);
0501   hMinvPi0EB_->setAxisTitle("Inv Mass [GeV] ", 1);
0502 
0503   hMinvPi0EE_ = ibooker.book1D("Pi0InvmassEE", "Pi0 Invariant Mass in EE", 100, 0., 0.5);
0504   hMinvPi0EE_->setAxisTitle("Inv Mass [GeV] ", 1);
0505 
0506   hMinvEtaEB_ = ibooker.book1D("EtaInvmassEB", "Eta Invariant Mass in EB", 100, 0., 0.85);
0507   hMinvEtaEB_->setAxisTitle("Inv Mass [GeV] ", 1);
0508 
0509   hMinvEtaEE_ = ibooker.book1D("EtaInvmassEE", "Eta Invariant Mass in EE", 100, 0., 0.85);
0510   hMinvEtaEE_->setAxisTitle("Inv Mass [GeV] ", 1);
0511 
0512   hPt1Pi0EB_ = ibooker.book1D("Pt1Pi0EB", "Pt 1st most energetic Pi0 photon in EB", 100, 0., 20.);
0513   hPt1Pi0EB_->setAxisTitle("1st photon Pt [GeV] ", 1);
0514 
0515   hPt1Pi0EE_ = ibooker.book1D("Pt1Pi0EE", "Pt 1st most energetic Pi0 photon in EE", 100, 0., 20.);
0516   hPt1Pi0EE_->setAxisTitle("1st photon Pt [GeV] ", 1);
0517 
0518   hPt1EtaEB_ = ibooker.book1D("Pt1EtaEB", "Pt 1st most energetic Eta photon in EB", 100, 0., 20.);
0519   hPt1EtaEB_->setAxisTitle("1st photon Pt [GeV] ", 1);
0520 
0521   hPt1EtaEE_ = ibooker.book1D("Pt1EtaEE", "Pt 1st most energetic Eta photon in EE", 100, 0., 20.);
0522   hPt1EtaEE_->setAxisTitle("1st photon Pt [GeV] ", 1);
0523 
0524   hPt2Pi0EB_ = ibooker.book1D("Pt2Pi0EB", "Pt 2nd most energetic Pi0 photon in EB", 100, 0., 20.);
0525   hPt2Pi0EB_->setAxisTitle("2nd photon Pt [GeV] ", 1);
0526 
0527   hPt2Pi0EE_ = ibooker.book1D("Pt2Pi0EE", "Pt 2nd most energetic Pi0 photon in EE", 100, 0., 20.);
0528   hPt2Pi0EE_->setAxisTitle("2nd photon Pt [GeV] ", 1);
0529 
0530   hPt2EtaEB_ = ibooker.book1D("Pt2EtaEB", "Pt 2nd most energetic Eta photon in EB", 100, 0., 20.);
0531   hPt2EtaEB_->setAxisTitle("2nd photon Pt [GeV] ", 1);
0532 
0533   hPt2EtaEE_ = ibooker.book1D("Pt2EtaEE", "Pt 2nd most energetic Eta photon in EE", 100, 0., 20.);
0534   hPt2EtaEE_->setAxisTitle("2nd photon Pt [GeV] ", 1);
0535 
0536   hPtPi0EB_ = ibooker.book1D("PtPi0EB", "Pi0 Pt in EB", 100, 0., 20.);
0537   hPtPi0EB_->setAxisTitle("Pi0 Pt [GeV] ", 1);
0538 
0539   hPtPi0EE_ = ibooker.book1D("PtPi0EE", "Pi0 Pt in EE", 100, 0., 20.);
0540   hPtPi0EE_->setAxisTitle("Pi0 Pt [GeV] ", 1);
0541 
0542   hPtEtaEB_ = ibooker.book1D("PtEtaEB", "Eta Pt in EB", 100, 0., 20.);
0543   hPtEtaEB_->setAxisTitle("Eta Pt [GeV] ", 1);
0544 
0545   hPtEtaEE_ = ibooker.book1D("PtEtaEE", "Eta Pt in EE", 100, 0., 20.);
0546   hPtEtaEE_->setAxisTitle("Eta Pt [GeV] ", 1);
0547 
0548   hIsoPi0EB_ = ibooker.book1D("IsoPi0EB", "Pi0 Iso in EB", 50, 0., 1.);
0549   hIsoPi0EB_->setAxisTitle("Pi0 Iso", 1);
0550 
0551   hIsoPi0EE_ = ibooker.book1D("IsoPi0EE", "Pi0 Iso in EE", 50, 0., 1.);
0552   hIsoPi0EE_->setAxisTitle("Pi0 Iso", 1);
0553 
0554   hIsoEtaEB_ = ibooker.book1D("IsoEtaEB", "Eta Iso in EB", 50, 0., 1.);
0555   hIsoEtaEB_->setAxisTitle("Eta Iso", 1);
0556 
0557   hIsoEtaEE_ = ibooker.book1D("IsoEtaEE", "Eta Iso in EE", 50, 0., 1.);
0558   hIsoEtaEE_->setAxisTitle("Eta Iso", 1);
0559 
0560   hS4S91Pi0EB_ = ibooker.book1D("S4S91Pi0EB", "S4S9 1st most energetic Pi0 photon in EB", 50, 0., 1.);
0561   hS4S91Pi0EB_->setAxisTitle("S4S9 of the 1st Pi0 Photon ", 1);
0562 
0563   hS4S91Pi0EE_ = ibooker.book1D("S4S91Pi0EE", "S4S9 1st most energetic Pi0 photon in EE", 50, 0., 1.);
0564   hS4S91Pi0EE_->setAxisTitle("S4S9 of the 1st Pi0 Photon ", 1);
0565 
0566   hS4S91EtaEB_ = ibooker.book1D("S4S91EtaEB", "S4S9 1st most energetic Eta photon in EB", 50, 0., 1.);
0567   hS4S91EtaEB_->setAxisTitle("S4S9 of the 1st Eta Photon ", 1);
0568 
0569   hS4S91EtaEE_ = ibooker.book1D("S4S91EtaEE", "S4S9 1st most energetic Eta photon in EE", 50, 0., 1.);
0570   hS4S91EtaEE_->setAxisTitle("S4S9 of the 1st Eta Photon ", 1);
0571 
0572   hS4S92Pi0EB_ = ibooker.book1D("S4S92Pi0EB", "S4S9 2nd most energetic Pi0 photon in EB", 50, 0., 1.);
0573   hS4S92Pi0EB_->setAxisTitle("S4S9 of the 2nd Pi0 Photon", 1);
0574 
0575   hS4S92Pi0EE_ = ibooker.book1D("S4S92Pi0EE", "S4S9 2nd most energetic Pi0 photon in EE", 50, 0., 1.);
0576   hS4S92Pi0EE_->setAxisTitle("S4S9 of the 2nd Pi0 Photon", 1);
0577 
0578   hS4S92EtaEB_ = ibooker.book1D("S4S92EtaEB", "S4S9 2nd most energetic Pi0 photon in EB", 50, 0., 1.);
0579   hS4S92EtaEB_->setAxisTitle("S4S9 of the 2nd Eta Photon", 1);
0580 
0581   hS4S92EtaEE_ = ibooker.book1D("S4S92EtaEE", "S4S9 2nd most energetic Pi0 photon in EE", 50, 0., 1.);
0582   hS4S92EtaEE_->setAxisTitle("S4S9 of the 2nd Eta Photon", 1);
0583 }
0584 
0585 //-------------------------------------------------------------
0586 void DQMSourcePi0::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0587   if (eventCounter_ % prescaleFactor_)
0588     return;
0589   eventCounter_++;
0590 
0591   auto const &theCaloTopology = iSetup.getHandle(caloTopoToken_);
0592 
0593   std::vector<EcalRecHit> seeds;
0594   seeds.clear();
0595 
0596   std::vector<EBDetId> usedXtals;
0597   usedXtals.clear();
0598 
0599   detIdEBRecHits.clear();  //// EBDetId
0600   EBRecHits.clear();       /// EcalRecHit
0601 
0602   edm::Handle<EcalRecHitCollection> rhEBpi0;
0603   edm::Handle<EcalRecHitCollection> rhEBeta;
0604   edm::Handle<EcalRecHitCollection> rhEEpi0;
0605   edm::Handle<EcalRecHitCollection> rhEEeta;
0606 
0607   if (isMonEBpi0_)
0608     iEvent.getByToken(productMonitoredEBpi0_, rhEBpi0);
0609   if (isMonEBeta_)
0610     iEvent.getByToken(productMonitoredEBeta_, rhEBeta);
0611   if (isMonEEpi0_)
0612     iEvent.getByToken(productMonitoredEEpi0_, rhEEpi0);
0613   if (isMonEEeta_)
0614     iEvent.getByToken(productMonitoredEEeta_, rhEEeta);
0615 
0616   // Initialize the Position Calc
0617 
0618   //edm::ESHandle<CaloGeometry> geoHandle;
0619   //iSetup.get<CaloGeometryRecord>().get(geoHandle);
0620   const auto &geoHandle = iSetup.getHandle(caloGeomToken_);
0621   const CaloSubdetectorGeometry *geometry_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0622   const CaloSubdetectorGeometry *geometryEE_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
0623   const CaloSubdetectorGeometry *geometryES_p = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
0624 
0625   const CaloSubdetectorTopology *topology_p = theCaloTopology->getSubdetectorTopology(DetId::Ecal, EcalBarrel);
0626   const CaloSubdetectorTopology *topology_ee = theCaloTopology->getSubdetectorTopology(DetId::Ecal, EcalEndcap);
0627 
0628   EcalRecHitCollection::const_iterator itb;
0629 
0630   // fill EB pi0 histos
0631   if (isMonEBpi0_) {
0632     if (rhEBpi0.isValid() && (!rhEBpi0->empty())) {
0633       const EcalRecHitCollection *hitCollection_p = rhEBpi0.product();
0634       float etot = 0;
0635       for (itb = rhEBpi0->begin(); itb != rhEBpi0->end(); ++itb) {
0636         EBDetId id(itb->id());
0637         double energy = itb->energy();
0638         if (energy < seleXtalMinEnergy_)
0639           continue;
0640 
0641         EBDetId det = itb->id();
0642 
0643         detIdEBRecHits.push_back(det);
0644         EBRecHits.push_back(*itb);
0645 
0646         if (energy > clusSeedThr_)
0647           seeds.push_back(*itb);
0648 
0649         hiPhiDistrEBpi0_->Fill(id.iphi());
0650         hiEtaDistrEBpi0_->Fill(id.ieta());
0651         hRechitEnergyEBpi0_->Fill(itb->energy());
0652 
0653         etot += itb->energy();
0654       }  // Eb rechits
0655 
0656       hNRecHitsEBpi0_->Fill(rhEBpi0->size());
0657       hMeanRecHitEnergyEBpi0_->Fill(etot / rhEBpi0->size());
0658       hEventEnergyEBpi0_->Fill(etot);
0659 
0660       //      cout << " EB RH Pi0 collection: #, mean rh_e, event E
0661       //      "<<rhEBpi0->size()<<" "<<etot/rhEBpi0->size()<<" "<<etot<<endl;
0662 
0663       // Pi0 maker
0664 
0665       // cout<< " RH coll size: "<<rhEBpi0->size()<<endl;
0666       // cout<< " Pi0 seeds: "<<seeds.size()<<endl;
0667 
0668       int nClus;
0669       std::vector<float> eClus;
0670       std::vector<float> etClus;
0671       std::vector<float> etaClus;
0672       std::vector<float> thetaClus;
0673       std::vector<float> phiClus;
0674       std::vector<EBDetId> max_hit;
0675 
0676       std::vector<std::vector<EcalRecHit>> RecHitsCluster;
0677       std::vector<std::vector<EcalRecHit>> RecHitsCluster5x5;
0678       std::vector<float> s4s9Clus;
0679       std::vector<float> s9s25Clus;
0680 
0681       nClus = 0;
0682 
0683       // Make own simple clusters (3x3, 5x5 or clusPhiSize_ x clusEtaSize_)
0684       std::sort(seeds.begin(), seeds.end(), ecalRecHitGreater);
0685 
0686       for (std::vector<EcalRecHit>::iterator itseed = seeds.begin(); itseed != seeds.end(); itseed++) {
0687         EBDetId seed_id = itseed->id();
0688         std::vector<EBDetId>::const_iterator usedIds;
0689 
0690         bool seedAlreadyUsed = false;
0691         for (usedIds = usedXtals.begin(); usedIds != usedXtals.end(); usedIds++) {
0692           if (*usedIds == seed_id) {
0693             seedAlreadyUsed = true;
0694             // cout<< " Seed with energy "<<itseed->energy()<<" was used
0695             // !"<<endl;
0696             break;
0697           }
0698         }
0699         if (seedAlreadyUsed)
0700           continue;
0701         std::vector<DetId> clus_v = topology_p->getWindow(seed_id, clusEtaSize_, clusPhiSize_);
0702         std::vector<std::pair<DetId, float>> clus_used;
0703         // Reject the seed if not able to build the cluster around it correctly
0704         // if(clus_v.size() < clusEtaSize_*clusPhiSize_){cout<<" Not enough
0705         // RecHits "<<endl; continue;}
0706         std::vector<EcalRecHit> RecHitsInWindow;
0707         std::vector<EcalRecHit> RecHitsInWindow5x5;
0708 
0709         double simple_energy = 0;
0710 
0711         for (std::vector<DetId>::iterator det = clus_v.begin(); det != clus_v.end(); det++) {
0712           EBDetId EBdet = *det;
0713           //      cout<<" det "<< EBdet<<" ieta "<<EBdet.ieta()<<" iphi
0714           //      "<<EBdet.iphi()<<endl;
0715           bool HitAlreadyUsed = false;
0716           for (usedIds = usedXtals.begin(); usedIds != usedXtals.end(); usedIds++) {
0717             if (*usedIds == *det) {
0718               HitAlreadyUsed = true;
0719               break;
0720             }
0721           }
0722           if (HitAlreadyUsed)
0723             continue;
0724 
0725           std::vector<EBDetId>::iterator itdet = find(detIdEBRecHits.begin(), detIdEBRecHits.end(), EBdet);
0726           if (itdet == detIdEBRecHits.end())
0727             continue;
0728 
0729           int nn = int(itdet - detIdEBRecHits.begin());
0730           usedXtals.push_back(*det);
0731           RecHitsInWindow.push_back(EBRecHits[nn]);
0732           clus_used.push_back(std::make_pair(*det, 1));
0733           simple_energy = simple_energy + EBRecHits[nn].energy();
0734         }
0735 
0736         if (simple_energy <= 0)
0737           continue;
0738 
0739         math::XYZPoint clus_pos =
0740             posCalculator_.Calculate_Location(clus_used, hitCollection_p, geometry_p, geometryES_p);
0741         // cout<< "       Simple Clustering: Total energy for this simple
0742         // cluster : "<<simple_energy<<endl; cout<< "       Simple Clustering:
0743         // eta phi : "<<clus_pos.eta()<<" "<<clus_pos.phi()<<endl; cout<< "
0744         // Simple Clustering: x y z : "<<clus_pos.x()<<" "<<clus_pos.y()<<"
0745         // "<<clus_pos.z()<<endl;
0746 
0747         float theta_s = 2. * atan(exp(-clus_pos.eta()));
0748         //  float p0x_s = simple_energy * sin(theta_s) *
0749         // cos(clus_pos.phi()); float p0y_s = simple_energy * sin(theta_s) *
0750         // sin(clus_pos.phi());
0751         //    float p0z_s = simple_energy * cos(theta_s);
0752         // float et_s = sqrt( p0x_s*p0x_s + p0y_s*p0y_s);
0753 
0754         float et_s = simple_energy * sin(theta_s);
0755         // cout << "       Simple Clustering: E,Et,px,py,pz: "<<simple_energy<<"
0756         // "<<et_s<<" "<<p0x_s<<" "<<p0y_s<<" "<<endl;
0757 
0758         // Compute S4/S9 variable
0759         // We are not sure to have 9 RecHits so need to check eta and phi:
0760         /// check s4s9
0761         float s4s9_tmp[4];
0762         for (int i = 0; i < 4; i++)
0763           s4s9_tmp[i] = 0;
0764 
0765         int seed_ieta = seed_id.ieta();
0766         int seed_iphi = seed_id.iphi();
0767 
0768         convxtalid(seed_iphi, seed_ieta);
0769 
0770         float e3x3 = 0;
0771         float e5x5 = 0;
0772 
0773         for (unsigned int j = 0; j < RecHitsInWindow.size(); j++) {
0774           EBDetId det = (EBDetId)RecHitsInWindow[j].id();
0775 
0776           int ieta = det.ieta();
0777           int iphi = det.iphi();
0778 
0779           convxtalid(iphi, ieta);
0780 
0781           float en = RecHitsInWindow[j].energy();
0782 
0783           int dx = diff_neta_s(seed_ieta, ieta);
0784           int dy = diff_nphi_s(seed_iphi, iphi);
0785 
0786           if (dx <= 0 && dy <= 0)
0787             s4s9_tmp[0] += en;
0788           if (dx >= 0 && dy <= 0)
0789             s4s9_tmp[1] += en;
0790           if (dx <= 0 && dy >= 0)
0791             s4s9_tmp[2] += en;
0792           if (dx >= 0 && dy >= 0)
0793             s4s9_tmp[3] += en;
0794 
0795           if (std::abs(dx) <= 1 && std::abs(dy) <= 1)
0796             e3x3 += en;
0797           if (std::abs(dx) <= 2 && std::abs(dy) <= 2)
0798             e5x5 += en;
0799         }
0800 
0801         if (e3x3 <= 0)
0802           continue;
0803 
0804         float s4s9_max = *std::max_element(s4s9_tmp, s4s9_tmp + 4) / e3x3;
0805 
0806         /// calculate e5x5
0807         std::vector<DetId> clus_v5x5 = topology_p->getWindow(seed_id, 5, 5);
0808         for (std::vector<DetId>::const_iterator idItr = clus_v5x5.begin(); idItr != clus_v5x5.end(); idItr++) {
0809           EBDetId det = *idItr;
0810 
0811           std::vector<EBDetId>::iterator itdet0 = find(usedXtals.begin(), usedXtals.end(), det);
0812 
0813           /// already clustered
0814           if (itdet0 != usedXtals.end())
0815             continue;
0816 
0817           // inside collections
0818           std::vector<EBDetId>::iterator itdet = find(detIdEBRecHits.begin(), detIdEBRecHits.end(), det);
0819           if (itdet == detIdEBRecHits.end())
0820             continue;
0821 
0822           int nn = int(itdet - detIdEBRecHits.begin());
0823 
0824           RecHitsInWindow5x5.push_back(EBRecHits[nn]);
0825           e5x5 += EBRecHits[nn].energy();
0826         }
0827 
0828         if (e5x5 <= 0)
0829           continue;
0830 
0831         eClus.push_back(simple_energy);
0832         etClus.push_back(et_s);
0833         etaClus.push_back(clus_pos.eta());
0834         thetaClus.push_back(theta_s);
0835         phiClus.push_back(clus_pos.phi());
0836         s4s9Clus.push_back(s4s9_max);
0837         s9s25Clus.push_back(e3x3 / e5x5);
0838         RecHitsCluster.push_back(RecHitsInWindow);
0839         RecHitsCluster5x5.push_back(RecHitsInWindow5x5);
0840 
0841         //      std::cout<<" EB pi0 cluster (n,nxt,e,et eta,phi,s4s9)
0842         //"<<nClus<<" "<<int(RecHitsInWindow.size())<<" "<<eClus[nClus]<<" "<<"
0843         //"<<etClus[nClus]<<" "<<etaClus[nClus]<<" "<<phiClus[nClus]<<"
0844         //"<<s4s9Clus[nClus]<<std::endl;
0845 
0846         nClus++;
0847       }
0848 
0849       // cout<< " Pi0 clusters: "<<nClus<<endl;
0850 
0851       // Selection, based on Simple clustering
0852       // pi0 candidates
0853       int npi0_s = 0;
0854 
0855       //      if (nClus <= 1) return;
0856       for (Int_t i = 0; i < nClus; i++) {
0857         for (Int_t j = i + 1; j < nClus; j++) {
0858           //      cout<<" i "<<i<<"  etClus[i] "<<etClus[i]<<" j "<<j<<"
0859           //      etClus[j] "<<etClus[j]<<endl;
0860           if (etClus[i] > selePtGamma_ && etClus[j] > selePtGamma_ && s4s9Clus[i] > seleS4S9Gamma_ &&
0861               s4s9Clus[j] > seleS4S9Gamma_) {
0862             float p0x = etClus[i] * cos(phiClus[i]);
0863             float p1x = etClus[j] * cos(phiClus[j]);
0864             float p0y = etClus[i] * sin(phiClus[i]);
0865             float p1y = etClus[j] * sin(phiClus[j]);
0866             float p0z = eClus[i] * cos(thetaClus[i]);
0867             float p1z = eClus[j] * cos(thetaClus[j]);
0868 
0869             float pt_pair = sqrt((p0x + p1x) * (p0x + p1x) + (p0y + p1y) * (p0y + p1y));
0870 
0871             if (pt_pair < selePtPi0_)
0872               continue;
0873 
0874             float m_inv = sqrt((eClus[i] + eClus[j]) * (eClus[i] + eClus[j]) - (p0x + p1x) * (p0x + p1x) -
0875                                (p0y + p1y) * (p0y + p1y) - (p0z + p1z) * (p0z + p1z));
0876             if ((m_inv < seleMinvMaxPi0_) && (m_inv > seleMinvMinPi0_)) {
0877               // New Loop on cluster to measure isolation:
0878               std::vector<int> IsoClus;
0879               IsoClus.clear();
0880               float Iso = 0;
0881               TVector3 pairVect = TVector3((p0x + p1x), (p0y + p1y), (p0z + p1z));
0882               for (Int_t k = 0; k < nClus; k++) {
0883                 if (etClus[k] < ptMinForIsolation_)
0884                   continue;
0885 
0886                 if (k == i || k == j)
0887                   continue;
0888                 TVector3 ClusVect =
0889                     TVector3(etClus[k] * cos(phiClus[k]), etClus[k] * sin(phiClus[k]), eClus[k] * cos(thetaClus[k]));
0890 
0891                 float dretacl = fabs(etaClus[k] - pairVect.Eta());
0892                 float drcl = ClusVect.DeltaR(pairVect);
0893                 //      cout<< "   Iso: k, E, drclpi0, detaclpi0, dphiclpi0
0894                 //      "<<k<<" "<<eClus[k]<<" "<<drclpi0<<"
0895                 //      "<<dretaclpi0<<endl;
0896                 if ((drcl < selePi0BeltDR_) && (dretacl < selePi0BeltDeta_)) {
0897                   //              cout<< "   ... good iso cluster #: "<<k<<"
0898                   //              etClus[k] "<<etClus[k] <<endl;
0899                   Iso = Iso + etClus[k];
0900                   IsoClus.push_back(k);
0901                 }
0902               }
0903 
0904               //      cout<<"  Iso/pt_pi0 "<<Iso/pt_pi0<<endl;
0905               if (Iso / pt_pair < selePi0Iso_) {
0906                 // for(unsigned int
0907                 // Rec=0;Rec<RecHitsCluster[i].size();Rec++)pi0EBRecHitCollection->push_back(RecHitsCluster[i][Rec]);
0908                 // for(unsigned int
0909                 // Rec2=0;Rec2<RecHitsCluster[j].size();Rec2++)pi0EBRecHitCollection->push_back(RecHitsCluster[j][Rec2]);
0910 
0911                 hMinvPi0EB_->Fill(m_inv);
0912                 hPt1Pi0EB_->Fill(etClus[i]);
0913                 hPt2Pi0EB_->Fill(etClus[j]);
0914                 hPtPi0EB_->Fill(pt_pair);
0915                 hIsoPi0EB_->Fill(Iso / pt_pair);
0916                 hS4S91Pi0EB_->Fill(s4s9Clus[i]);
0917                 hS4S92Pi0EB_->Fill(s4s9Clus[j]);
0918 
0919                 //      cout <<"  EB Simple Clustering: pi0 Candidate
0920                 // pt, eta, phi, Iso, m_inv, i, j :   "<<pt_pair<<"
0921                 //"<<pairVect.Eta()<<" "<<pairVect.Phi()<<" "<<Iso<<"
0922                 //"<<m_inv<<" "<<i<<" "<<j<<" "<<endl;
0923 
0924                 npi0_s++;
0925               }
0926             }
0927           }
0928         }  // End of the "j" loop over Simple Clusters
0929       }    // End of the "i" loop over Simple Clusters
0930 
0931       //      cout<<"  (Simple Clustering) EB Pi0 candidates #: "<<npi0_s<<endl;
0932 
0933     }  // rhEBpi0.valid() ends
0934 
0935   }  //  isMonEBpi0 ends
0936 
0937   //------------------ End of pi0 in EB --------------------------//
0938 
0939   // fill EB eta histos
0940   if (isMonEBeta_) {
0941     if (rhEBeta.isValid() && (!rhEBeta->empty())) {
0942       const EcalRecHitCollection *hitCollection_p = rhEBeta.product();
0943       float etot = 0;
0944       for (itb = rhEBeta->begin(); itb != rhEBeta->end(); ++itb) {
0945         EBDetId id(itb->id());
0946         double energy = itb->energy();
0947         if (energy < seleXtalMinEnergy_)
0948           continue;
0949 
0950         EBDetId det = itb->id();
0951 
0952         detIdEBRecHits.push_back(det);
0953         EBRecHits.push_back(*itb);
0954 
0955         if (energy > clusSeedThr_)
0956           seeds.push_back(*itb);
0957 
0958         hiPhiDistrEBeta_->Fill(id.iphi());
0959         hiEtaDistrEBeta_->Fill(id.ieta());
0960         hRechitEnergyEBeta_->Fill(itb->energy());
0961 
0962         etot += itb->energy();
0963       }  // Eb rechits
0964 
0965       hNRecHitsEBeta_->Fill(rhEBeta->size());
0966       hMeanRecHitEnergyEBeta_->Fill(etot / rhEBeta->size());
0967       hEventEnergyEBeta_->Fill(etot);
0968 
0969       //      cout << " EB RH Eta collection: #, mean rh_e, event E
0970       //      "<<rhEBeta->size()<<" "<<etot/rhEBeta->size()<<" "<<etot<<endl;
0971 
0972       // Eta maker
0973 
0974       // cout<< " RH coll size: "<<rhEBeta->size()<<endl;
0975       // cout<< " Eta seeds: "<<seeds.size()<<endl;
0976 
0977       int nClus;
0978       std::vector<float> eClus;
0979       std::vector<float> etClus;
0980       std::vector<float> etaClus;
0981       std::vector<float> thetaClus;
0982       std::vector<float> phiClus;
0983       std::vector<EBDetId> max_hit;
0984 
0985       std::vector<std::vector<EcalRecHit>> RecHitsCluster;
0986       std::vector<std::vector<EcalRecHit>> RecHitsCluster5x5;
0987       std::vector<float> s4s9Clus;
0988       std::vector<float> s9s25Clus;
0989 
0990       nClus = 0;
0991 
0992       // Make own simple clusters (3x3, 5x5 or clusPhiSize_ x clusEtaSize_)
0993       std::sort(seeds.begin(), seeds.end(), ecalRecHitGreater);
0994 
0995       for (std::vector<EcalRecHit>::iterator itseed = seeds.begin(); itseed != seeds.end(); itseed++) {
0996         EBDetId seed_id = itseed->id();
0997         std::vector<EBDetId>::const_iterator usedIds;
0998 
0999         bool seedAlreadyUsed = false;
1000         for (usedIds = usedXtals.begin(); usedIds != usedXtals.end(); usedIds++) {
1001           if (*usedIds == seed_id) {
1002             seedAlreadyUsed = true;
1003             // cout<< " Seed with energy "<<itseed->energy()<<" was used
1004             // !"<<endl;
1005             break;
1006           }
1007         }
1008         if (seedAlreadyUsed)
1009           continue;
1010         std::vector<DetId> clus_v = topology_p->getWindow(seed_id, clusEtaSize_, clusPhiSize_);
1011         std::vector<std::pair<DetId, float>> clus_used;
1012         // Reject the seed if not able to build the cluster around it correctly
1013         // if(clus_v.size() < clusEtaSize_*clusPhiSize_){cout<<" Not enough
1014         // RecHits "<<endl; continue;}
1015         std::vector<EcalRecHit> RecHitsInWindow;
1016         std::vector<EcalRecHit> RecHitsInWindow5x5;
1017 
1018         double simple_energy = 0;
1019 
1020         for (std::vector<DetId>::iterator det = clus_v.begin(); det != clus_v.end(); det++) {
1021           EBDetId EBdet = *det;
1022           //      cout<<" det "<< EBdet<<" ieta "<<EBdet.ieta()<<" iphi
1023           //      "<<EBdet.iphi()<<endl;
1024           bool HitAlreadyUsed = false;
1025           for (usedIds = usedXtals.begin(); usedIds != usedXtals.end(); usedIds++) {
1026             if (*usedIds == *det) {
1027               HitAlreadyUsed = true;
1028               break;
1029             }
1030           }
1031           if (HitAlreadyUsed)
1032             continue;
1033 
1034           std::vector<EBDetId>::iterator itdet = find(detIdEBRecHits.begin(), detIdEBRecHits.end(), EBdet);
1035           if (itdet == detIdEBRecHits.end())
1036             continue;
1037 
1038           int nn = int(itdet - detIdEBRecHits.begin());
1039           usedXtals.push_back(*det);
1040           RecHitsInWindow.push_back(EBRecHits[nn]);
1041           clus_used.push_back(std::make_pair(*det, 1));
1042           simple_energy = simple_energy + EBRecHits[nn].energy();
1043         }
1044 
1045         if (simple_energy <= 0)
1046           continue;
1047 
1048         math::XYZPoint clus_pos =
1049             posCalculator_.Calculate_Location(clus_used, hitCollection_p, geometry_p, geometryES_p);
1050         // cout<< "       Simple Clustering: Total energy for this simple
1051         // cluster : "<<simple_energy<<endl; cout<< "       Simple Clustering:
1052         // eta phi : "<<clus_pos.eta()<<" "<<clus_pos.phi()<<endl; cout<< "
1053         // Simple Clustering: x y z : "<<clus_pos.x()<<" "<<clus_pos.y()<<"
1054         // "<<clus_pos.z()<<endl;
1055 
1056         float theta_s = 2. * atan(exp(-clus_pos.eta()));
1057         //  float p0x_s = simple_energy * sin(theta_s) *
1058         // cos(clus_pos.phi()); float p0y_s = simple_energy * sin(theta_s) *
1059         // sin(clus_pos.phi());
1060         //    float p0z_s = simple_energy * cos(theta_s);
1061         // float et_s = sqrt( p0x_s*p0x_s + p0y_s*p0y_s);
1062 
1063         float et_s = simple_energy * sin(theta_s);
1064         // cout << "       Simple Clustering: E,Et,px,py,pz: "<<simple_energy<<"
1065         // "<<et_s<<" "<<p0x_s<<" "<<p0y_s<<" "<<endl;
1066 
1067         // Compute S4/S9 variable
1068         // We are not sure to have 9 RecHits so need to check eta and phi:
1069         /// check s4s9
1070         float s4s9_tmp[4];
1071         for (int i = 0; i < 4; i++)
1072           s4s9_tmp[i] = 0;
1073 
1074         int seed_ieta = seed_id.ieta();
1075         int seed_iphi = seed_id.iphi();
1076 
1077         convxtalid(seed_iphi, seed_ieta);
1078 
1079         float e3x3 = 0;
1080         float e5x5 = 0;
1081 
1082         for (unsigned int j = 0; j < RecHitsInWindow.size(); j++) {
1083           EBDetId det = (EBDetId)RecHitsInWindow[j].id();
1084 
1085           int ieta = det.ieta();
1086           int iphi = det.iphi();
1087 
1088           convxtalid(iphi, ieta);
1089 
1090           float en = RecHitsInWindow[j].energy();
1091 
1092           int dx = diff_neta_s(seed_ieta, ieta);
1093           int dy = diff_nphi_s(seed_iphi, iphi);
1094 
1095           if (dx <= 0 && dy <= 0)
1096             s4s9_tmp[0] += en;
1097           if (dx >= 0 && dy <= 0)
1098             s4s9_tmp[1] += en;
1099           if (dx <= 0 && dy >= 0)
1100             s4s9_tmp[2] += en;
1101           if (dx >= 0 && dy >= 0)
1102             s4s9_tmp[3] += en;
1103 
1104           if (std::abs(dx) <= 1 && std::abs(dy) <= 1)
1105             e3x3 += en;
1106           if (std::abs(dx) <= 2 && std::abs(dy) <= 2)
1107             e5x5 += en;
1108         }
1109 
1110         if (e3x3 <= 0)
1111           continue;
1112 
1113         float s4s9_max = *std::max_element(s4s9_tmp, s4s9_tmp + 4) / e3x3;
1114 
1115         /// calculate e5x5
1116         std::vector<DetId> clus_v5x5 = topology_p->getWindow(seed_id, 5, 5);
1117         for (std::vector<DetId>::const_iterator idItr = clus_v5x5.begin(); idItr != clus_v5x5.end(); idItr++) {
1118           EBDetId det = *idItr;
1119 
1120           std::vector<EBDetId>::iterator itdet0 = find(usedXtals.begin(), usedXtals.end(), det);
1121 
1122           /// already clustered
1123           if (itdet0 != usedXtals.end())
1124             continue;
1125 
1126           // inside collections
1127           std::vector<EBDetId>::iterator itdet = find(detIdEBRecHits.begin(), detIdEBRecHits.end(), det);
1128           if (itdet == detIdEBRecHits.end())
1129             continue;
1130 
1131           int nn = int(itdet - detIdEBRecHits.begin());
1132 
1133           RecHitsInWindow5x5.push_back(EBRecHits[nn]);
1134           e5x5 += EBRecHits[nn].energy();
1135         }
1136 
1137         if (e5x5 <= 0)
1138           continue;
1139 
1140         eClus.push_back(simple_energy);
1141         etClus.push_back(et_s);
1142         etaClus.push_back(clus_pos.eta());
1143         thetaClus.push_back(theta_s);
1144         phiClus.push_back(clus_pos.phi());
1145         s4s9Clus.push_back(s4s9_max);
1146         s9s25Clus.push_back(e3x3 / e5x5);
1147         RecHitsCluster.push_back(RecHitsInWindow);
1148         RecHitsCluster5x5.push_back(RecHitsInWindow5x5);
1149 
1150         //      std::cout<<" EB Eta cluster (n,nxt,e,et eta,phi,s4s9)
1151         //"<<nClus<<" "<<int(RecHitsInWindow.size())<<" "<<eClus[nClus]<<" "<<"
1152         //"<<etClus[nClus]<<" "<<etaClus[nClus]<<" "<<phiClus[nClus]<<"
1153         //"<<s4s9Clus[nClus]<<std::endl;
1154 
1155         nClus++;
1156       }
1157 
1158       // cout<< " Eta clusters: "<<nClus<<endl;
1159 
1160       // Selection, based on Simple clustering
1161       // eta candidates
1162       int npi0_s = 0;
1163 
1164       //      if (nClus <= 1) return;
1165       for (Int_t i = 0; i < nClus; i++) {
1166         for (Int_t j = i + 1; j < nClus; j++) {
1167           //      cout<<" i "<<i<<"  etClus[i] "<<etClus[i]<<" j "<<j<<"
1168           //      etClus[j] "<<etClus[j]<<endl;
1169           if (etClus[i] > selePtGammaEta_ && etClus[j] > selePtGammaEta_ && s4s9Clus[i] > seleS4S9GammaEta_ &&
1170               s4s9Clus[j] > seleS4S9GammaEta_) {
1171             float p0x = etClus[i] * cos(phiClus[i]);
1172             float p1x = etClus[j] * cos(phiClus[j]);
1173             float p0y = etClus[i] * sin(phiClus[i]);
1174             float p1y = etClus[j] * sin(phiClus[j]);
1175             float p0z = eClus[i] * cos(thetaClus[i]);
1176             float p1z = eClus[j] * cos(thetaClus[j]);
1177 
1178             float pt_pair = sqrt((p0x + p1x) * (p0x + p1x) + (p0y + p1y) * (p0y + p1y));
1179 
1180             if (pt_pair < selePtEta_)
1181               continue;
1182 
1183             float m_inv = sqrt((eClus[i] + eClus[j]) * (eClus[i] + eClus[j]) - (p0x + p1x) * (p0x + p1x) -
1184                                (p0y + p1y) * (p0y + p1y) - (p0z + p1z) * (p0z + p1z));
1185             if ((m_inv < seleMinvMaxEta_) && (m_inv > seleMinvMinEta_)) {
1186               // New Loop on cluster to measure isolation:
1187               std::vector<int> IsoClus;
1188               IsoClus.clear();
1189               float Iso = 0;
1190               TVector3 pairVect = TVector3((p0x + p1x), (p0y + p1y), (p0z + p1z));
1191               for (Int_t k = 0; k < nClus; k++) {
1192                 if (etClus[k] < ptMinForIsolationEta_)
1193                   continue;
1194 
1195                 if (k == i || k == j)
1196                   continue;
1197                 TVector3 ClusVect =
1198                     TVector3(etClus[k] * cos(phiClus[k]), etClus[k] * sin(phiClus[k]), eClus[k] * cos(thetaClus[k]));
1199 
1200                 float dretacl = fabs(etaClus[k] - pairVect.Eta());
1201                 float drcl = ClusVect.DeltaR(pairVect);
1202                 //      cout<< "   Iso: k, E, drclpi0, detaclpi0, dphiclpi0
1203                 //      "<<k<<" "<<eClus[k]<<" "<<drclpi0<<"
1204                 //      "<<dretaclpi0<<endl;
1205                 if ((drcl < seleEtaBeltDR_) && (dretacl < seleEtaBeltDeta_)) {
1206                   //              cout<< "   ... good iso cluster #: "<<k<<"
1207                   //              etClus[k] "<<etClus[k] <<endl;
1208                   Iso = Iso + etClus[k];
1209                   IsoClus.push_back(k);
1210                 }
1211               }
1212 
1213               //      cout<<"  Iso/pt_pi0 "<<Iso/pt_pi0<<endl;
1214               if (Iso / pt_pair < seleEtaIso_) {
1215                 // for(unsigned int
1216                 // Rec=0;Rec<RecHitsCluster[i].size();Rec++)pi0EBRecHitCollection->push_back(RecHitsCluster[i][Rec]);
1217                 // for(unsigned int
1218                 // Rec2=0;Rec2<RecHitsCluster[j].size();Rec2++)pi0EBRecHitCollection->push_back(RecHitsCluster[j][Rec2]);
1219 
1220                 hMinvEtaEB_->Fill(m_inv);
1221                 hPt1EtaEB_->Fill(etClus[i]);
1222                 hPt2EtaEB_->Fill(etClus[j]);
1223                 hPtEtaEB_->Fill(pt_pair);
1224                 hIsoEtaEB_->Fill(Iso / pt_pair);
1225                 hS4S91EtaEB_->Fill(s4s9Clus[i]);
1226                 hS4S92EtaEB_->Fill(s4s9Clus[j]);
1227 
1228                 //      cout <<"  EB Simple Clustering: Eta Candidate
1229                 // pt, eta, phi, Iso, m_inv, i, j :   "<<pt_pair<<"
1230                 //"<<pairVect.Eta()<<" "<<pairVect.Phi()<<" "<<Iso<<"
1231                 //"<<m_inv<<" "<<i<<" "<<j<<" "<<endl;
1232 
1233                 npi0_s++;
1234               }
1235             }
1236           }
1237         }  // End of the "j" loop over Simple Clusters
1238       }    // End of the "i" loop over Simple Clusters
1239 
1240       //      cout<<"  (Simple Clustering) EB Eta candidates #: "<<npi0_s<<endl;
1241 
1242     }  // rhEBeta.valid() ends
1243 
1244   }  //  isMonEBeta ends
1245 
1246   //------------------ End of Eta in EB --------------------------//
1247 
1248   //----------------- End of the EB --------------------------//
1249 
1250   //----------------- Start of the EE --------------------//
1251 
1252   // fill pi0 EE histos
1253   if (isMonEEpi0_) {
1254     if (rhEEpi0.isValid() && (!rhEEpi0->empty())) {
1255       const EcalRecHitCollection *hitCollection_ee = rhEEpi0.product();
1256       float etot = 0;
1257 
1258       detIdEERecHits.clear();  //// EEDetId
1259       EERecHits.clear();       /// EcalRecHit
1260 
1261       std::vector<EcalRecHit> seedsEndCap;
1262       seedsEndCap.clear();
1263 
1264       std::vector<EEDetId> usedXtalsEndCap;
1265       usedXtalsEndCap.clear();
1266 
1267       ////make seeds.
1268       EERecHitCollection::const_iterator ite;
1269       for (ite = rhEEpi0->begin(); ite != rhEEpi0->end(); ite++) {
1270         double energy = ite->energy();
1271         if (energy < seleXtalMinEnergyEndCap_)
1272           continue;
1273 
1274         EEDetId det = ite->id();
1275         EEDetId id(ite->id());
1276 
1277         detIdEERecHits.push_back(det);
1278         EERecHits.push_back(*ite);
1279 
1280         hiXDistrEEpi0_->Fill(id.ix());
1281         hiYDistrEEpi0_->Fill(id.iy());
1282         hRechitEnergyEEpi0_->Fill(ite->energy());
1283 
1284         if (energy > clusSeedThrEndCap_)
1285           seedsEndCap.push_back(*ite);
1286 
1287         etot += ite->energy();
1288       }  // EE rechits
1289 
1290       hNRecHitsEEpi0_->Fill(rhEEpi0->size());
1291       hMeanRecHitEnergyEEpi0_->Fill(etot / rhEEpi0->size());
1292       hEventEnergyEEpi0_->Fill(etot);
1293 
1294       //      cout << " EE RH Pi0 collection: #, mean rh_e, event E
1295       //      "<<rhEEpi0->size()<<" "<<etot/rhEEpi0->size()<<" "<<etot<<endl;
1296 
1297       int nClusEndCap;
1298       std::vector<float> eClusEndCap;
1299       std::vector<float> etClusEndCap;
1300       std::vector<float> etaClusEndCap;
1301       std::vector<float> thetaClusEndCap;
1302       std::vector<float> phiClusEndCap;
1303       std::vector<std::vector<EcalRecHit>> RecHitsClusterEndCap;
1304       std::vector<std::vector<EcalRecHit>> RecHitsCluster5x5EndCap;
1305       std::vector<float> s4s9ClusEndCap;
1306       std::vector<float> s9s25ClusEndCap;
1307 
1308       nClusEndCap = 0;
1309 
1310       // Make own simple clusters (3x3, 5x5 or clusPhiSize_ x clusEtaSize_)
1311       std::sort(seedsEndCap.begin(), seedsEndCap.end(), ecalRecHitGreater);
1312 
1313       for (std::vector<EcalRecHit>::iterator itseed = seedsEndCap.begin(); itseed != seedsEndCap.end(); itseed++) {
1314         EEDetId seed_id = itseed->id();
1315         std::vector<EEDetId>::const_iterator usedIds;
1316 
1317         bool seedAlreadyUsed = false;
1318         for (usedIds = usedXtalsEndCap.begin(); usedIds != usedXtalsEndCap.end(); usedIds++) {
1319           if (*usedIds == seed_id) {
1320             seedAlreadyUsed = true;
1321             break;
1322           }
1323         }
1324 
1325         if (seedAlreadyUsed)
1326           continue;
1327         std::vector<DetId> clus_v = topology_ee->getWindow(seed_id, clusEtaSize_, clusPhiSize_);
1328         std::vector<std::pair<DetId, float>> clus_used;
1329 
1330         std::vector<EcalRecHit> RecHitsInWindow;
1331         std::vector<EcalRecHit> RecHitsInWindow5x5;
1332 
1333         float simple_energy = 0;
1334 
1335         for (std::vector<DetId>::iterator det = clus_v.begin(); det != clus_v.end(); det++) {
1336           EEDetId EEdet = *det;
1337 
1338           bool HitAlreadyUsed = false;
1339           for (usedIds = usedXtalsEndCap.begin(); usedIds != usedXtalsEndCap.end(); usedIds++) {
1340             if (*usedIds == *det) {
1341               HitAlreadyUsed = true;
1342               break;
1343             }
1344           }
1345 
1346           if (HitAlreadyUsed)
1347             continue;
1348 
1349           std::vector<EEDetId>::iterator itdet = find(detIdEERecHits.begin(), detIdEERecHits.end(), EEdet);
1350           if (itdet == detIdEERecHits.end())
1351             continue;
1352 
1353           int nn = int(itdet - detIdEERecHits.begin());
1354           usedXtalsEndCap.push_back(*det);
1355           RecHitsInWindow.push_back(EERecHits[nn]);
1356           clus_used.push_back(std::make_pair(*det, 1));
1357           simple_energy = simple_energy + EERecHits[nn].energy();
1358         }
1359 
1360         if (simple_energy <= 0)
1361           continue;
1362 
1363         math::XYZPoint clus_pos =
1364             posCalculator_.Calculate_Location(clus_used, hitCollection_ee, geometryEE_p, geometryES_p);
1365 
1366         float theta_s = 2. * atan(exp(-clus_pos.eta()));
1367         float et_s = simple_energy * sin(theta_s);
1368         //  float p0x_s = simple_energy * sin(theta_s) *
1369         // cos(clus_pos.phi()); float p0y_s = simple_energy * sin(theta_s) *
1370         // sin(clus_pos.phi()); float et_s = sqrt( p0x_s*p0x_s + p0y_s*p0y_s);
1371 
1372         // Compute S4/S9 variable
1373         // We are not sure to have 9 RecHits so need to check eta and phi:
1374         float s4s9_tmp[4];
1375         for (int i = 0; i < 4; i++)
1376           s4s9_tmp[i] = 0;
1377 
1378         int ixSeed = seed_id.ix();
1379         int iySeed = seed_id.iy();
1380         float e3x3 = 0;
1381         float e5x5 = 0;
1382 
1383         for (unsigned int j = 0; j < RecHitsInWindow.size(); j++) {
1384           EEDetId det_this = (EEDetId)RecHitsInWindow[j].id();
1385           int dx = ixSeed - det_this.ix();
1386           int dy = iySeed - det_this.iy();
1387 
1388           float en = RecHitsInWindow[j].energy();
1389 
1390           if (dx <= 0 && dy <= 0)
1391             s4s9_tmp[0] += en;
1392           if (dx >= 0 && dy <= 0)
1393             s4s9_tmp[1] += en;
1394           if (dx <= 0 && dy >= 0)
1395             s4s9_tmp[2] += en;
1396           if (dx >= 0 && dy >= 0)
1397             s4s9_tmp[3] += en;
1398 
1399           if (std::abs(dx) <= 1 && std::abs(dy) <= 1)
1400             e3x3 += en;
1401           if (std::abs(dx) <= 2 && std::abs(dy) <= 2)
1402             e5x5 += en;
1403         }
1404 
1405         if (e3x3 <= 0)
1406           continue;
1407 
1408         eClusEndCap.push_back(simple_energy);
1409         etClusEndCap.push_back(et_s);
1410         etaClusEndCap.push_back(clus_pos.eta());
1411         thetaClusEndCap.push_back(theta_s);
1412         phiClusEndCap.push_back(clus_pos.phi());
1413         s4s9ClusEndCap.push_back(*std::max_element(s4s9_tmp, s4s9_tmp + 4) / e3x3);
1414         s9s25ClusEndCap.push_back(e3x3 / e5x5);
1415         RecHitsClusterEndCap.push_back(RecHitsInWindow);
1416         RecHitsCluster5x5EndCap.push_back(RecHitsInWindow5x5);
1417 
1418         //  std::cout<<" EE pi0 cluster (n,nxt,e,et eta,phi,s4s9)
1419         //"<<nClusEndCap<<" "<<int(RecHitsInWindow.size())<<"
1420         //"<<eClusEndCap[nClusEndCap]<<" "<<" "<<etClusEndCap[nClusEndCap]<<"
1421         //"<<etaClusEndCap[nClusEndCap]<<" "<<phiClusEndCap[nClusEndCap]<<"
1422         //"<<s4s9ClusEndCap[nClusEndCap]<<std::endl;
1423 
1424         nClusEndCap++;
1425       }
1426 
1427       // Selection, based on Simple clustering
1428       // pi0 candidates
1429       int npi0_se = 0;
1430 
1431       for (Int_t i = 0; i < nClusEndCap; i++) {
1432         for (Int_t j = i + 1; j < nClusEndCap; j++) {
1433           if (etClusEndCap[i] > selePtGammaEndCap_ && etClusEndCap[j] > selePtGammaEndCap_ &&
1434               s4s9ClusEndCap[i] > seleS4S9GammaEndCap_ && s4s9ClusEndCap[j] > seleS4S9GammaEndCap_) {
1435             float p0x = etClusEndCap[i] * cos(phiClusEndCap[i]);
1436             float p1x = etClusEndCap[j] * cos(phiClusEndCap[j]);
1437             float p0y = etClusEndCap[i] * sin(phiClusEndCap[i]);
1438             float p1y = etClusEndCap[j] * sin(phiClusEndCap[j]);
1439             float p0z = eClusEndCap[i] * cos(thetaClusEndCap[i]);
1440             float p1z = eClusEndCap[j] * cos(thetaClusEndCap[j]);
1441 
1442             float pt_pair = sqrt((p0x + p1x) * (p0x + p1x) + (p0y + p1y) * (p0y + p1y));
1443             if (pt_pair < selePtPi0EndCap_)
1444               continue;
1445             float m_inv = sqrt((eClusEndCap[i] + eClusEndCap[j]) * (eClusEndCap[i] + eClusEndCap[j]) -
1446                                (p0x + p1x) * (p0x + p1x) - (p0y + p1y) * (p0y + p1y) - (p0z + p1z) * (p0z + p1z));
1447 
1448             if ((m_inv < seleMinvMaxPi0EndCap_) && (m_inv > seleMinvMinPi0EndCap_)) {
1449               // New Loop on cluster to measure isolation:
1450               std::vector<int> IsoClus;
1451               IsoClus.clear();
1452               float Iso = 0;
1453               TVector3 pairVect = TVector3((p0x + p1x), (p0y + p1y), (p0z + p1z));
1454               for (Int_t k = 0; k < nClusEndCap; k++) {
1455                 if (etClusEndCap[k] < ptMinForIsolationEndCap_)
1456                   continue;
1457 
1458                 if (k == i || k == j)
1459                   continue;
1460 
1461                 TVector3 clusVect = TVector3(etClusEndCap[k] * cos(phiClusEndCap[k]),
1462                                              etClusEndCap[k] * sin(phiClusEndCap[k]),
1463                                              eClusEndCap[k] * cos(thetaClusEndCap[k]));
1464                 float dretacl = fabs(etaClusEndCap[k] - pairVect.Eta());
1465                 float drcl = clusVect.DeltaR(pairVect);
1466 
1467                 if (drcl < selePi0BeltDREndCap_ && dretacl < selePi0BeltDetaEndCap_) {
1468                   Iso = Iso + etClusEndCap[k];
1469                   IsoClus.push_back(k);
1470                 }
1471               }
1472 
1473               if (Iso / pt_pair < selePi0IsoEndCap_) {
1474                 // cout <<"  EE Simple Clustering: pi0 Candidate pt, eta, phi,
1475                 // Iso, m_inv, i, j :   "<<pt_pair<<" "<<pairVect.Eta()<<"
1476                 // "<<pairVect.Phi()<<" "<<Iso<<" "<<m_inv<<" "<<i<<" "<<j<<"
1477                 // "<<endl;
1478 
1479                 hMinvPi0EE_->Fill(m_inv);
1480                 hPt1Pi0EE_->Fill(etClusEndCap[i]);
1481                 hPt2Pi0EE_->Fill(etClusEndCap[j]);
1482                 hPtPi0EE_->Fill(pt_pair);
1483                 hIsoPi0EE_->Fill(Iso / pt_pair);
1484                 hS4S91Pi0EE_->Fill(s4s9ClusEndCap[i]);
1485                 hS4S92Pi0EE_->Fill(s4s9ClusEndCap[j]);
1486 
1487                 npi0_se++;
1488               }
1489             }
1490           }
1491         }  // End of the "j" loop over Simple Clusters
1492       }    // End of the "i" loop over Simple Clusters
1493 
1494       //      cout<<"  (Simple Clustering) EE Pi0 candidates #:
1495       //      "<<npi0_se<<endl;
1496 
1497     }  // rhEEpi0
1498   }    // isMonEEpi0
1499 
1500   //================End of Pi0 endcap=======================//
1501 
1502   //================ Eta in EE===============================//
1503 
1504   // fill pi0 EE histos
1505   if (isMonEEeta_) {
1506     if (rhEEeta.isValid() && (!rhEEeta->empty())) {
1507       const EcalRecHitCollection *hitCollection_ee = rhEEeta.product();
1508       float etot = 0;
1509 
1510       detIdEERecHits.clear();  //// EEDetId
1511       EERecHits.clear();       /// EcalRecHit
1512 
1513       std::vector<EcalRecHit> seedsEndCap;
1514       seedsEndCap.clear();
1515 
1516       std::vector<EEDetId> usedXtalsEndCap;
1517       usedXtalsEndCap.clear();
1518 
1519       ////make seeds.
1520       EERecHitCollection::const_iterator ite;
1521       for (ite = rhEEeta->begin(); ite != rhEEeta->end(); ite++) {
1522         double energy = ite->energy();
1523         if (energy < seleXtalMinEnergyEndCap_)
1524           continue;
1525 
1526         EEDetId det = ite->id();
1527         EEDetId id(ite->id());
1528 
1529         detIdEERecHits.push_back(det);
1530         EERecHits.push_back(*ite);
1531 
1532         hiXDistrEEeta_->Fill(id.ix());
1533         hiYDistrEEeta_->Fill(id.iy());
1534         hRechitEnergyEEeta_->Fill(ite->energy());
1535 
1536         if (energy > clusSeedThrEndCap_)
1537           seedsEndCap.push_back(*ite);
1538 
1539         etot += ite->energy();
1540       }  // EE rechits
1541 
1542       hNRecHitsEEeta_->Fill(rhEEeta->size());
1543       hMeanRecHitEnergyEEeta_->Fill(etot / rhEEeta->size());
1544       hEventEnergyEEeta_->Fill(etot);
1545 
1546       //      cout << " EE RH Eta collection: #, mean rh_e, event E
1547       //      "<<rhEEeta->size()<<" "<<etot/rhEEeta->size()<<" "<<etot<<endl;
1548 
1549       int nClusEndCap;
1550       std::vector<float> eClusEndCap;
1551       std::vector<float> etClusEndCap;
1552       std::vector<float> etaClusEndCap;
1553       std::vector<float> thetaClusEndCap;
1554       std::vector<float> phiClusEndCap;
1555       std::vector<std::vector<EcalRecHit>> RecHitsClusterEndCap;
1556       std::vector<std::vector<EcalRecHit>> RecHitsCluster5x5EndCap;
1557       std::vector<float> s4s9ClusEndCap;
1558       std::vector<float> s9s25ClusEndCap;
1559 
1560       nClusEndCap = 0;
1561 
1562       // Make own simple clusters (3x3, 5x5 or clusPhiSize_ x clusEtaSize_)
1563       std::sort(seedsEndCap.begin(), seedsEndCap.end(), ecalRecHitGreater);
1564 
1565       for (std::vector<EcalRecHit>::iterator itseed = seedsEndCap.begin(); itseed != seedsEndCap.end(); itseed++) {
1566         EEDetId seed_id = itseed->id();
1567         std::vector<EEDetId>::const_iterator usedIds;
1568 
1569         bool seedAlreadyUsed = false;
1570         for (usedIds = usedXtalsEndCap.begin(); usedIds != usedXtalsEndCap.end(); usedIds++) {
1571           if (*usedIds == seed_id) {
1572             seedAlreadyUsed = true;
1573             break;
1574           }
1575         }
1576 
1577         if (seedAlreadyUsed)
1578           continue;
1579         std::vector<DetId> clus_v = topology_ee->getWindow(seed_id, clusEtaSize_, clusPhiSize_);
1580         std::vector<std::pair<DetId, float>> clus_used;
1581 
1582         std::vector<EcalRecHit> RecHitsInWindow;
1583         std::vector<EcalRecHit> RecHitsInWindow5x5;
1584 
1585         float simple_energy = 0;
1586 
1587         for (std::vector<DetId>::iterator det = clus_v.begin(); det != clus_v.end(); det++) {
1588           EEDetId EEdet = *det;
1589 
1590           bool HitAlreadyUsed = false;
1591           for (usedIds = usedXtalsEndCap.begin(); usedIds != usedXtalsEndCap.end(); usedIds++) {
1592             if (*usedIds == *det) {
1593               HitAlreadyUsed = true;
1594               break;
1595             }
1596           }
1597 
1598           if (HitAlreadyUsed)
1599             continue;
1600 
1601           std::vector<EEDetId>::iterator itdet = find(detIdEERecHits.begin(), detIdEERecHits.end(), EEdet);
1602           if (itdet == detIdEERecHits.end())
1603             continue;
1604 
1605           int nn = int(itdet - detIdEERecHits.begin());
1606           usedXtalsEndCap.push_back(*det);
1607           RecHitsInWindow.push_back(EERecHits[nn]);
1608           clus_used.push_back(std::make_pair(*det, 1));
1609           simple_energy = simple_energy + EERecHits[nn].energy();
1610         }
1611 
1612         if (simple_energy <= 0)
1613           continue;
1614 
1615         math::XYZPoint clus_pos =
1616             posCalculator_.Calculate_Location(clus_used, hitCollection_ee, geometryEE_p, geometryES_p);
1617 
1618         float theta_s = 2. * atan(exp(-clus_pos.eta()));
1619         float et_s = simple_energy * sin(theta_s);
1620         //  float p0x_s = simple_energy * sin(theta_s) *
1621         // cos(clus_pos.phi()); float p0y_s = simple_energy * sin(theta_s) *
1622         // sin(clus_pos.phi()); float et_s = sqrt( p0x_s*p0x_s + p0y_s*p0y_s);
1623 
1624         // Compute S4/S9 variable
1625         // We are not sure to have 9 RecHits so need to check eta and phi:
1626         float s4s9_tmp[4];
1627         for (int i = 0; i < 4; i++)
1628           s4s9_tmp[i] = 0;
1629 
1630         int ixSeed = seed_id.ix();
1631         int iySeed = seed_id.iy();
1632         float e3x3 = 0;
1633         float e5x5 = 0;
1634 
1635         for (unsigned int j = 0; j < RecHitsInWindow.size(); j++) {
1636           EEDetId det_this = (EEDetId)RecHitsInWindow[j].id();
1637           int dx = ixSeed - det_this.ix();
1638           int dy = iySeed - det_this.iy();
1639 
1640           float en = RecHitsInWindow[j].energy();
1641 
1642           if (dx <= 0 && dy <= 0)
1643             s4s9_tmp[0] += en;
1644           if (dx >= 0 && dy <= 0)
1645             s4s9_tmp[1] += en;
1646           if (dx <= 0 && dy >= 0)
1647             s4s9_tmp[2] += en;
1648           if (dx >= 0 && dy >= 0)
1649             s4s9_tmp[3] += en;
1650 
1651           if (std::abs(dx) <= 1 && std::abs(dy) <= 1)
1652             e3x3 += en;
1653           if (std::abs(dx) <= 2 && std::abs(dy) <= 2)
1654             e5x5 += en;
1655         }
1656 
1657         if (e3x3 <= 0)
1658           continue;
1659 
1660         eClusEndCap.push_back(simple_energy);
1661         etClusEndCap.push_back(et_s);
1662         etaClusEndCap.push_back(clus_pos.eta());
1663         thetaClusEndCap.push_back(theta_s);
1664         phiClusEndCap.push_back(clus_pos.phi());
1665         s4s9ClusEndCap.push_back(*std::max_element(s4s9_tmp, s4s9_tmp + 4) / e3x3);
1666         s9s25ClusEndCap.push_back(e3x3 / e5x5);
1667         RecHitsClusterEndCap.push_back(RecHitsInWindow);
1668         RecHitsCluster5x5EndCap.push_back(RecHitsInWindow5x5);
1669 
1670         //  std::cout<<" EE Eta cluster (n,nxt,e,et eta,phi,s4s9)
1671         //"<<nClusEndCap<<" "<<int(RecHitsInWindow.size())<<"
1672         //"<<eClusEndCap[nClusEndCap]<<" "<<" "<<etClusEndCap[nClusEndCap]<<"
1673         //"<<etaClusEndCap[nClusEndCap]<<" "<<phiClusEndCap[nClusEndCap]<<"
1674         //"<<s4s9ClusEndCap[nClusEndCap]<<std::endl;
1675 
1676         nClusEndCap++;
1677       }
1678 
1679       // Selection, based on Simple clustering
1680       // pi0 candidates
1681       int npi0_se = 0;
1682 
1683       for (Int_t i = 0; i < nClusEndCap; i++) {
1684         for (Int_t j = i + 1; j < nClusEndCap; j++) {
1685           if (etClusEndCap[i] > selePtGammaEtaEndCap_ && etClusEndCap[j] > selePtGammaEtaEndCap_ &&
1686               s4s9ClusEndCap[i] > seleS4S9GammaEtaEndCap_ && s4s9ClusEndCap[j] > seleS4S9GammaEtaEndCap_) {
1687             float p0x = etClusEndCap[i] * cos(phiClusEndCap[i]);
1688             float p1x = etClusEndCap[j] * cos(phiClusEndCap[j]);
1689             float p0y = etClusEndCap[i] * sin(phiClusEndCap[i]);
1690             float p1y = etClusEndCap[j] * sin(phiClusEndCap[j]);
1691             float p0z = eClusEndCap[i] * cos(thetaClusEndCap[i]);
1692             float p1z = eClusEndCap[j] * cos(thetaClusEndCap[j]);
1693 
1694             float pt_pair = sqrt((p0x + p1x) * (p0x + p1x) + (p0y + p1y) * (p0y + p1y));
1695             if (pt_pair < selePtEtaEndCap_)
1696               continue;
1697             float m_inv = sqrt((eClusEndCap[i] + eClusEndCap[j]) * (eClusEndCap[i] + eClusEndCap[j]) -
1698                                (p0x + p1x) * (p0x + p1x) - (p0y + p1y) * (p0y + p1y) - (p0z + p1z) * (p0z + p1z));
1699 
1700             if ((m_inv < seleMinvMaxEtaEndCap_) && (m_inv > seleMinvMinEtaEndCap_)) {
1701               // New Loop on cluster to measure isolation:
1702               std::vector<int> IsoClus;
1703               IsoClus.clear();
1704               float Iso = 0;
1705               TVector3 pairVect = TVector3((p0x + p1x), (p0y + p1y), (p0z + p1z));
1706               for (Int_t k = 0; k < nClusEndCap; k++) {
1707                 if (etClusEndCap[k] < ptMinForIsolationEtaEndCap_)
1708                   continue;
1709 
1710                 if (k == i || k == j)
1711                   continue;
1712 
1713                 TVector3 clusVect = TVector3(etClusEndCap[k] * cos(phiClusEndCap[k]),
1714                                              etClusEndCap[k] * sin(phiClusEndCap[k]),
1715                                              eClusEndCap[k] * cos(thetaClusEndCap[k]));
1716                 float dretacl = fabs(etaClusEndCap[k] - pairVect.Eta());
1717                 float drcl = clusVect.DeltaR(pairVect);
1718 
1719                 if (drcl < seleEtaBeltDREndCap_ && dretacl < seleEtaBeltDetaEndCap_) {
1720                   Iso = Iso + etClusEndCap[k];
1721                   IsoClus.push_back(k);
1722                 }
1723               }
1724 
1725               if (Iso / pt_pair < seleEtaIsoEndCap_) {
1726                 //  cout <<"  EE Simple Clustering: Eta Candidate pt, eta,
1727                 // phi, Iso, m_inv, i, j :   "<<pt_pair<<" "<<pairVect.Eta()<<"
1728                 //"<<pairVect.Phi()<<" "<<Iso<<" "<<m_inv<<" "<<i<<" "<<j<<"
1729                 //"<<endl;
1730 
1731                 hMinvEtaEE_->Fill(m_inv);
1732                 hPt1EtaEE_->Fill(etClusEndCap[i]);
1733                 hPt2EtaEE_->Fill(etClusEndCap[j]);
1734                 hPtEtaEE_->Fill(pt_pair);
1735                 hIsoEtaEE_->Fill(Iso / pt_pair);
1736                 hS4S91EtaEE_->Fill(s4s9ClusEndCap[i]);
1737                 hS4S92EtaEE_->Fill(s4s9ClusEndCap[j]);
1738 
1739                 npi0_se++;
1740               }
1741             }
1742           }
1743         }  // End of the "j" loop over Simple Clusters
1744       }    // End of the "i" loop over Simple Clusters
1745 
1746       //      cout<<"  (Simple Clustering) EE Eta candidates #:
1747       //      "<<npi0_se<<endl;
1748 
1749     }  // rhEEeta
1750   }    // isMonEEeta
1751 
1752   //================End of Pi0 endcap=======================//
1753 
1754   ////==============End of endcap ===============///
1755 }
1756 
1757 void DQMSourcePi0::convxtalid(Int_t &nphi, Int_t &neta) {
1758   // Barrel only
1759   // Output nphi 0...359; neta 0...84; nside=+1 (for eta>0), or 0 (for eta<0).
1760   // neta will be [-85,-1] , or [0,84], the minus sign indicates the z<0 side.
1761 
1762   if (neta > 0)
1763     neta -= 1;
1764   if (nphi > 359)
1765     nphi = nphi - 360;
1766 
1767 }  // end of convxtalid
1768 
1769 int DQMSourcePi0::diff_neta_s(Int_t neta1, Int_t neta2) {
1770   Int_t mdiff;
1771   mdiff = (neta1 - neta2);
1772   return mdiff;
1773 }
1774 
1775 // Calculate the distance in xtals taking into account the periodicity of the
1776 // Barrel
1777 int DQMSourcePi0::diff_nphi_s(Int_t nphi1, Int_t nphi2) {
1778   Int_t mdiff;
1779   if (std::abs(nphi1 - nphi2) < (360 - std::abs(nphi1 - nphi2))) {
1780     mdiff = nphi1 - nphi2;
1781   } else {
1782     mdiff = 360 - std::abs(nphi1 - nphi2);
1783     if (nphi1 > nphi2)
1784       mdiff = -mdiff;
1785   }
1786   return mdiff;
1787 }
1788 
1789 DEFINE_FWK_MODULE(DQMSourcePi0);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!