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
0011
0012 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0013 #include "DQMServices/Core/interface/DQMStore.h"
0014
0015
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
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
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
0062 MonitorElement *hiPhiDistrEBpi0_;
0063
0064
0065 MonitorElement *hiXDistrEEpi0_;
0066
0067
0068 MonitorElement *hiPhiDistrEBeta_;
0069
0070
0071 MonitorElement *hiXDistrEEeta_;
0072
0073
0074 MonitorElement *hiEtaDistrEBpi0_;
0075
0076
0077 MonitorElement *hiYDistrEEpi0_;
0078
0079
0080 MonitorElement *hiEtaDistrEBeta_;
0081
0082
0083 MonitorElement *hiYDistrEEeta_;
0084
0085
0086 MonitorElement *hRechitEnergyEBpi0_;
0087
0088
0089 MonitorElement *hRechitEnergyEEpi0_;
0090
0091
0092 MonitorElement *hRechitEnergyEBeta_;
0093
0094
0095 MonitorElement *hRechitEnergyEEeta_;
0096
0097
0098 MonitorElement *hEventEnergyEBpi0_;
0099
0100
0101 MonitorElement *hEventEnergyEEpi0_;
0102
0103
0104 MonitorElement *hEventEnergyEBeta_;
0105
0106
0107 MonitorElement *hEventEnergyEEeta_;
0108
0109
0110 MonitorElement *hNRecHitsEBpi0_;
0111
0112
0113 MonitorElement *hNRecHitsEEpi0_;
0114
0115
0116 MonitorElement *hNRecHitsEBeta_;
0117
0118
0119 MonitorElement *hNRecHitsEEeta_;
0120
0121
0122 MonitorElement *hMeanRecHitEnergyEBpi0_;
0123
0124
0125 MonitorElement *hMeanRecHitEnergyEEpi0_;
0126
0127
0128 MonitorElement *hMeanRecHitEnergyEBeta_;
0129
0130
0131 MonitorElement *hMeanRecHitEnergyEEeta_;
0132
0133
0134 MonitorElement *hMinvPi0EB_;
0135
0136
0137 MonitorElement *hMinvPi0EE_;
0138
0139
0140 MonitorElement *hMinvEtaEB_;
0141
0142
0143 MonitorElement *hMinvEtaEE_;
0144
0145
0146 MonitorElement *hPt1Pi0EB_;
0147
0148
0149 MonitorElement *hPt1Pi0EE_;
0150
0151
0152 MonitorElement *hPt1EtaEB_;
0153
0154
0155 MonitorElement *hPt1EtaEE_;
0156
0157
0158 MonitorElement *hPt2Pi0EB_;
0159
0160
0161 MonitorElement *hPt2Pi0EE_;
0162
0163
0164 MonitorElement *hPt2EtaEB_;
0165
0166
0167 MonitorElement *hPt2EtaEE_;
0168
0169
0170 MonitorElement *hPtPi0EB_;
0171
0172
0173 MonitorElement *hPtPi0EE_;
0174
0175
0176 MonitorElement *hPtEtaEB_;
0177
0178
0179 MonitorElement *hPtEtaEE_;
0180
0181
0182 MonitorElement *hIsoPi0EB_;
0183
0184
0185 MonitorElement *hIsoPi0EE_;
0186
0187
0188 MonitorElement *hIsoEtaEB_;
0189
0190
0191 MonitorElement *hIsoEtaEE_;
0192
0193
0194 MonitorElement *hS4S91Pi0EB_;
0195
0196
0197 MonitorElement *hS4S91Pi0EE_;
0198
0199
0200 MonitorElement *hS4S91EtaEB_;
0201
0202
0203 MonitorElement *hS4S91EtaEE_;
0204
0205
0206 MonitorElement *hS4S92Pi0EB_;
0207
0208
0209 MonitorElement *hS4S92Pi0EE_;
0210
0211
0212 MonitorElement *hS4S92EtaEB_;
0213
0214
0215 MonitorElement *hS4S92EtaEE_;
0216
0217
0218 edm::EDGetTokenT<EcalRecHitCollection> productMonitoredEBpi0_;
0219 edm::EDGetTokenT<EcalRecHitCollection> productMonitoredEBeta_;
0220
0221
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
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
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
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
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
0300 unsigned int prescaleFactor_;
0301
0302
0303 std::string folderName_;
0304
0305
0306 bool saveToFile_;
0307
0308
0309 bool isMonEBpi0_;
0310 bool isMonEBeta_;
0311 bool isMonEEpi0_;
0312 bool isMonEEeta_;
0313
0314
0315 std::string fileName_;
0316 };
0317
0318 #define TWOPI 6.283185308
0319
0320
0321
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
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
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
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
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
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
0412 ibooker.setCurrentFolder(folderName_);
0413
0414
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();
0600 EBRecHits.clear();
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
0617
0618
0619
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
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 }
0655
0656 hNRecHitsEBpi0_->Fill(rhEBpi0->size());
0657 hMeanRecHitEnergyEBpi0_->Fill(etot / rhEBpi0->size());
0658 hEventEnergyEBpi0_->Fill(etot);
0659
0660
0661
0662
0663
0664
0665
0666
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
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
0695
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
0704
0705
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
0714
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
0742
0743
0744
0745
0746
0747 float theta_s = 2. * atan(exp(-clus_pos.eta()));
0748
0749
0750
0751
0752
0753
0754 float et_s = simple_energy * sin(theta_s);
0755
0756
0757
0758
0759
0760
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
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
0814 if (itdet0 != usedXtals.end())
0815 continue;
0816
0817
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
0842
0843
0844
0845
0846 nClus++;
0847 }
0848
0849
0850
0851
0852
0853 int npi0_s = 0;
0854
0855
0856 for (Int_t i = 0; i < nClus; i++) {
0857 for (Int_t j = i + 1; j < nClus; j++) {
0858
0859
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
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
0894
0895
0896 if ((drcl < selePi0BeltDR_) && (dretacl < selePi0BeltDeta_)) {
0897
0898
0899 Iso = Iso + etClus[k];
0900 IsoClus.push_back(k);
0901 }
0902 }
0903
0904
0905 if (Iso / pt_pair < selePi0Iso_) {
0906
0907
0908
0909
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
0920
0921
0922
0923
0924 npi0_s++;
0925 }
0926 }
0927 }
0928 }
0929 }
0930
0931
0932
0933 }
0934
0935 }
0936
0937
0938
0939
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 }
0964
0965 hNRecHitsEBeta_->Fill(rhEBeta->size());
0966 hMeanRecHitEnergyEBeta_->Fill(etot / rhEBeta->size());
0967 hEventEnergyEBeta_->Fill(etot);
0968
0969
0970
0971
0972
0973
0974
0975
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
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
1004
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
1013
1014
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
1023
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
1051
1052
1053
1054
1055
1056 float theta_s = 2. * atan(exp(-clus_pos.eta()));
1057
1058
1059
1060
1061
1062
1063 float et_s = simple_energy * sin(theta_s);
1064
1065
1066
1067
1068
1069
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
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
1123 if (itdet0 != usedXtals.end())
1124 continue;
1125
1126
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
1151
1152
1153
1154
1155 nClus++;
1156 }
1157
1158
1159
1160
1161
1162 int npi0_s = 0;
1163
1164
1165 for (Int_t i = 0; i < nClus; i++) {
1166 for (Int_t j = i + 1; j < nClus; j++) {
1167
1168
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
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
1203
1204
1205 if ((drcl < seleEtaBeltDR_) && (dretacl < seleEtaBeltDeta_)) {
1206
1207
1208 Iso = Iso + etClus[k];
1209 IsoClus.push_back(k);
1210 }
1211 }
1212
1213
1214 if (Iso / pt_pair < seleEtaIso_) {
1215
1216
1217
1218
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
1229
1230
1231
1232
1233 npi0_s++;
1234 }
1235 }
1236 }
1237 }
1238 }
1239
1240
1241
1242 }
1243
1244 }
1245
1246
1247
1248
1249
1250
1251
1252
1253 if (isMonEEpi0_) {
1254 if (rhEEpi0.isValid() && (!rhEEpi0->empty())) {
1255 const EcalRecHitCollection *hitCollection_ee = rhEEpi0.product();
1256 float etot = 0;
1257
1258 detIdEERecHits.clear();
1259 EERecHits.clear();
1260
1261 std::vector<EcalRecHit> seedsEndCap;
1262 seedsEndCap.clear();
1263
1264 std::vector<EEDetId> usedXtalsEndCap;
1265 usedXtalsEndCap.clear();
1266
1267
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 }
1289
1290 hNRecHitsEEpi0_->Fill(rhEEpi0->size());
1291 hMeanRecHitEnergyEEpi0_->Fill(etot / rhEEpi0->size());
1292 hEventEnergyEEpi0_->Fill(etot);
1293
1294
1295
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
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
1369
1370
1371
1372
1373
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
1419
1420
1421
1422
1423
1424 nClusEndCap++;
1425 }
1426
1427
1428
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
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
1475
1476
1477
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 }
1492 }
1493
1494
1495
1496
1497 }
1498 }
1499
1500
1501
1502
1503
1504
1505 if (isMonEEeta_) {
1506 if (rhEEeta.isValid() && (!rhEEeta->empty())) {
1507 const EcalRecHitCollection *hitCollection_ee = rhEEeta.product();
1508 float etot = 0;
1509
1510 detIdEERecHits.clear();
1511 EERecHits.clear();
1512
1513 std::vector<EcalRecHit> seedsEndCap;
1514 seedsEndCap.clear();
1515
1516 std::vector<EEDetId> usedXtalsEndCap;
1517 usedXtalsEndCap.clear();
1518
1519
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 }
1541
1542 hNRecHitsEEeta_->Fill(rhEEeta->size());
1543 hMeanRecHitEnergyEEeta_->Fill(etot / rhEEeta->size());
1544 hEventEnergyEEeta_->Fill(etot);
1545
1546
1547
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
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
1621
1622
1623
1624
1625
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
1671
1672
1673
1674
1675
1676 nClusEndCap++;
1677 }
1678
1679
1680
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
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
1727
1728
1729
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 }
1744 }
1745
1746
1747
1748
1749 }
1750 }
1751
1752
1753
1754
1755 }
1756
1757 void DQMSourcePi0::convxtalid(Int_t &nphi, Int_t &neta) {
1758
1759
1760
1761
1762 if (neta > 0)
1763 neta -= 1;
1764 if (nphi > 359)
1765 nphi = nphi - 360;
1766
1767 }
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
1776
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);