Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-11 04:32:44 

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 
0854       //      if (nClus <= 1) return;
0855       for (Int_t i = 0; i < nClus; i++) {
0856         for (Int_t j = i + 1; j < nClus; j++) {
0857           //      cout<<" i "<<i<<"  etClus[i] "<<etClus[i]<<" j "<<j<<"
0858           //      etClus[j] "<<etClus[j]<<endl;
0859           if (etClus[i] > selePtGamma_ && etClus[j] > selePtGamma_ && s4s9Clus[i] > seleS4S9Gamma_ &&
0860               s4s9Clus[j] > seleS4S9Gamma_) {
0861             float p0x = etClus[i] * cos(phiClus[i]);
0862             float p1x = etClus[j] * cos(phiClus[j]);
0863             float p0y = etClus[i] * sin(phiClus[i]);
0864             float p1y = etClus[j] * sin(phiClus[j]);
0865             float p0z = eClus[i] * cos(thetaClus[i]);
0866             float p1z = eClus[j] * cos(thetaClus[j]);
0867 
0868             float pt_pair = sqrt((p0x + p1x) * (p0x + p1x) + (p0y + p1y) * (p0y + p1y));
0869 
0870             if (pt_pair < selePtPi0_)
0871               continue;
0872 
0873             float m_inv = sqrt((eClus[i] + eClus[j]) * (eClus[i] + eClus[j]) - (p0x + p1x) * (p0x + p1x) -
0874                                (p0y + p1y) * (p0y + p1y) - (p0z + p1z) * (p0z + p1z));
0875             if ((m_inv < seleMinvMaxPi0_) && (m_inv > seleMinvMinPi0_)) {
0876               // New Loop on cluster to measure isolation:
0877               std::vector<int> IsoClus;
0878               IsoClus.clear();
0879               float Iso = 0;
0880               TVector3 pairVect = TVector3((p0x + p1x), (p0y + p1y), (p0z + p1z));
0881               for (Int_t k = 0; k < nClus; k++) {
0882                 if (etClus[k] < ptMinForIsolation_)
0883                   continue;
0884 
0885                 if (k == i || k == j)
0886                   continue;
0887                 TVector3 ClusVect =
0888                     TVector3(etClus[k] * cos(phiClus[k]), etClus[k] * sin(phiClus[k]), eClus[k] * cos(thetaClus[k]));
0889 
0890                 float dretacl = fabs(etaClus[k] - pairVect.Eta());
0891                 float drcl = ClusVect.DeltaR(pairVect);
0892                 //      cout<< "   Iso: k, E, drclpi0, detaclpi0, dphiclpi0
0893                 //      "<<k<<" "<<eClus[k]<<" "<<drclpi0<<"
0894                 //      "<<dretaclpi0<<endl;
0895                 if ((drcl < selePi0BeltDR_) && (dretacl < selePi0BeltDeta_)) {
0896                   //              cout<< "   ... good iso cluster #: "<<k<<"
0897                   //              etClus[k] "<<etClus[k] <<endl;
0898                   Iso = Iso + etClus[k];
0899                   IsoClus.push_back(k);
0900                 }
0901               }
0902 
0903               //      cout<<"  Iso/pt_pi0 "<<Iso/pt_pi0<<endl;
0904               if (Iso / pt_pair < selePi0Iso_) {
0905                 // for(unsigned int
0906                 // Rec=0;Rec<RecHitsCluster[i].size();Rec++)pi0EBRecHitCollection->push_back(RecHitsCluster[i][Rec]);
0907                 // for(unsigned int
0908                 // Rec2=0;Rec2<RecHitsCluster[j].size();Rec2++)pi0EBRecHitCollection->push_back(RecHitsCluster[j][Rec2]);
0909 
0910                 hMinvPi0EB_->Fill(m_inv);
0911                 hPt1Pi0EB_->Fill(etClus[i]);
0912                 hPt2Pi0EB_->Fill(etClus[j]);
0913                 hPtPi0EB_->Fill(pt_pair);
0914                 hIsoPi0EB_->Fill(Iso / pt_pair);
0915                 hS4S91Pi0EB_->Fill(s4s9Clus[i]);
0916                 hS4S92Pi0EB_->Fill(s4s9Clus[j]);
0917 
0918                 //      cout <<"  EB Simple Clustering: pi0 Candidate
0919                 // pt, eta, phi, Iso, m_inv, i, j :   "<<pt_pair<<"
0920                 //"<<pairVect.Eta()<<" "<<pairVect.Phi()<<" "<<Iso<<"
0921                 //"<<m_inv<<" "<<i<<" "<<j<<" "<<endl;
0922               }
0923             }
0924           }
0925         }  // End of the "j" loop over Simple Clusters
0926       }  // End of the "i" loop over Simple Clusters
0927 
0928     }  // rhEBpi0.valid() ends
0929 
0930   }  //  isMonEBpi0 ends
0931 
0932   //------------------ End of pi0 in EB --------------------------//
0933 
0934   // fill EB eta histos
0935   if (isMonEBeta_) {
0936     if (rhEBeta.isValid() && (!rhEBeta->empty())) {
0937       const EcalRecHitCollection *hitCollection_p = rhEBeta.product();
0938       float etot = 0;
0939       for (itb = rhEBeta->begin(); itb != rhEBeta->end(); ++itb) {
0940         EBDetId id(itb->id());
0941         double energy = itb->energy();
0942         if (energy < seleXtalMinEnergy_)
0943           continue;
0944 
0945         EBDetId det = itb->id();
0946 
0947         detIdEBRecHits.push_back(det);
0948         EBRecHits.push_back(*itb);
0949 
0950         if (energy > clusSeedThr_)
0951           seeds.push_back(*itb);
0952 
0953         hiPhiDistrEBeta_->Fill(id.iphi());
0954         hiEtaDistrEBeta_->Fill(id.ieta());
0955         hRechitEnergyEBeta_->Fill(itb->energy());
0956 
0957         etot += itb->energy();
0958       }  // Eb rechits
0959 
0960       hNRecHitsEBeta_->Fill(rhEBeta->size());
0961       hMeanRecHitEnergyEBeta_->Fill(etot / rhEBeta->size());
0962       hEventEnergyEBeta_->Fill(etot);
0963 
0964       //      cout << " EB RH Eta collection: #, mean rh_e, event E
0965       //      "<<rhEBeta->size()<<" "<<etot/rhEBeta->size()<<" "<<etot<<endl;
0966 
0967       // Eta maker
0968 
0969       // cout<< " RH coll size: "<<rhEBeta->size()<<endl;
0970       // cout<< " Eta seeds: "<<seeds.size()<<endl;
0971 
0972       int nClus;
0973       std::vector<float> eClus;
0974       std::vector<float> etClus;
0975       std::vector<float> etaClus;
0976       std::vector<float> thetaClus;
0977       std::vector<float> phiClus;
0978       std::vector<EBDetId> max_hit;
0979 
0980       std::vector<std::vector<EcalRecHit>> RecHitsCluster;
0981       std::vector<std::vector<EcalRecHit>> RecHitsCluster5x5;
0982       std::vector<float> s4s9Clus;
0983       std::vector<float> s9s25Clus;
0984 
0985       nClus = 0;
0986 
0987       // Make own simple clusters (3x3, 5x5 or clusPhiSize_ x clusEtaSize_)
0988       std::sort(seeds.begin(), seeds.end(), ecalRecHitGreater);
0989 
0990       for (std::vector<EcalRecHit>::iterator itseed = seeds.begin(); itseed != seeds.end(); itseed++) {
0991         EBDetId seed_id = itseed->id();
0992         std::vector<EBDetId>::const_iterator usedIds;
0993 
0994         bool seedAlreadyUsed = false;
0995         for (usedIds = usedXtals.begin(); usedIds != usedXtals.end(); usedIds++) {
0996           if (*usedIds == seed_id) {
0997             seedAlreadyUsed = true;
0998             // cout<< " Seed with energy "<<itseed->energy()<<" was used
0999             // !"<<endl;
1000             break;
1001           }
1002         }
1003         if (seedAlreadyUsed)
1004           continue;
1005         std::vector<DetId> clus_v = topology_p->getWindow(seed_id, clusEtaSize_, clusPhiSize_);
1006         std::vector<std::pair<DetId, float>> clus_used;
1007         // Reject the seed if not able to build the cluster around it correctly
1008         // if(clus_v.size() < clusEtaSize_*clusPhiSize_){cout<<" Not enough
1009         // RecHits "<<endl; continue;}
1010         std::vector<EcalRecHit> RecHitsInWindow;
1011         std::vector<EcalRecHit> RecHitsInWindow5x5;
1012 
1013         double simple_energy = 0;
1014 
1015         for (std::vector<DetId>::iterator det = clus_v.begin(); det != clus_v.end(); det++) {
1016           EBDetId EBdet = *det;
1017           //      cout<<" det "<< EBdet<<" ieta "<<EBdet.ieta()<<" iphi
1018           //      "<<EBdet.iphi()<<endl;
1019           bool HitAlreadyUsed = false;
1020           for (usedIds = usedXtals.begin(); usedIds != usedXtals.end(); usedIds++) {
1021             if (*usedIds == *det) {
1022               HitAlreadyUsed = true;
1023               break;
1024             }
1025           }
1026           if (HitAlreadyUsed)
1027             continue;
1028 
1029           std::vector<EBDetId>::iterator itdet = find(detIdEBRecHits.begin(), detIdEBRecHits.end(), EBdet);
1030           if (itdet == detIdEBRecHits.end())
1031             continue;
1032 
1033           int nn = int(itdet - detIdEBRecHits.begin());
1034           usedXtals.push_back(*det);
1035           RecHitsInWindow.push_back(EBRecHits[nn]);
1036           clus_used.push_back(std::make_pair(*det, 1));
1037           simple_energy = simple_energy + EBRecHits[nn].energy();
1038         }
1039 
1040         if (simple_energy <= 0)
1041           continue;
1042 
1043         math::XYZPoint clus_pos =
1044             posCalculator_.Calculate_Location(clus_used, hitCollection_p, geometry_p, geometryES_p);
1045         // cout<< "       Simple Clustering: Total energy for this simple
1046         // cluster : "<<simple_energy<<endl; cout<< "       Simple Clustering:
1047         // eta phi : "<<clus_pos.eta()<<" "<<clus_pos.phi()<<endl; cout<< "
1048         // Simple Clustering: x y z : "<<clus_pos.x()<<" "<<clus_pos.y()<<"
1049         // "<<clus_pos.z()<<endl;
1050 
1051         float theta_s = 2. * atan(exp(-clus_pos.eta()));
1052         //  float p0x_s = simple_energy * sin(theta_s) *
1053         // cos(clus_pos.phi()); float p0y_s = simple_energy * sin(theta_s) *
1054         // sin(clus_pos.phi());
1055         //    float p0z_s = simple_energy * cos(theta_s);
1056         // float et_s = sqrt( p0x_s*p0x_s + p0y_s*p0y_s);
1057 
1058         float et_s = simple_energy * sin(theta_s);
1059         // cout << "       Simple Clustering: E,Et,px,py,pz: "<<simple_energy<<"
1060         // "<<et_s<<" "<<p0x_s<<" "<<p0y_s<<" "<<endl;
1061 
1062         // Compute S4/S9 variable
1063         // We are not sure to have 9 RecHits so need to check eta and phi:
1064         /// check s4s9
1065         float s4s9_tmp[4];
1066         for (int i = 0; i < 4; i++)
1067           s4s9_tmp[i] = 0;
1068 
1069         int seed_ieta = seed_id.ieta();
1070         int seed_iphi = seed_id.iphi();
1071 
1072         convxtalid(seed_iphi, seed_ieta);
1073 
1074         float e3x3 = 0;
1075         float e5x5 = 0;
1076 
1077         for (unsigned int j = 0; j < RecHitsInWindow.size(); j++) {
1078           EBDetId det = (EBDetId)RecHitsInWindow[j].id();
1079 
1080           int ieta = det.ieta();
1081           int iphi = det.iphi();
1082 
1083           convxtalid(iphi, ieta);
1084 
1085           float en = RecHitsInWindow[j].energy();
1086 
1087           int dx = diff_neta_s(seed_ieta, ieta);
1088           int dy = diff_nphi_s(seed_iphi, iphi);
1089 
1090           if (dx <= 0 && dy <= 0)
1091             s4s9_tmp[0] += en;
1092           if (dx >= 0 && dy <= 0)
1093             s4s9_tmp[1] += en;
1094           if (dx <= 0 && dy >= 0)
1095             s4s9_tmp[2] += en;
1096           if (dx >= 0 && dy >= 0)
1097             s4s9_tmp[3] += en;
1098 
1099           if (std::abs(dx) <= 1 && std::abs(dy) <= 1)
1100             e3x3 += en;
1101           if (std::abs(dx) <= 2 && std::abs(dy) <= 2)
1102             e5x5 += en;
1103         }
1104 
1105         if (e3x3 <= 0)
1106           continue;
1107 
1108         float s4s9_max = *std::max_element(s4s9_tmp, s4s9_tmp + 4) / e3x3;
1109 
1110         /// calculate e5x5
1111         std::vector<DetId> clus_v5x5 = topology_p->getWindow(seed_id, 5, 5);
1112         for (std::vector<DetId>::const_iterator idItr = clus_v5x5.begin(); idItr != clus_v5x5.end(); idItr++) {
1113           EBDetId det = *idItr;
1114 
1115           std::vector<EBDetId>::iterator itdet0 = find(usedXtals.begin(), usedXtals.end(), det);
1116 
1117           /// already clustered
1118           if (itdet0 != usedXtals.end())
1119             continue;
1120 
1121           // inside collections
1122           std::vector<EBDetId>::iterator itdet = find(detIdEBRecHits.begin(), detIdEBRecHits.end(), det);
1123           if (itdet == detIdEBRecHits.end())
1124             continue;
1125 
1126           int nn = int(itdet - detIdEBRecHits.begin());
1127 
1128           RecHitsInWindow5x5.push_back(EBRecHits[nn]);
1129           e5x5 += EBRecHits[nn].energy();
1130         }
1131 
1132         if (e5x5 <= 0)
1133           continue;
1134 
1135         eClus.push_back(simple_energy);
1136         etClus.push_back(et_s);
1137         etaClus.push_back(clus_pos.eta());
1138         thetaClus.push_back(theta_s);
1139         phiClus.push_back(clus_pos.phi());
1140         s4s9Clus.push_back(s4s9_max);
1141         s9s25Clus.push_back(e3x3 / e5x5);
1142         RecHitsCluster.push_back(RecHitsInWindow);
1143         RecHitsCluster5x5.push_back(RecHitsInWindow5x5);
1144 
1145         //      std::cout<<" EB Eta cluster (n,nxt,e,et eta,phi,s4s9)
1146         //"<<nClus<<" "<<int(RecHitsInWindow.size())<<" "<<eClus[nClus]<<" "<<"
1147         //"<<etClus[nClus]<<" "<<etaClus[nClus]<<" "<<phiClus[nClus]<<"
1148         //"<<s4s9Clus[nClus]<<std::endl;
1149 
1150         nClus++;
1151       }
1152 
1153       // cout<< " Eta clusters: "<<nClus<<endl;
1154 
1155       // Selection, based on Simple clustering
1156       // eta candidates
1157 
1158       //      if (nClus <= 1) return;
1159       for (Int_t i = 0; i < nClus; i++) {
1160         for (Int_t j = i + 1; j < nClus; j++) {
1161           //      cout<<" i "<<i<<"  etClus[i] "<<etClus[i]<<" j "<<j<<"
1162           //      etClus[j] "<<etClus[j]<<endl;
1163           if (etClus[i] > selePtGammaEta_ && etClus[j] > selePtGammaEta_ && s4s9Clus[i] > seleS4S9GammaEta_ &&
1164               s4s9Clus[j] > seleS4S9GammaEta_) {
1165             float p0x = etClus[i] * cos(phiClus[i]);
1166             float p1x = etClus[j] * cos(phiClus[j]);
1167             float p0y = etClus[i] * sin(phiClus[i]);
1168             float p1y = etClus[j] * sin(phiClus[j]);
1169             float p0z = eClus[i] * cos(thetaClus[i]);
1170             float p1z = eClus[j] * cos(thetaClus[j]);
1171 
1172             float pt_pair = sqrt((p0x + p1x) * (p0x + p1x) + (p0y + p1y) * (p0y + p1y));
1173 
1174             if (pt_pair < selePtEta_)
1175               continue;
1176 
1177             float m_inv = sqrt((eClus[i] + eClus[j]) * (eClus[i] + eClus[j]) - (p0x + p1x) * (p0x + p1x) -
1178                                (p0y + p1y) * (p0y + p1y) - (p0z + p1z) * (p0z + p1z));
1179             if ((m_inv < seleMinvMaxEta_) && (m_inv > seleMinvMinEta_)) {
1180               // New Loop on cluster to measure isolation:
1181               std::vector<int> IsoClus;
1182               IsoClus.clear();
1183               float Iso = 0;
1184               TVector3 pairVect = TVector3((p0x + p1x), (p0y + p1y), (p0z + p1z));
1185               for (Int_t k = 0; k < nClus; k++) {
1186                 if (etClus[k] < ptMinForIsolationEta_)
1187                   continue;
1188 
1189                 if (k == i || k == j)
1190                   continue;
1191                 TVector3 ClusVect =
1192                     TVector3(etClus[k] * cos(phiClus[k]), etClus[k] * sin(phiClus[k]), eClus[k] * cos(thetaClus[k]));
1193 
1194                 float dretacl = fabs(etaClus[k] - pairVect.Eta());
1195                 float drcl = ClusVect.DeltaR(pairVect);
1196                 //      cout<< "   Iso: k, E, drclpi0, detaclpi0, dphiclpi0
1197                 //      "<<k<<" "<<eClus[k]<<" "<<drclpi0<<"
1198                 //      "<<dretaclpi0<<endl;
1199                 if ((drcl < seleEtaBeltDR_) && (dretacl < seleEtaBeltDeta_)) {
1200                   //              cout<< "   ... good iso cluster #: "<<k<<"
1201                   //              etClus[k] "<<etClus[k] <<endl;
1202                   Iso = Iso + etClus[k];
1203                   IsoClus.push_back(k);
1204                 }
1205               }
1206 
1207               //      cout<<"  Iso/pt_pi0 "<<Iso/pt_pi0<<endl;
1208               if (Iso / pt_pair < seleEtaIso_) {
1209                 // for(unsigned int
1210                 // Rec=0;Rec<RecHitsCluster[i].size();Rec++)pi0EBRecHitCollection->push_back(RecHitsCluster[i][Rec]);
1211                 // for(unsigned int
1212                 // Rec2=0;Rec2<RecHitsCluster[j].size();Rec2++)pi0EBRecHitCollection->push_back(RecHitsCluster[j][Rec2]);
1213 
1214                 hMinvEtaEB_->Fill(m_inv);
1215                 hPt1EtaEB_->Fill(etClus[i]);
1216                 hPt2EtaEB_->Fill(etClus[j]);
1217                 hPtEtaEB_->Fill(pt_pair);
1218                 hIsoEtaEB_->Fill(Iso / pt_pair);
1219                 hS4S91EtaEB_->Fill(s4s9Clus[i]);
1220                 hS4S92EtaEB_->Fill(s4s9Clus[j]);
1221 
1222                 //      cout <<"  EB Simple Clustering: Eta Candidate
1223                 // pt, eta, phi, Iso, m_inv, i, j :   "<<pt_pair<<"
1224                 //"<<pairVect.Eta()<<" "<<pairVect.Phi()<<" "<<Iso<<"
1225                 //"<<m_inv<<" "<<i<<" "<<j<<" "<<endl;
1226               }
1227             }
1228           }
1229         }  // End of the "j" loop over Simple Clusters
1230       }  // End of the "i" loop over Simple Clusters
1231 
1232     }  // rhEBeta.valid() ends
1233 
1234   }  //  isMonEBeta ends
1235 
1236   //------------------ End of Eta in EB --------------------------//
1237 
1238   //----------------- End of the EB --------------------------//
1239 
1240   //----------------- Start of the EE --------------------//
1241 
1242   // fill pi0 EE histos
1243   if (isMonEEpi0_) {
1244     if (rhEEpi0.isValid() && (!rhEEpi0->empty())) {
1245       const EcalRecHitCollection *hitCollection_ee = rhEEpi0.product();
1246       float etot = 0;
1247 
1248       detIdEERecHits.clear();  //// EEDetId
1249       EERecHits.clear();       /// EcalRecHit
1250 
1251       std::vector<EcalRecHit> seedsEndCap;
1252       seedsEndCap.clear();
1253 
1254       std::vector<EEDetId> usedXtalsEndCap;
1255       usedXtalsEndCap.clear();
1256 
1257       ////make seeds.
1258       EERecHitCollection::const_iterator ite;
1259       for (ite = rhEEpi0->begin(); ite != rhEEpi0->end(); ite++) {
1260         double energy = ite->energy();
1261         if (energy < seleXtalMinEnergyEndCap_)
1262           continue;
1263 
1264         EEDetId det = ite->id();
1265         EEDetId id(ite->id());
1266 
1267         detIdEERecHits.push_back(det);
1268         EERecHits.push_back(*ite);
1269 
1270         hiXDistrEEpi0_->Fill(id.ix());
1271         hiYDistrEEpi0_->Fill(id.iy());
1272         hRechitEnergyEEpi0_->Fill(ite->energy());
1273 
1274         if (energy > clusSeedThrEndCap_)
1275           seedsEndCap.push_back(*ite);
1276 
1277         etot += ite->energy();
1278       }  // EE rechits
1279 
1280       hNRecHitsEEpi0_->Fill(rhEEpi0->size());
1281       hMeanRecHitEnergyEEpi0_->Fill(etot / rhEEpi0->size());
1282       hEventEnergyEEpi0_->Fill(etot);
1283 
1284       //      cout << " EE RH Pi0 collection: #, mean rh_e, event E
1285       //      "<<rhEEpi0->size()<<" "<<etot/rhEEpi0->size()<<" "<<etot<<endl;
1286 
1287       int nClusEndCap;
1288       std::vector<float> eClusEndCap;
1289       std::vector<float> etClusEndCap;
1290       std::vector<float> etaClusEndCap;
1291       std::vector<float> thetaClusEndCap;
1292       std::vector<float> phiClusEndCap;
1293       std::vector<std::vector<EcalRecHit>> RecHitsClusterEndCap;
1294       std::vector<std::vector<EcalRecHit>> RecHitsCluster5x5EndCap;
1295       std::vector<float> s4s9ClusEndCap;
1296       std::vector<float> s9s25ClusEndCap;
1297 
1298       nClusEndCap = 0;
1299 
1300       // Make own simple clusters (3x3, 5x5 or clusPhiSize_ x clusEtaSize_)
1301       std::sort(seedsEndCap.begin(), seedsEndCap.end(), ecalRecHitGreater);
1302 
1303       for (std::vector<EcalRecHit>::iterator itseed = seedsEndCap.begin(); itseed != seedsEndCap.end(); itseed++) {
1304         EEDetId seed_id = itseed->id();
1305         std::vector<EEDetId>::const_iterator usedIds;
1306 
1307         bool seedAlreadyUsed = false;
1308         for (usedIds = usedXtalsEndCap.begin(); usedIds != usedXtalsEndCap.end(); usedIds++) {
1309           if (*usedIds == seed_id) {
1310             seedAlreadyUsed = true;
1311             break;
1312           }
1313         }
1314 
1315         if (seedAlreadyUsed)
1316           continue;
1317         std::vector<DetId> clus_v = topology_ee->getWindow(seed_id, clusEtaSize_, clusPhiSize_);
1318         std::vector<std::pair<DetId, float>> clus_used;
1319 
1320         std::vector<EcalRecHit> RecHitsInWindow;
1321         std::vector<EcalRecHit> RecHitsInWindow5x5;
1322 
1323         float simple_energy = 0;
1324 
1325         for (std::vector<DetId>::iterator det = clus_v.begin(); det != clus_v.end(); det++) {
1326           EEDetId EEdet = *det;
1327 
1328           bool HitAlreadyUsed = false;
1329           for (usedIds = usedXtalsEndCap.begin(); usedIds != usedXtalsEndCap.end(); usedIds++) {
1330             if (*usedIds == *det) {
1331               HitAlreadyUsed = true;
1332               break;
1333             }
1334           }
1335 
1336           if (HitAlreadyUsed)
1337             continue;
1338 
1339           std::vector<EEDetId>::iterator itdet = find(detIdEERecHits.begin(), detIdEERecHits.end(), EEdet);
1340           if (itdet == detIdEERecHits.end())
1341             continue;
1342 
1343           int nn = int(itdet - detIdEERecHits.begin());
1344           usedXtalsEndCap.push_back(*det);
1345           RecHitsInWindow.push_back(EERecHits[nn]);
1346           clus_used.push_back(std::make_pair(*det, 1));
1347           simple_energy = simple_energy + EERecHits[nn].energy();
1348         }
1349 
1350         if (simple_energy <= 0)
1351           continue;
1352 
1353         math::XYZPoint clus_pos =
1354             posCalculator_.Calculate_Location(clus_used, hitCollection_ee, geometryEE_p, geometryES_p);
1355 
1356         float theta_s = 2. * atan(exp(-clus_pos.eta()));
1357         float et_s = simple_energy * sin(theta_s);
1358         //  float p0x_s = simple_energy * sin(theta_s) *
1359         // cos(clus_pos.phi()); float p0y_s = simple_energy * sin(theta_s) *
1360         // sin(clus_pos.phi()); float et_s = sqrt( p0x_s*p0x_s + p0y_s*p0y_s);
1361 
1362         // Compute S4/S9 variable
1363         // We are not sure to have 9 RecHits so need to check eta and phi:
1364         float s4s9_tmp[4];
1365         for (int i = 0; i < 4; i++)
1366           s4s9_tmp[i] = 0;
1367 
1368         int ixSeed = seed_id.ix();
1369         int iySeed = seed_id.iy();
1370         float e3x3 = 0;
1371         float e5x5 = 0;
1372 
1373         for (unsigned int j = 0; j < RecHitsInWindow.size(); j++) {
1374           EEDetId det_this = (EEDetId)RecHitsInWindow[j].id();
1375           int dx = ixSeed - det_this.ix();
1376           int dy = iySeed - det_this.iy();
1377 
1378           float en = RecHitsInWindow[j].energy();
1379 
1380           if (dx <= 0 && dy <= 0)
1381             s4s9_tmp[0] += en;
1382           if (dx >= 0 && dy <= 0)
1383             s4s9_tmp[1] += en;
1384           if (dx <= 0 && dy >= 0)
1385             s4s9_tmp[2] += en;
1386           if (dx >= 0 && dy >= 0)
1387             s4s9_tmp[3] += en;
1388 
1389           if (std::abs(dx) <= 1 && std::abs(dy) <= 1)
1390             e3x3 += en;
1391           if (std::abs(dx) <= 2 && std::abs(dy) <= 2)
1392             e5x5 += en;
1393         }
1394 
1395         if (e3x3 <= 0)
1396           continue;
1397 
1398         eClusEndCap.push_back(simple_energy);
1399         etClusEndCap.push_back(et_s);
1400         etaClusEndCap.push_back(clus_pos.eta());
1401         thetaClusEndCap.push_back(theta_s);
1402         phiClusEndCap.push_back(clus_pos.phi());
1403         s4s9ClusEndCap.push_back(*std::max_element(s4s9_tmp, s4s9_tmp + 4) / e3x3);
1404         s9s25ClusEndCap.push_back(e3x3 / e5x5);
1405         RecHitsClusterEndCap.push_back(RecHitsInWindow);
1406         RecHitsCluster5x5EndCap.push_back(RecHitsInWindow5x5);
1407 
1408         //  std::cout<<" EE pi0 cluster (n,nxt,e,et eta,phi,s4s9)
1409         //"<<nClusEndCap<<" "<<int(RecHitsInWindow.size())<<"
1410         //"<<eClusEndCap[nClusEndCap]<<" "<<" "<<etClusEndCap[nClusEndCap]<<"
1411         //"<<etaClusEndCap[nClusEndCap]<<" "<<phiClusEndCap[nClusEndCap]<<"
1412         //"<<s4s9ClusEndCap[nClusEndCap]<<std::endl;
1413 
1414         nClusEndCap++;
1415       }
1416 
1417       // Selection, based on Simple clustering
1418       // pi0 candidates
1419 
1420       for (Int_t i = 0; i < nClusEndCap; i++) {
1421         for (Int_t j = i + 1; j < nClusEndCap; j++) {
1422           if (etClusEndCap[i] > selePtGammaEndCap_ && etClusEndCap[j] > selePtGammaEndCap_ &&
1423               s4s9ClusEndCap[i] > seleS4S9GammaEndCap_ && s4s9ClusEndCap[j] > seleS4S9GammaEndCap_) {
1424             float p0x = etClusEndCap[i] * cos(phiClusEndCap[i]);
1425             float p1x = etClusEndCap[j] * cos(phiClusEndCap[j]);
1426             float p0y = etClusEndCap[i] * sin(phiClusEndCap[i]);
1427             float p1y = etClusEndCap[j] * sin(phiClusEndCap[j]);
1428             float p0z = eClusEndCap[i] * cos(thetaClusEndCap[i]);
1429             float p1z = eClusEndCap[j] * cos(thetaClusEndCap[j]);
1430 
1431             float pt_pair = sqrt((p0x + p1x) * (p0x + p1x) + (p0y + p1y) * (p0y + p1y));
1432             if (pt_pair < selePtPi0EndCap_)
1433               continue;
1434             float m_inv = sqrt((eClusEndCap[i] + eClusEndCap[j]) * (eClusEndCap[i] + eClusEndCap[j]) -
1435                                (p0x + p1x) * (p0x + p1x) - (p0y + p1y) * (p0y + p1y) - (p0z + p1z) * (p0z + p1z));
1436 
1437             if ((m_inv < seleMinvMaxPi0EndCap_) && (m_inv > seleMinvMinPi0EndCap_)) {
1438               // New Loop on cluster to measure isolation:
1439               std::vector<int> IsoClus;
1440               IsoClus.clear();
1441               float Iso = 0;
1442               TVector3 pairVect = TVector3((p0x + p1x), (p0y + p1y), (p0z + p1z));
1443               for (Int_t k = 0; k < nClusEndCap; k++) {
1444                 if (etClusEndCap[k] < ptMinForIsolationEndCap_)
1445                   continue;
1446 
1447                 if (k == i || k == j)
1448                   continue;
1449 
1450                 TVector3 clusVect = TVector3(etClusEndCap[k] * cos(phiClusEndCap[k]),
1451                                              etClusEndCap[k] * sin(phiClusEndCap[k]),
1452                                              eClusEndCap[k] * cos(thetaClusEndCap[k]));
1453                 float dretacl = fabs(etaClusEndCap[k] - pairVect.Eta());
1454                 float drcl = clusVect.DeltaR(pairVect);
1455 
1456                 if (drcl < selePi0BeltDREndCap_ && dretacl < selePi0BeltDetaEndCap_) {
1457                   Iso = Iso + etClusEndCap[k];
1458                   IsoClus.push_back(k);
1459                 }
1460               }
1461 
1462               if (Iso / pt_pair < selePi0IsoEndCap_) {
1463                 // cout <<"  EE Simple Clustering: pi0 Candidate pt, eta, phi,
1464                 // Iso, m_inv, i, j :   "<<pt_pair<<" "<<pairVect.Eta()<<"
1465                 // "<<pairVect.Phi()<<" "<<Iso<<" "<<m_inv<<" "<<i<<" "<<j<<"
1466                 // "<<endl;
1467 
1468                 hMinvPi0EE_->Fill(m_inv);
1469                 hPt1Pi0EE_->Fill(etClusEndCap[i]);
1470                 hPt2Pi0EE_->Fill(etClusEndCap[j]);
1471                 hPtPi0EE_->Fill(pt_pair);
1472                 hIsoPi0EE_->Fill(Iso / pt_pair);
1473                 hS4S91Pi0EE_->Fill(s4s9ClusEndCap[i]);
1474                 hS4S92Pi0EE_->Fill(s4s9ClusEndCap[j]);
1475               }
1476             }
1477           }
1478         }  // End of the "j" loop over Simple Clusters
1479       }  // End of the "i" loop over Simple Clusters
1480 
1481     }  // rhEEpi0
1482   }  // isMonEEpi0
1483 
1484   //================End of Pi0 endcap=======================//
1485 
1486   //================ Eta in EE===============================//
1487 
1488   // fill pi0 EE histos
1489   if (isMonEEeta_) {
1490     if (rhEEeta.isValid() && (!rhEEeta->empty())) {
1491       const EcalRecHitCollection *hitCollection_ee = rhEEeta.product();
1492       float etot = 0;
1493 
1494       detIdEERecHits.clear();  //// EEDetId
1495       EERecHits.clear();       /// EcalRecHit
1496 
1497       std::vector<EcalRecHit> seedsEndCap;
1498       seedsEndCap.clear();
1499 
1500       std::vector<EEDetId> usedXtalsEndCap;
1501       usedXtalsEndCap.clear();
1502 
1503       ////make seeds.
1504       EERecHitCollection::const_iterator ite;
1505       for (ite = rhEEeta->begin(); ite != rhEEeta->end(); ite++) {
1506         double energy = ite->energy();
1507         if (energy < seleXtalMinEnergyEndCap_)
1508           continue;
1509 
1510         EEDetId det = ite->id();
1511         EEDetId id(ite->id());
1512 
1513         detIdEERecHits.push_back(det);
1514         EERecHits.push_back(*ite);
1515 
1516         hiXDistrEEeta_->Fill(id.ix());
1517         hiYDistrEEeta_->Fill(id.iy());
1518         hRechitEnergyEEeta_->Fill(ite->energy());
1519 
1520         if (energy > clusSeedThrEndCap_)
1521           seedsEndCap.push_back(*ite);
1522 
1523         etot += ite->energy();
1524       }  // EE rechits
1525 
1526       hNRecHitsEEeta_->Fill(rhEEeta->size());
1527       hMeanRecHitEnergyEEeta_->Fill(etot / rhEEeta->size());
1528       hEventEnergyEEeta_->Fill(etot);
1529 
1530       //      cout << " EE RH Eta collection: #, mean rh_e, event E
1531       //      "<<rhEEeta->size()<<" "<<etot/rhEEeta->size()<<" "<<etot<<endl;
1532 
1533       int nClusEndCap;
1534       std::vector<float> eClusEndCap;
1535       std::vector<float> etClusEndCap;
1536       std::vector<float> etaClusEndCap;
1537       std::vector<float> thetaClusEndCap;
1538       std::vector<float> phiClusEndCap;
1539       std::vector<std::vector<EcalRecHit>> RecHitsClusterEndCap;
1540       std::vector<std::vector<EcalRecHit>> RecHitsCluster5x5EndCap;
1541       std::vector<float> s4s9ClusEndCap;
1542       std::vector<float> s9s25ClusEndCap;
1543 
1544       nClusEndCap = 0;
1545 
1546       // Make own simple clusters (3x3, 5x5 or clusPhiSize_ x clusEtaSize_)
1547       std::sort(seedsEndCap.begin(), seedsEndCap.end(), ecalRecHitGreater);
1548 
1549       for (std::vector<EcalRecHit>::iterator itseed = seedsEndCap.begin(); itseed != seedsEndCap.end(); itseed++) {
1550         EEDetId seed_id = itseed->id();
1551         std::vector<EEDetId>::const_iterator usedIds;
1552 
1553         bool seedAlreadyUsed = false;
1554         for (usedIds = usedXtalsEndCap.begin(); usedIds != usedXtalsEndCap.end(); usedIds++) {
1555           if (*usedIds == seed_id) {
1556             seedAlreadyUsed = true;
1557             break;
1558           }
1559         }
1560 
1561         if (seedAlreadyUsed)
1562           continue;
1563         std::vector<DetId> clus_v = topology_ee->getWindow(seed_id, clusEtaSize_, clusPhiSize_);
1564         std::vector<std::pair<DetId, float>> clus_used;
1565 
1566         std::vector<EcalRecHit> RecHitsInWindow;
1567         std::vector<EcalRecHit> RecHitsInWindow5x5;
1568 
1569         float simple_energy = 0;
1570 
1571         for (std::vector<DetId>::iterator det = clus_v.begin(); det != clus_v.end(); det++) {
1572           EEDetId EEdet = *det;
1573 
1574           bool HitAlreadyUsed = false;
1575           for (usedIds = usedXtalsEndCap.begin(); usedIds != usedXtalsEndCap.end(); usedIds++) {
1576             if (*usedIds == *det) {
1577               HitAlreadyUsed = true;
1578               break;
1579             }
1580           }
1581 
1582           if (HitAlreadyUsed)
1583             continue;
1584 
1585           std::vector<EEDetId>::iterator itdet = find(detIdEERecHits.begin(), detIdEERecHits.end(), EEdet);
1586           if (itdet == detIdEERecHits.end())
1587             continue;
1588 
1589           int nn = int(itdet - detIdEERecHits.begin());
1590           usedXtalsEndCap.push_back(*det);
1591           RecHitsInWindow.push_back(EERecHits[nn]);
1592           clus_used.push_back(std::make_pair(*det, 1));
1593           simple_energy = simple_energy + EERecHits[nn].energy();
1594         }
1595 
1596         if (simple_energy <= 0)
1597           continue;
1598 
1599         math::XYZPoint clus_pos =
1600             posCalculator_.Calculate_Location(clus_used, hitCollection_ee, geometryEE_p, geometryES_p);
1601 
1602         float theta_s = 2. * atan(exp(-clus_pos.eta()));
1603         float et_s = simple_energy * sin(theta_s);
1604         //  float p0x_s = simple_energy * sin(theta_s) *
1605         // cos(clus_pos.phi()); float p0y_s = simple_energy * sin(theta_s) *
1606         // sin(clus_pos.phi()); float et_s = sqrt( p0x_s*p0x_s + p0y_s*p0y_s);
1607 
1608         // Compute S4/S9 variable
1609         // We are not sure to have 9 RecHits so need to check eta and phi:
1610         float s4s9_tmp[4];
1611         for (int i = 0; i < 4; i++)
1612           s4s9_tmp[i] = 0;
1613 
1614         int ixSeed = seed_id.ix();
1615         int iySeed = seed_id.iy();
1616         float e3x3 = 0;
1617         float e5x5 = 0;
1618 
1619         for (unsigned int j = 0; j < RecHitsInWindow.size(); j++) {
1620           EEDetId det_this = (EEDetId)RecHitsInWindow[j].id();
1621           int dx = ixSeed - det_this.ix();
1622           int dy = iySeed - det_this.iy();
1623 
1624           float en = RecHitsInWindow[j].energy();
1625 
1626           if (dx <= 0 && dy <= 0)
1627             s4s9_tmp[0] += en;
1628           if (dx >= 0 && dy <= 0)
1629             s4s9_tmp[1] += en;
1630           if (dx <= 0 && dy >= 0)
1631             s4s9_tmp[2] += en;
1632           if (dx >= 0 && dy >= 0)
1633             s4s9_tmp[3] += en;
1634 
1635           if (std::abs(dx) <= 1 && std::abs(dy) <= 1)
1636             e3x3 += en;
1637           if (std::abs(dx) <= 2 && std::abs(dy) <= 2)
1638             e5x5 += en;
1639         }
1640 
1641         if (e3x3 <= 0)
1642           continue;
1643 
1644         eClusEndCap.push_back(simple_energy);
1645         etClusEndCap.push_back(et_s);
1646         etaClusEndCap.push_back(clus_pos.eta());
1647         thetaClusEndCap.push_back(theta_s);
1648         phiClusEndCap.push_back(clus_pos.phi());
1649         s4s9ClusEndCap.push_back(*std::max_element(s4s9_tmp, s4s9_tmp + 4) / e3x3);
1650         s9s25ClusEndCap.push_back(e3x3 / e5x5);
1651         RecHitsClusterEndCap.push_back(RecHitsInWindow);
1652         RecHitsCluster5x5EndCap.push_back(RecHitsInWindow5x5);
1653 
1654         //  std::cout<<" EE Eta cluster (n,nxt,e,et eta,phi,s4s9)
1655         //"<<nClusEndCap<<" "<<int(RecHitsInWindow.size())<<"
1656         //"<<eClusEndCap[nClusEndCap]<<" "<<" "<<etClusEndCap[nClusEndCap]<<"
1657         //"<<etaClusEndCap[nClusEndCap]<<" "<<phiClusEndCap[nClusEndCap]<<"
1658         //"<<s4s9ClusEndCap[nClusEndCap]<<std::endl;
1659 
1660         nClusEndCap++;
1661       }
1662 
1663       // Selection, based on Simple clustering
1664       // pi0 candidates
1665 
1666       for (Int_t i = 0; i < nClusEndCap; i++) {
1667         for (Int_t j = i + 1; j < nClusEndCap; j++) {
1668           if (etClusEndCap[i] > selePtGammaEtaEndCap_ && etClusEndCap[j] > selePtGammaEtaEndCap_ &&
1669               s4s9ClusEndCap[i] > seleS4S9GammaEtaEndCap_ && s4s9ClusEndCap[j] > seleS4S9GammaEtaEndCap_) {
1670             float p0x = etClusEndCap[i] * cos(phiClusEndCap[i]);
1671             float p1x = etClusEndCap[j] * cos(phiClusEndCap[j]);
1672             float p0y = etClusEndCap[i] * sin(phiClusEndCap[i]);
1673             float p1y = etClusEndCap[j] * sin(phiClusEndCap[j]);
1674             float p0z = eClusEndCap[i] * cos(thetaClusEndCap[i]);
1675             float p1z = eClusEndCap[j] * cos(thetaClusEndCap[j]);
1676 
1677             float pt_pair = sqrt((p0x + p1x) * (p0x + p1x) + (p0y + p1y) * (p0y + p1y));
1678             if (pt_pair < selePtEtaEndCap_)
1679               continue;
1680             float m_inv = sqrt((eClusEndCap[i] + eClusEndCap[j]) * (eClusEndCap[i] + eClusEndCap[j]) -
1681                                (p0x + p1x) * (p0x + p1x) - (p0y + p1y) * (p0y + p1y) - (p0z + p1z) * (p0z + p1z));
1682 
1683             if ((m_inv < seleMinvMaxEtaEndCap_) && (m_inv > seleMinvMinEtaEndCap_)) {
1684               // New Loop on cluster to measure isolation:
1685               std::vector<int> IsoClus;
1686               IsoClus.clear();
1687               float Iso = 0;
1688               TVector3 pairVect = TVector3((p0x + p1x), (p0y + p1y), (p0z + p1z));
1689               for (Int_t k = 0; k < nClusEndCap; k++) {
1690                 if (etClusEndCap[k] < ptMinForIsolationEtaEndCap_)
1691                   continue;
1692 
1693                 if (k == i || k == j)
1694                   continue;
1695 
1696                 TVector3 clusVect = TVector3(etClusEndCap[k] * cos(phiClusEndCap[k]),
1697                                              etClusEndCap[k] * sin(phiClusEndCap[k]),
1698                                              eClusEndCap[k] * cos(thetaClusEndCap[k]));
1699                 float dretacl = fabs(etaClusEndCap[k] - pairVect.Eta());
1700                 float drcl = clusVect.DeltaR(pairVect);
1701 
1702                 if (drcl < seleEtaBeltDREndCap_ && dretacl < seleEtaBeltDetaEndCap_) {
1703                   Iso = Iso + etClusEndCap[k];
1704                   IsoClus.push_back(k);
1705                 }
1706               }
1707 
1708               if (Iso / pt_pair < seleEtaIsoEndCap_) {
1709                 //  cout <<"  EE Simple Clustering: Eta Candidate pt, eta,
1710                 // phi, Iso, m_inv, i, j :   "<<pt_pair<<" "<<pairVect.Eta()<<"
1711                 //"<<pairVect.Phi()<<" "<<Iso<<" "<<m_inv<<" "<<i<<" "<<j<<"
1712                 //"<<endl;
1713 
1714                 hMinvEtaEE_->Fill(m_inv);
1715                 hPt1EtaEE_->Fill(etClusEndCap[i]);
1716                 hPt2EtaEE_->Fill(etClusEndCap[j]);
1717                 hPtEtaEE_->Fill(pt_pair);
1718                 hIsoEtaEE_->Fill(Iso / pt_pair);
1719                 hS4S91EtaEE_->Fill(s4s9ClusEndCap[i]);
1720                 hS4S92EtaEE_->Fill(s4s9ClusEndCap[j]);
1721               }
1722             }
1723           }
1724         }  // End of the "j" loop over Simple Clusters
1725       }  // End of the "i" loop over Simple Clusters
1726 
1727     }  // rhEEeta
1728   }  // isMonEEeta
1729 
1730   //================End of Pi0 endcap=======================//
1731 
1732   ////==============End of endcap ===============///
1733 }
1734 
1735 void DQMSourcePi0::convxtalid(Int_t &nphi, Int_t &neta) {
1736   // Barrel only
1737   // Output nphi 0...359; neta 0...84; nside=+1 (for eta>0), or 0 (for eta<0).
1738   // neta will be [-85,-1] , or [0,84], the minus sign indicates the z<0 side.
1739 
1740   if (neta > 0)
1741     neta -= 1;
1742   if (nphi > 359)
1743     nphi = nphi - 360;
1744 
1745 }  // end of convxtalid
1746 
1747 int DQMSourcePi0::diff_neta_s(Int_t neta1, Int_t neta2) {
1748   Int_t mdiff;
1749   mdiff = (neta1 - neta2);
1750   return mdiff;
1751 }
1752 
1753 // Calculate the distance in xtals taking into account the periodicity of the
1754 // Barrel
1755 int DQMSourcePi0::diff_nphi_s(Int_t nphi1, Int_t nphi2) {
1756   Int_t mdiff;
1757   if (std::abs(nphi1 - nphi2) < (360 - std::abs(nphi1 - nphi2))) {
1758     mdiff = nphi1 - nphi2;
1759   } else {
1760     mdiff = 360 - std::abs(nphi1 - nphi2);
1761     if (nphi1 > nphi2)
1762       mdiff = -mdiff;
1763   }
1764   return mdiff;
1765 }
1766 
1767 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!