Type1MET/data/FWLiteFunc.C

0001 #if !defined(__CINT__) && !defined(__MAKECINT__)
0002 #include "DataFormats/FWLite/interface/Handle.h"
0003 #include "DataFormats/FWLite/interface/Event.h"
0004 #include "DataFormats/METReco/interface/CaloMET.h"
0005 #include "DataFormats/METReco/interface/METCollection.h"
0006 #include "DataFormats/METReco/interface/MET.h"
0007 #include "DataFormats/METReco/interface/CaloMETCollection.h"
0008 #include "DataFormats/METReco/interface/SpecificCaloMETData.h"
0009 #include "DataFormats/JetReco/interface/CaloJet.h"
0010 #include "DataFormats/TrackReco/interface/Track.h"
0011 #include "DataFormats/MuonReco/interface/Muon.h"
0012 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0013 #endif
0014
0015 namespace reco {
0016   class Muon;
0017   class CaloMET;
0018 };
0019 #include "DataFormats/Math/interface/Point3D.h"
0020 #include "DataFormats/Math/interface/LorentzVector.h"
0021 #include "RecoMET/METAlgorithms/interface/MuonMETInfo.h"
0022
0023 #include <vector>
0024 #include "TMath.h"
0025 #include "TFile.h"
0026
0027
0028 //----------------------------------------------------------------------------
0029 void correctMETforMuon(double& metx, double& mety, double bfield, int muonCharge,
0030                math::XYZTLorentzVector muonP4, math::XYZPoint muonVertex,
0031                MuonMETInfo& muonMETInfo) {
0032
0033
0034
0035   double mu_p     = muonP4.P();
0036   double mu_pt    = muonP4.Pt();
0037   double mu_phi   = muonP4.Phi();
0038   double mu_eta   = muonP4.Eta();
0039   double mu_vz    = muonVertex.z()/100.;
0040   double mu_pz    = muonP4.Pz();
0041
0042
0043   //add in the muon's pt
0044   metx -= mu_pt*cos(mu_phi);
0045   mety -= mu_pt*sin(mu_phi);
0046
0047   double ecalPhi, ecalTheta;
0048   double hcalPhi, hcalTheta;
0049   double hoPhi, hoTheta;
0050
0051
0052   //should always be false for FWLite
0053   //unless you want to supply co-ordinates at
0054   //the calorimeter sub-detectors yourself
0055   if(muonMETInfo.useTkAssociatorPositions) {
0056     ecalPhi   = muonMETInfo.ecalPos.Phi();
0057     ecalTheta = muonMETInfo.ecalPos.Theta();
0058     hcalPhi   = muonMETInfo.hcalPos.Phi();
0059     hcalTheta = muonMETInfo.hcalPos.Theta();
0060     hoPhi     = muonMETInfo.hoPos.Phi();
0061     hoTheta   = muonMETInfo.hoPos.Theta();
0062   } else {
0063
0064     /*
0065       use the analytical solution for the
0066       intersection of a helix with a cylinder
0067       to find the positions of the muon
0068       at the various calo surfaces
0069     */
0070     double rEcal = 1.290;
0071     double rHcal = 1.9;
0072     double rHo   = 3.82;
0073     if(abs(mu_eta) > 0.3) rHo = 4.07;
0074     //distance from the center of detector to face of Ecal
0075     double zFaceEcal = 3.209;
0076     if(mu_eta < 0 ) zFaceEcal = -1*zFaceEcal;
0077     //distance from the center of detector to face of Hcal
0078     double zFaceHcal = 3.88;
0079     if(mu_eta < 0 ) zFaceHcal = -1*zFaceHcal;
0080
0081     //now we have to get Phi
0082     //bending radius of the muon (units are meters)
0083     double bendr = mu_pt*1000/(300*bfield);
0084
0085     double tb_ecal = TMath::ACos(1-rEcal*rEcal/(2*bendr*bendr)); //helix time interval parameter
0086     double tb_hcal = TMath::ACos(1-rHcal*rHcal/(2*bendr*bendr)); //helix time interval parameter
0087     double tb_ho   = TMath::ACos(1-rHo*rHo/(2*bendr*bendr));     //helix time interval parameter
0088     double xEcal,yEcal,zEcal;
0089     double xHcal,yHcal,zHcal;
0090     double xHo, yHo,zHo;
0091     //Ecal
0092     //in the barrel and if not a looper
0093     if(fabs(mu_pz*bendr*tb_ecal/mu_pt+mu_vz) < fabs(zFaceEcal) && rEcal < 2*bendr) {
0094       xEcal = bendr*(TMath::Sin(tb_ecal+mu_phi)-TMath::Sin(mu_phi));
0095       yEcal = bendr*(-TMath::Cos(tb_ecal+mu_phi)+TMath::Cos(mu_phi));
0096       zEcal = bendr*tb_ecal*mu_pz/mu_pt + mu_vz;
0097     } else { //endcap
0098       if(mu_pz > 0) {
0099         double te_ecal = (fabs(zFaceEcal) - mu_vz)*mu_pt/(bendr*mu_pz);
0100         xEcal = bendr*(TMath::Sin(te_ecal+mu_phi) - TMath::Sin(mu_phi));
0101         yEcal = bendr*(-TMath::Cos(te_ecal+mu_phi) + TMath::Cos(mu_phi));
0102         zEcal = fabs(zFaceEcal);
0103       } else {
0104         double te_ecal = -(fabs(zFaceEcal) + mu_vz)*mu_pt/(bendr*mu_pz);
0105         xEcal = bendr*(TMath::Sin(te_ecal+mu_phi) - TMath::Sin(mu_phi));
0106         yEcal = bendr*(-TMath::Cos(te_ecal+mu_phi) + TMath::Cos(mu_phi));
0107         zEcal = -fabs(zFaceEcal);
0108       }
0109     }
0110
0111     //Hcal
0112     if(fabs(mu_pz*bendr*tb_hcal/mu_pt+mu_vz) < fabs(zFaceHcal) && rEcal < 2*bendr) { //in the barrel
0113       xHcal = bendr*(TMath::Sin(tb_hcal+mu_phi)-TMath::Sin(mu_phi));
0114       yHcal = bendr*(-TMath::Cos(tb_hcal+mu_phi)+TMath::Cos(mu_phi));
0115       zHcal = bendr*tb_hcal*mu_pz/mu_pt + mu_vz;
0116     } else { //endcap
0117       if(mu_pz > 0) {
0118         double te_hcal = (fabs(zFaceHcal) - mu_vz)*mu_pt/(bendr*mu_pz);
0119         xHcal = bendr*(TMath::Sin(te_hcal+mu_phi) - TMath::Sin(mu_phi));
0120         yHcal = bendr*(-TMath::Cos(te_hcal+mu_phi) + TMath::Cos(mu_phi));
0121         zHcal = fabs(zFaceHcal);
0122       } else {
0123         double te_hcal = -(fabs(zFaceHcal) + mu_vz)*mu_pt/(bendr*mu_pz);
0124         xHcal = bendr*(TMath::Sin(te_hcal+mu_phi) - TMath::Sin(mu_phi));
0125         yHcal = bendr*(-TMath::Cos(te_hcal+mu_phi) + TMath::Cos(mu_phi));
0126         zHcal = -fabs(zFaceHcal);
0127       }
0128     }
0129
0130     //Ho - just the barrel
0131     xHo = bendr*(TMath::Sin(tb_ho+mu_phi)-TMath::Sin(mu_phi));
0132     yHo = bendr*(-TMath::Cos(tb_ho+mu_phi)+TMath::Cos(mu_phi));
0133     zHo = bendr*tb_ho*mu_pz/mu_pt + mu_vz;
0134
0135     ecalTheta = TMath::ACos(zEcal/sqrt(pow(xEcal,2) + pow(yEcal,2)+pow(zEcal,2)));
0136     ecalPhi   = atan2(yEcal,xEcal);
0137     hcalTheta = TMath::ACos(zHcal/sqrt(pow(xHcal,2) + pow(yHcal,2)+pow(zHcal,2)));
0138     hcalPhi   = atan2(yHcal,xHcal);
0139     hoTheta   = TMath::ACos(zHo/sqrt(pow(xHo,2) + pow(yHo,2)+pow(zHo,2)));
0140     hoPhi     = atan2(yHo,xHo);
0141
0142     //2d radius in x-y plane
0143     double r2dEcal = sqrt(pow(xEcal,2)+pow(yEcal,2));
0144     double r2dHcal = sqrt(pow(xHcal,2)+pow(yHcal,2));
0145     double r2dHo   = sqrt(pow(xHo,2)  +pow(yHo,2));
0146
0147     /*
0148       the above prescription is for right handed helicies only
0149       Positively charged muons trace a left handed helix
0150       so we correct for that
0151     */
0152     if(muonCharge > 0) {
0153
0154       //Ecal
0155       double dphi = mu_phi - ecalPhi;
0156       if(fabs(dphi) > TMath::Pi())
0157         dphi = 2*TMath::Pi() - fabs(dphi);
0158       ecalPhi = mu_phi - fabs(dphi);
0159       if(fabs(ecalPhi) > TMath::Pi()) {
0160         double temp = 2*TMath::Pi() - fabs(ecalPhi);
0161         ecalPhi = -1*temp*ecalPhi/fabs(ecalPhi);
0162       }
0163       xEcal = r2dEcal*TMath::Cos(ecalPhi);
0164       yEcal = r2dEcal*TMath::Sin(ecalPhi);
0165
0166       //Hcal
0167       dphi = mu_phi - hcalPhi;
0168       if(fabs(dphi) > TMath::Pi())
0169         dphi = 2*TMath::Pi() - fabs(dphi);
0170       hcalPhi = mu_phi - fabs(dphi);
0171       if(fabs(hcalPhi) > TMath::Pi()) {
0172         double temp = 2*TMath::Pi() - fabs(hcalPhi);
0173         hcalPhi = -1*temp*hcalPhi/fabs(hcalPhi);
0174       }
0175       xHcal = r2dHcal*TMath::Cos(hcalPhi);
0176       yHcal = r2dHcal*TMath::Sin(hcalPhi);
0177
0178       //Ho
0179       dphi = mu_phi - hoPhi;
0180       if(fabs(dphi) > TMath::Pi())
0181         dphi = 2*TMath::Pi() - fabs(dphi);
0182       hoPhi = mu_phi - fabs(dphi);
0183       if(fabs(hoPhi) > TMath::Pi()) {
0184         double temp = 2*TMath::Pi() - fabs(hoPhi);
0185         hoPhi = -1*temp*hoPhi/fabs(hoPhi);
0186       }
0187       xHo = r2dHo*TMath::Cos(hoPhi);
0188       yHo = r2dHo*TMath::Sin(hoPhi);
0189
0190     }
0191   }
0192
0193   //being safe...
0194   if(!muonMETInfo.useHO) muonMETInfo.hoE = 0;
0195
0196   //for isolated muons
0197   if(!muonMETInfo.useAverage) {
0198
0199     double mu_Ex =  muonMETInfo.ecalE*sin(ecalTheta)*cos(ecalPhi)
0200       + muonMETInfo.hcalE*sin(hcalTheta)*cos(hcalPhi)
0201       + muonMETInfo.hoE*sin(hoTheta)*cos(hoPhi);
0202     double mu_Ey =  muonMETInfo.ecalE*sin(ecalTheta)*sin(ecalPhi)
0203       + muonMETInfo.hcalE*sin(hcalTheta)*sin(hcalPhi)
0204       + muonMETInfo.hoE*sin(hoTheta)*sin(hoPhi);
0205
0206     metx += mu_Ex;
0207     mety += mu_Ey;
0208
0209   } else { //non-isolated muons - derive the correction
0210
0211     //dE/dx in matter for iron:
0212     //-(11.4 + 0.96*fabs(log(p0*2.8)) + 0.033*p0*(1.0 - pow(p0, -0.33)) )*1e-3
0213     //from http://cmslxr.fnal.gov/lxr/source/TrackPropagation/SteppingHelixPropagator/src/SteppingHelixPropagator.ccyes,
0214     //line ~1100
0215     //normalisation is at 50 GeV
0216     double dEdx_normalization = -(11.4 + 0.96*fabs(log(50*2.8)) + 0.033*50*(1.0 - pow(50, -0.33)) )*1e-3;
0217     double dEdx_numerator     = -(11.4 + 0.96*fabs(log(mu_p*2.8)) + 0.033*mu_p*(1.0 - pow(mu_p, -0.33)) )*1e-3;
0218
0219     double temp = 0.0;
0220
0221     if(muonMETInfo.useHO) {
0222       //for the Towers, with HO
0223       if(fabs(mu_eta) < 0.2)
0224         temp = 2.75*(1-0.00003*mu_p);
0225       if(fabs(mu_eta) > 0.2 && fabs(mu_eta) < 1.0)
0226         temp = (2.38+0.0144*fabs(mu_eta))*(1-0.0003*mu_p);
0227       if(fabs(mu_eta) > 1.0 && fabs(mu_eta) < 1.3)
0228         temp = 7.413-5.12*fabs(mu_eta);
0229       if(fabs(mu_eta) > 1.3)
0230         temp = 2.084-0.743*fabs(mu_eta);
0231     } else {
0232       if(fabs(mu_eta) < 1.0)
0233         temp = 2.33*(1-0.0004*mu_p);
0234       if(fabs(mu_eta) > 1.0 && fabs(mu_eta) < 1.3)
0235         temp = (7.413-5.12*fabs(mu_eta))*(1-0.0003*mu_p);
0236       if(fabs(mu_eta) > 1.3)
0237         temp = 2.084-0.743*fabs(mu_eta);
0238     }
0239
0240     double dep = temp*dEdx_normalization/dEdx_numerator;
0241     if(dep < 0.5) dep = 0.0;
0242     //use the average phi of the 3 subdetectors
0243     if(fabs(mu_eta) < 1.3) {
0244       metx += dep*cos((ecalPhi+hcalPhi+hoPhi)/3);
0245       mety += dep*sin((ecalPhi+hcalPhi+hoPhi)/3);
0246     } else {
0247       metx += dep*cos( (ecalPhi+hcalPhi)/2);
0248       mety += dep*cos( (ecalPhi+hcalPhi)/2);
0249     }
0250   }
0251
0252
0253 }
0254 //----------------------------------------------------------------------------
0255
0256
0257 void correctMETforMuon(double& metx, double&mety, double bfield,
0258                const reco::Muon* muon, MuonMETInfo& muonMETInfo){
0259
0260   //The tracker has better resolution for pt < 200 GeV
0261   math::XYZTLorentzVector mup4;
0262   if(muon->combinedMuon()->pt() < 200) {
0263     mup4 = math::XYZTLorentzVector(muon->innerTrack()->px(), muon->innerTrack()->py(),
0264                    muon->innerTrack()->pz(), muon->innerTrack()->p());
0265   } else {
0266     mup4 = math::XYZTLorentzVector(muon->combinedMuon()->px(), muon->combinedMuon()->py(),
0267                    muon->combinedMuon()->pz(), muon->combinedMuon()->p());
0268   }
0269
0270   correctMETforMuon(metx, mety, bfield, muon->charge(),
0271             mup4, muon->vertex(),
0272             muonMETInfo);
0273 }
0274
0275 //----------------------------------------------------------------------------
0276
0277 void print_data() {
0278   TFile file("0AC27F01-D260-DD11-8207-0018F3D096EE.root");
0279
0280   fwlite::Event ev(&file);
0281
0282   for( ev.toBegin();
0283        ! ev.atEnd();
0284        ++ev) {
0285     fwlite::Handle<std::vector<reco::Muon> > muonCollection;
0286     muonCollection.getByLabel(ev,"muons");
0287
0288     fwlite::Handle<std::vector<reco::CaloMET> > metCollection;
0289     metCollection.getByLabel(ev,"met");
0290
0291     //this is how one would get the B field in the FW:
0292     //edm::ESHandle<MagneticField> magneticField;
0293     //magneticField.get<IdealMagneticFieldRecord>().get(magneticField);
0294     //how does one do it FWLite?
0295
0296     reco::CaloMETCollection::const_iterator metIt = metCollection->begin();
0297     double metx = metIt->px();
0298     double mety = metIt->py();
0299
0300     //now can access data
0301     for(reco::MuonCollection::const_iterator it = muonCollection->begin();
0302     it != muonCollection->end(); it++) {
0303
0304       if(it->isGlobalMuon() ==0 ) continue;
0305
0306
0307       std::cout << "Before Correcting for the muon of pt: " << it->pt()
0308         << " the MET is: " << sqrt(pow(metx,2)+pow(mety,2)) << endl;
0309
0310       MuonMETInfo muonMETinfo;
0311       muonMETinfo.ecalE = it->isEnergyValid() ? it->calEnergy().emS9 : 0.0;
0312       muonMETinfo.hcalE = it->isEnergyValid() ? it->calEnergy().hadS9: 0.0;
0313       muonMETinfo.hoE   = it->isEnergyValid() ? it->calEnergy().hoS9 : 0.0;
0314       //we don't have the positions at the calo surfaces
0315       //from the TrackAssociator, so use our analytical
0316       //calculation
0317       muonMETinfo.useTkAssociatorPositions = false;
0318
0319       //hardwire the Bfield for now...should get this from
0320       //the EventSetup
0321
0322       correctMETforMuon(metx,mety,3.8112,&(*it), muonMETinfo);
0323
0324       std::cout << "After correcting for the muon, the MET is: "
0325         << sqrt(pow(metx,2)+pow(mety,2)) << endl;
0326
0327     }
0328
0329   }
0330 }