EcalCalibAlgos/src/ZeeKinematicTools.cc

0001
0002 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0003
0004 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0005 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0006 #include "DataFormats/DetId/interface/DetId.h"
0007 #include "CondFormats/EcalObjects/interface/EcalIntercalibConstants.h"
0008 #include "CondFormats/DataRecord/interface/EcalIntercalibConstantsRcd.h"
0009 #include "Calibration/Tools/interface/calibXMLwriter.h"
0010 #include "Calibration/Tools/interface/CalibrationCluster.h"
0011 #include "Calibration/Tools/interface/HouseholderDecomposition.h"
0012 #include "Calibration/Tools/interface/MinL3Algorithm.h"
0013 #include "Calibration/EcalCalibAlgos/interface/ZeeKinematicTools.h"
0014 #include "Calibration/Tools/interface/ZIterativeAlgorithmWithFit.h"
0015 #include "DataFormats/EgammaCandidates/interface/Electron.h"
0016 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0017 #include "DataFormats/TrackReco/interface/Track.h"
0018 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0019 #include "DataFormats/TrackReco/interface/TrackExtraFwd.h"
0020 #include "TFile.h"
0021 #include "TH1.h"
0022 #include "TH2.h"
0023 #include "TF1.h"
0024 #include "TRandom.h"
0025
0026 #include <iostream>
0027 #include <string>
0028 #include <stdexcept>
0029 #include <vector>
0030
0031 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0032
0033 ZeeKinematicTools::ZeeKinematicTools() {}
0034
0035 ZeeKinematicTools::~ZeeKinematicTools() {}
0036
0037 //--------------------------------------------
0038
0039 float ZeeKinematicTools::cosThetaElectrons_SC(const std::pair<calib::CalibElectron*, calib::CalibElectron*>& aZCandidate,
0040                                               float ele1EnergyCorrection,
0041                                               float ele2EnergyCorrection) {
0042   float theta1 = 2. * atan(exp(-aZCandidate.first->getRecoElectron()->superCluster()->eta()));
0043   float phi1 = aZCandidate.first->getRecoElectron()->superCluster()->phi();
0044
0045   float x1 = aZCandidate.first->getRecoElectron()->superCluster()->energy() * sin(theta1) * cos(phi1);
0046   float y1 = aZCandidate.first->getRecoElectron()->superCluster()->energy() * sin(theta1) * sin(phi1);
0047   float z1 = aZCandidate.first->getRecoElectron()->superCluster()->energy() * cos(theta1);
0048   float mod1 = sqrt(x1 * x1 + y1 * y1 + z1 * z1);
0049
0050   float theta2 = 2. * atan(exp(-aZCandidate.second->getRecoElectron()->superCluster()->eta()));
0051   float phi2 = aZCandidate.second->getRecoElectron()->superCluster()->phi();
0052
0053   float x2 = aZCandidate.second->getRecoElectron()->superCluster()->energy() * sin(theta2) * cos(phi2);
0054   float y2 = aZCandidate.second->getRecoElectron()->superCluster()->energy() * sin(theta2) * sin(phi2);
0055   float z2 = aZCandidate.second->getRecoElectron()->superCluster()->energy() * cos(theta2);
0056   float mod2 = sqrt(x2 * x2 + y2 * y2 + z2 * z2);
0057
0058   return (x1 * x2 + y1 * y2 + z1 * z2) / (mod1 * mod2);
0059 }
0060
0061 //--------------------------------------------
0062
0063 float ZeeKinematicTools::cosThetaElectrons_TK(const std::pair<calib::CalibElectron*, calib::CalibElectron*>& aZCandidate,
0064                                               float ele1EnergyCorrection,
0065                                               float ele2EnergyCorrection) {
0066   float theta1 = 2. * atan(exp(-aZCandidate.first->getRecoElectron()->eta()));
0067   float phi1 = aZCandidate.first->getRecoElectron()->phi();
0068
0069   float x1 = aZCandidate.first->getRecoElectron()->superCluster()->energy() * sin(theta1) * cos(phi1);
0070   float y1 = aZCandidate.first->getRecoElectron()->superCluster()->energy() * sin(theta1) * sin(phi1);
0071   float z1 = aZCandidate.first->getRecoElectron()->superCluster()->energy() * cos(theta1);
0072   float mod1 = sqrt(x1 * x1 + y1 * y1 + z1 * z1);
0073
0074   float theta2 = 2. * atan(exp(-aZCandidate.second->getRecoElectron()->eta()));
0075   float phi2 = aZCandidate.second->getRecoElectron()->phi();
0076
0077   float x2 = aZCandidate.second->getRecoElectron()->superCluster()->energy() * sin(theta2) * cos(phi2);
0078   float y2 = aZCandidate.second->getRecoElectron()->superCluster()->energy() * sin(theta2) * sin(phi2);
0079   float z2 = aZCandidate.second->getRecoElectron()->superCluster()->energy() * cos(theta2);
0080   float mod2 = sqrt(x2 * x2 + y2 * y2 + z2 * z2);
0081
0082   return (x1 * x2 + y1 * y2 + z1 * z2) / (mod1 * mod2);
0083 }
0084
0085 //--------------------------------------------
0086
0087 float ZeeKinematicTools::calculateZMassWithCorrectedElectrons_noTK(
0088     const std::pair<calib::CalibElectron*, calib::CalibElectron*>& aZCandidate,
0089     float ele1EnergyCorrection,
0090     float ele2EnergyCorrection) {
0091   return ZIterativeAlgorithmWithFit::invMassCalc(
0092       aZCandidate.first->getRecoElectron()->superCluster()->energy() / ele1EnergyCorrection,
0093       aZCandidate.first->getRecoElectron()->superCluster()->eta(),
0094       aZCandidate.first->getRecoElectron()->superCluster()->phi(),
0095       aZCandidate.second->getRecoElectron()->superCluster()->energy() / ele2EnergyCorrection,
0096       aZCandidate.second->getRecoElectron()->superCluster()->eta(),
0097       aZCandidate.second->getRecoElectron()->superCluster()->phi());
0098 }
0099
0100 //--------------------------------------------
0101
0102 float ZeeKinematicTools::calculateZMassWithCorrectedElectrons_withTK(
0103     const std::pair<calib::CalibElectron*, calib::CalibElectron*>& aZCandidate,
0104     float ele1EnergyCorrection,
0105     float ele2EnergyCorrection) {
0106   return ZIterativeAlgorithmWithFit::invMassCalc(
0107       aZCandidate.first->getRecoElectron()->superCluster()->energy() / ele1EnergyCorrection,
0108       aZCandidate.first->getRecoElectron()->eta(),
0109       aZCandidate.first->getRecoElectron()->phi(),
0110       aZCandidate.second->getRecoElectron()->superCluster()->energy() / ele2EnergyCorrection,
0111       aZCandidate.second->getRecoElectron()->eta(),
0112       aZCandidate.second->getRecoElectron()->phi());
0113 }
0114
0115 //--------------------------------------------
0116
0117 float ZeeKinematicTools::calculateZMass_noTK(
0118     const std::pair<calib::CalibElectron*, calib::CalibElectron*>& aZCandidate) {
0119   return ZIterativeAlgorithmWithFit::invMassCalc(aZCandidate.first->getRecoElectron()->superCluster()->energy(),
0120                                                  aZCandidate.first->getRecoElectron()->superCluster()->eta(),
0121                                                  aZCandidate.first->getRecoElectron()->superCluster()->phi(),
0122                                                  aZCandidate.second->getRecoElectron()->superCluster()->energy(),
0123                                                  aZCandidate.second->getRecoElectron()->superCluster()->eta(),
0124                                                  aZCandidate.second->getRecoElectron()->superCluster()->phi());
0125 }
0126
0127 //--------------------------------------------
0128
0129 float ZeeKinematicTools::calculateZMass_withTK(
0130     const std::pair<calib::CalibElectron*, calib::CalibElectron*>& aZCandidate) {
0131   return ZIterativeAlgorithmWithFit::invMassCalc(aZCandidate.first->getRecoElectron()->superCluster()->energy(),
0132                                                  aZCandidate.first->getRecoElectron()->eta(),
0133                                                  aZCandidate.first->getRecoElectron()->phi(),
0134                                                  aZCandidate.second->getRecoElectron()->superCluster()->energy(),
0135                                                  aZCandidate.second->getRecoElectron()->eta(),
0136                                                  aZCandidate.second->getRecoElectron()->phi());
0137 }
0138
0139 //--------------------------------------------
0140
0141 float ZeeKinematicTools::calculateZRapidity(const std::pair<calib::CalibElectron*, calib::CalibElectron*>& aZCandidate) {
0142   TLorentzVector ele1LV(aZCandidate.first->getRecoElectron()->px(),
0143                         aZCandidate.first->getRecoElectron()->py(),
0144                         aZCandidate.first->getRecoElectron()->pz(),
0145                         aZCandidate.first->getRecoElectron()->superCluster()->energy());
0146
0147   TLorentzVector ele2LV(aZCandidate.second->getRecoElectron()->px(),
0148                         aZCandidate.second->getRecoElectron()->py(),
0149                         aZCandidate.second->getRecoElectron()->pz(),
0150                         aZCandidate.second->getRecoElectron()->superCluster()->energy());
0151
0152   return (ele1LV + ele2LV).Rapidity();
0153 }
0154
0155 //--------------------------------------------
0156
0157 float ZeeKinematicTools::calculateZEta(const std::pair<calib::CalibElectron*, calib::CalibElectron*>& aZCandidate) {
0158   TLorentzVector ele1LV(aZCandidate.first->getRecoElectron()->px(),
0159                         aZCandidate.first->getRecoElectron()->py(),
0160                         aZCandidate.first->getRecoElectron()->pz(),
0161                         aZCandidate.first->getRecoElectron()->superCluster()->energy());
0162
0163   TLorentzVector ele2LV(aZCandidate.second->getRecoElectron()->px(),
0164                         aZCandidate.second->getRecoElectron()->py(),
0165                         aZCandidate.second->getRecoElectron()->pz(),
0166                         aZCandidate.second->getRecoElectron()->superCluster()->energy());
0167
0168   return (ele1LV + ele2LV).Eta();
0169 }
0170
0171 //--------------------------------------------
0172
0173 float ZeeKinematicTools::calculateZTheta(const std::pair<calib::CalibElectron*, calib::CalibElectron*>& aZCandidate) {
0174   TLorentzVector ele1LV(aZCandidate.first->getRecoElectron()->px(),
0175                         aZCandidate.first->getRecoElectron()->py(),
0176                         aZCandidate.first->getRecoElectron()->pz(),
0177                         aZCandidate.first->getRecoElectron()->superCluster()->energy());
0178
0179   TLorentzVector ele2LV(aZCandidate.second->getRecoElectron()->px(),
0180                         aZCandidate.second->getRecoElectron()->py(),
0181                         aZCandidate.second->getRecoElectron()->pz(),
0182                         aZCandidate.second->getRecoElectron()->superCluster()->energy());
0183
0184   return (ele1LV + ele2LV).Theta();
0185 }
0186
0187 //--------------------------------------------
0188
0189 float ZeeKinematicTools::calculateZPhi(const std::pair<calib::CalibElectron*, calib::CalibElectron*>& aZCandidate) {
0190   TLorentzVector ele1LV(aZCandidate.first->getRecoElectron()->px(),
0191                         aZCandidate.first->getRecoElectron()->py(),
0192                         aZCandidate.first->getRecoElectron()->pz(),
0193                         aZCandidate.first->getRecoElectron()->superCluster()->energy());
0194
0195   TLorentzVector ele2LV(aZCandidate.second->getRecoElectron()->px(),
0196                         aZCandidate.second->getRecoElectron()->py(),
0197                         aZCandidate.second->getRecoElectron()->pz(),
0198                         aZCandidate.second->getRecoElectron()->superCluster()->energy());
0199
0200   return (ele1LV + ele2LV).Phi();
0201 }
0202
0203 //--------------------------------------------
0204
0205 float ZeeKinematicTools::calculateZPt(const std::pair<calib::CalibElectron*, calib::CalibElectron*>& aZCandidate) {
0206   TLorentzVector ele1LV(aZCandidate.first->getRecoElectron()->px(),
0207                         aZCandidate.first->getRecoElectron()->py(),
0208                         aZCandidate.first->getRecoElectron()->pz(),
0209                         aZCandidate.first->getRecoElectron()->superCluster()->energy());
0210
0211   TLorentzVector ele2LV(aZCandidate.second->getRecoElectron()->px(),
0212                         aZCandidate.second->getRecoElectron()->py(),
0213                         aZCandidate.second->getRecoElectron()->pz(),
0214                         aZCandidate.second->getRecoElectron()->superCluster()->energy());
0215
0216   return (ele1LV + ele2LV).Pt();
0217 }