Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 #include "DQM/EcalMonitorTasks/interface/PiZeroTask.h"
0002 
0003 namespace ecaldqm {
0004   PiZeroTask::PiZeroTask()
0005       : DQWorkerTask(),
0006         seleXtalMinEnergy_(0.f),
0007         clusSeedThr_(0.f),
0008         clusEtaSize_(0),
0009         clusPhiSize_(0),
0010         selePtGammaOne_(0.f),
0011         selePtGammaTwo_(0.f),
0012         seleS4S9GammaOne_(0.f),
0013         seleS4S9GammaTwo_(0.f),
0014         selePtPi0_(0.f),
0015         selePi0Iso_(0.f),
0016         selePi0BeltDR_(0.f),
0017         selePi0BeltDeta_(0.f),
0018         seleMinvMaxPi0_(0.f),
0019         seleMinvMinPi0_(0.f),
0020         posCalcParameters_(edm::ParameterSet()) {}
0021 
0022   void PiZeroTask::setParams(edm::ParameterSet const& params) {
0023     // Parameters needed for pi0 finding
0024     seleXtalMinEnergy_ = params.getParameter<double>("seleXtalMinEnergy");
0025 
0026     clusSeedThr_ = params.getParameter<double>("clusSeedThr");
0027     clusEtaSize_ = params.getParameter<int>("clusEtaSize");
0028     clusPhiSize_ = params.getParameter<int>("clusPhiSize");
0029 
0030     selePtGammaOne_ = params.getParameter<double>("selePtGammaOne");
0031     selePtGammaTwo_ = params.getParameter<double>("selePtGammaTwo");
0032     seleS4S9GammaOne_ = params.getParameter<double>("seleS4S9GammaOne");
0033     seleS4S9GammaTwo_ = params.getParameter<double>("seleS4S9GammaTwo");
0034     selePtPi0_ = params.getParameter<double>("selePtPi0");
0035     selePi0Iso_ = params.getParameter<double>("selePi0Iso");
0036     selePi0BeltDR_ = params.getParameter<double>("selePi0BeltDR");
0037     selePi0BeltDeta_ = params.getParameter<double>("selePi0BeltDeta");
0038     seleMinvMaxPi0_ = params.getParameter<double>("seleMinvMaxPi0");
0039     seleMinvMinPi0_ = params.getParameter<double>("seleMinvMinPi0");
0040 
0041     posCalcParameters_ = params.getParameter<edm::ParameterSet>("posCalcParameters");
0042   }
0043 
0044   bool PiZeroTask::filterRunType(short const* runType) {
0045     for (unsigned iFED(0); iFED != ecaldqm::nDCC; iFED++) {
0046       if (runType[iFED] == EcalDCCHeaderBlock::COSMIC || runType[iFED] == EcalDCCHeaderBlock::MTCC ||
0047           runType[iFED] == EcalDCCHeaderBlock::COSMICS_GLOBAL || runType[iFED] == EcalDCCHeaderBlock::PHYSICS_GLOBAL ||
0048           runType[iFED] == EcalDCCHeaderBlock::COSMICS_LOCAL || runType[iFED] == EcalDCCHeaderBlock::PHYSICS_LOCAL)
0049         return true;
0050     }
0051 
0052     return false;
0053   }
0054 
0055   void PiZeroTask::runOnEBRecHits(EcalRecHitCollection const& hits) {
0056     MESet& mePi0MinvEB(MEs_.at("Pi0MinvEB"));
0057     MESet& mePi0Pt1EB(MEs_.at("Pi0Pt1EB"));
0058     MESet& mePi0Pt2EB(MEs_.at("Pi0Pt2EB"));
0059     MESet& mePi0PtEB(MEs_.at("Pi0PtEB"));
0060     MESet& mePi0IsoEB(MEs_.at("Pi0IsoEB"));
0061 
0062     const CaloSubdetectorTopology* topology_p;
0063     const CaloSubdetectorGeometry* geometry_p = GetGeometry()->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0064     const CaloSubdetectorGeometry* geometryES_p = GetGeometry()->getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
0065 
0066     // Parameters for the position calculation:
0067     PositionCalc posCalculator_ = PositionCalc(posCalcParameters_);
0068 
0069     std::map<DetId, EcalRecHit> recHitsEB_map;
0070 
0071     std::vector<EcalRecHit> seeds;
0072     std::vector<EBDetId> usedXtals;
0073     seeds.clear();
0074     usedXtals.clear();
0075 
0076     int nClus = 0;
0077     std::vector<float> eClus;
0078     std::vector<float> etClus;
0079     std::vector<float> etaClus;
0080     std::vector<float> phiClus;
0081     std::vector<EBDetId> max_hit;
0082     std::vector<std::vector<EcalRecHit> > RecHitsCluster;
0083     std::vector<float> s4s9Clus;
0084 
0085     // Find cluster seeds in EB
0086     for (auto const& hit : hits) {
0087       EBDetId id(hit.id());
0088       double energy = hit.energy();
0089       if (energy > seleXtalMinEnergy_) {
0090         std::pair<DetId, EcalRecHit> map_entry(hit.id(), hit);
0091         recHitsEB_map.insert(map_entry);
0092       }
0093       if (energy > clusSeedThr_)
0094         seeds.push_back(hit);
0095     }  // EB rechits
0096 
0097     sort(seeds.begin(), seeds.end(), [](auto& x, auto& y) { return (x.energy() > y.energy()); });
0098     for (auto const& seed : seeds) {
0099       EBDetId seed_id = seed.id();
0100 
0101       bool seedAlreadyUsed = false;
0102       for (auto const& usedIds : usedXtals) {
0103         if (usedIds == seed_id) {
0104           seedAlreadyUsed = true;
0105           break;
0106         }
0107       }
0108       if (seedAlreadyUsed)
0109         continue;
0110       topology_p = GetTopology()->getSubdetectorTopology(DetId::Ecal, EcalBarrel);
0111       std::vector<DetId> clus_v = topology_p->getWindow(seed_id, clusEtaSize_, clusPhiSize_);
0112       std::vector<std::pair<DetId, float> > clus_used;
0113 
0114       std::vector<EcalRecHit> RecHitsInWindow;
0115 
0116       double simple_energy = 0;
0117 
0118       for (auto const& det : clus_v) {
0119         bool HitAlreadyUsed = false;
0120         for (auto const& usedIds : usedXtals) {
0121           if (usedIds == det) {
0122             HitAlreadyUsed = true;
0123             break;
0124           }
0125         }
0126         if (HitAlreadyUsed)
0127           continue;
0128         if (recHitsEB_map.find(det) != recHitsEB_map.end()) {
0129           std::map<DetId, EcalRecHit>::iterator aHit;
0130           aHit = recHitsEB_map.find(det);
0131           usedXtals.push_back(det);
0132           RecHitsInWindow.push_back(aHit->second);
0133           clus_used.push_back(std::pair<DetId, float>(det, 1.));
0134           simple_energy = simple_energy + aHit->second.energy();
0135         }
0136       }
0137 
0138       math::XYZPoint clus_pos = posCalculator_.Calculate_Location(clus_used, &hits, geometry_p, geometryES_p);
0139       float theta_s = 2. * atan(exp(-clus_pos.eta()));
0140       float p0x_s = simple_energy * sin(theta_s) * cos(clus_pos.phi());
0141       float p0y_s = simple_energy * sin(theta_s) * sin(clus_pos.phi());
0142       float et_s = sqrt(p0x_s * p0x_s + p0y_s * p0y_s);
0143 
0144       eClus.push_back(simple_energy);
0145       etClus.push_back(et_s);
0146       etaClus.push_back(clus_pos.eta());
0147       phiClus.push_back(clus_pos.phi());
0148       max_hit.push_back(seed_id);
0149       RecHitsCluster.push_back(RecHitsInWindow);
0150 
0151       // Compute S4/S9 variable
0152       // We are not sure to have 9 RecHits so need to check eta and phi:
0153       float s4s9_[4];
0154       for (int i = 0; i < 4; i++)
0155         s4s9_[i] = seed.energy();
0156       for (unsigned int j = 0; j < RecHitsInWindow.size(); j++) {
0157         if ((((EBDetId)RecHitsInWindow[j].id()).ieta() == seed_id.ieta() - 1 && seed_id.ieta() != 1) ||
0158             (seed_id.ieta() == 1 && (((EBDetId)RecHitsInWindow[j].id()).ieta() == seed_id.ieta() - 2))) {
0159           if (((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi() - 1 ||
0160               ((EBDetId)RecHitsInWindow[j].id()).iphi() - 360 == seed_id.iphi() - 1) {
0161             s4s9_[0] += RecHitsInWindow[j].energy();
0162           } else {
0163             if (((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi()) {
0164               s4s9_[0] += RecHitsInWindow[j].energy();
0165               s4s9_[1] += RecHitsInWindow[j].energy();
0166             } else {
0167               if (((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi() + 1 ||
0168                   ((EBDetId)RecHitsInWindow[j].id()).iphi() - 360 == seed_id.iphi() + 1) {
0169                 s4s9_[1] += RecHitsInWindow[j].energy();
0170               }
0171             }
0172           }
0173         } else {
0174           if (((EBDetId)RecHitsInWindow[j].id()).ieta() == seed_id.ieta()) {
0175             if (((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi() - 1 ||
0176                 ((EBDetId)RecHitsInWindow[j].id()).iphi() - 360 == seed_id.iphi() - 1) {
0177               s4s9_[0] += RecHitsInWindow[j].energy();
0178               s4s9_[3] += RecHitsInWindow[j].energy();
0179             } else {
0180               if (((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi() + 1 ||
0181                   ((EBDetId)RecHitsInWindow[j].id()).iphi() - 360 == seed_id.iphi() + 1) {
0182                 s4s9_[1] += RecHitsInWindow[j].energy();
0183                 s4s9_[2] += RecHitsInWindow[j].energy();
0184               }
0185             }
0186           } else {
0187             if ((((EBDetId)RecHitsInWindow[j].id()).ieta() == seed_id.ieta() + 1 && seed_id.ieta() != -1) ||
0188                 (seed_id.ieta() == -1 && (((EBDetId)RecHitsInWindow[j].id()).ieta() == seed_id.ieta() + 2))) {
0189               if (((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi() - 1 ||
0190                   ((EBDetId)RecHitsInWindow[j].id()).iphi() - 360 == seed_id.iphi() - 1) {
0191                 s4s9_[3] += RecHitsInWindow[j].energy();
0192               } else {
0193                 if (((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi()) {
0194                   s4s9_[2] += RecHitsInWindow[j].energy();
0195                   s4s9_[3] += RecHitsInWindow[j].energy();
0196                 } else {
0197                   if (((EBDetId)RecHitsInWindow[j].id()).iphi() == seed_id.iphi() + 1 ||
0198                       ((EBDetId)RecHitsInWindow[j].id()).iphi() - 360 == seed_id.iphi() + 1) {
0199                     s4s9_[2] += RecHitsInWindow[j].energy();
0200                   }
0201                 }
0202               }
0203             } else {
0204               edm::LogWarning("EcalDQM") << " (EBDetId)RecHitsInWindow[j].id()).ieta() "
0205                                          << ((EBDetId)RecHitsInWindow[j].id()).ieta() << " seed_id.ieta() "
0206                                          << seed_id.ieta() << "\n"
0207                                          << " Problem with S4 calculation\n";
0208               return;
0209             }
0210           }
0211         }
0212       }
0213       s4s9Clus.push_back(*std::max_element(s4s9_, s4s9_ + 4) / simple_energy);
0214       nClus++;
0215       if (nClus == MAXCLUS)
0216         return;
0217     }  //  End loop over seed clusters
0218 
0219     // Selection, based on simple clustering
0220     // pi0 candidates
0221     int npi0_s = 0;
0222 
0223     std::vector<EBDetId> scXtals;
0224     scXtals.clear();
0225 
0226     if (nClus <= 1)
0227       return;
0228     for (Int_t i = 0; i < nClus; i++) {
0229       for (Int_t j = i + 1; j < nClus; j++) {
0230         if (etClus[i] > selePtGammaOne_ && etClus[j] > selePtGammaTwo_ && s4s9Clus[i] > seleS4S9GammaOne_ &&
0231             s4s9Clus[j] > seleS4S9GammaTwo_) {
0232           float theta_0 = 2. * atan(exp(-etaClus[i]));
0233           float theta_1 = 2. * atan(exp(-etaClus[j]));
0234 
0235           float p0x = eClus[i] * sin(theta_0) * cos(phiClus[i]);
0236           float p1x = eClus[j] * sin(theta_1) * cos(phiClus[j]);
0237           float p0y = eClus[i] * sin(theta_0) * sin(phiClus[i]);
0238           float p1y = eClus[j] * sin(theta_1) * sin(phiClus[j]);
0239           float p0z = eClus[i] * cos(theta_0);
0240           float p1z = eClus[j] * cos(theta_1);
0241 
0242           float pt_pi0 = sqrt((p0x + p1x) * (p0x + p1x) + (p0y + p1y) * (p0y + p1y));
0243           if (pt_pi0 < selePtPi0_)
0244             continue;
0245           float m_inv = sqrt((eClus[i] + eClus[j]) * (eClus[i] + eClus[j]) - (p0x + p1x) * (p0x + p1x) -
0246                              (p0y + p1y) * (p0y + p1y) - (p0z + p1z) * (p0z + p1z));
0247           if ((m_inv < seleMinvMaxPi0_) && (m_inv > seleMinvMinPi0_)) {
0248             // New Loop on cluster to measure isolation:
0249             std::vector<int> IsoClus;
0250             IsoClus.clear();
0251             float Iso = 0;
0252             TVector3 pi0vect = TVector3((p0x + p1x), (p0y + p1y), (p0z + p1z));
0253             for (Int_t k = 0; k < nClus; k++) {
0254               if (k == i || k == j)
0255                 continue;
0256               TVector3 Clusvect = TVector3(eClus[k] * sin(2. * atan(exp(-etaClus[k]))) * cos(phiClus[k]),
0257                                            eClus[k] * sin(2. * atan(exp(-etaClus[k]))) * sin(phiClus[k]),
0258                                            eClus[k] * cos(2. * atan(exp(-etaClus[k]))));
0259               float dretaclpi0 = fabs(etaClus[k] - pi0vect.Eta());
0260               float drclpi0 = Clusvect.DeltaR(pi0vect);
0261 
0262               if ((drclpi0 < selePi0BeltDR_) && (dretaclpi0 < selePi0BeltDeta_)) {
0263                 Iso = Iso + etClus[k];
0264                 IsoClus.push_back(k);
0265               }
0266             }
0267 
0268             if (Iso / pt_pi0 < selePi0Iso_) {
0269               mePi0MinvEB.fill(getEcalDQMSetupObjects(), m_inv);
0270               mePi0Pt1EB.fill(getEcalDQMSetupObjects(), etClus[i]);
0271               mePi0Pt2EB.fill(getEcalDQMSetupObjects(), etClus[j]);
0272               mePi0PtEB.fill(getEcalDQMSetupObjects(), pt_pi0);
0273               mePi0IsoEB.fill(getEcalDQMSetupObjects(), Iso / pt_pi0);
0274 
0275               npi0_s++;
0276             }
0277 
0278             if (npi0_s == MAXPI0S)
0279               return;
0280           }  // pi0 inv mass window
0281         }  // pt and S4S9 cut
0282       }  // cluster "j" index loop
0283     }  // cluster "i" index loop
0284   }  // runonEBRecHits()
0285 
0286   DEFINE_ECALDQM_WORKER(PiZeroTask);
0287 }  // namespace ecaldqm
[ 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!