Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:19:25 

0001 #include <vector>
0002 #include "TMath.h"
0003 #include "Math/Point3D.h"
0004 #include "Math/GenVector/CoordinateSystemTags.h"
0005 #include "Math/PtEtaPhiE4D.h"
0006 #include "Math/PtEtaPhiM4D.h"
0007 #include "Math/LorentzVector.h"
0008 
0009 namespace math {
0010   typedef ROOT::Math::PositionVector3D<ROOT::Math::Cartesian3D<double> > XYZPointD;
0011   typedef XYZPointD XYZPoint;
0012   typedef ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > XYZTLorentzVectorD;
0013   typedef XYZTLorentzVectorD XYZTLorentzVector;
0014 }
0015 
0016 struct MuonMETInfo {
0017 
0018   float ecalE;
0019   float hcalE;
0020   float hoE;
0021   math::XYZPoint ecalPos;
0022   math::XYZPoint hcalPos;
0023   math::XYZPoint hoPos;
0024   bool useAverage;
0025   //if using in FWLite, this should be false
0026   bool useTkAssociatorPositions;
0027   bool useHO;
0028 
0029   MuonMETInfo():
0030   ecalE(0), hcalE(0), hoE(0),
0031   ecalPos(0,0,0),hcalPos(0,0,0), hoPos(0,0,0),
0032   useAverage(0), useTkAssociatorPositions(0),useHO(0){ }
0033 
0034 };
0035 
0036 
0037 
0038 //----------------------------------------------------------------------------
0039 void correctMETforMuon(double& metx, double& mety, double bfield, int muonCharge,
0040                math::XYZTLorentzVector muonP4, math::XYZPoint muonVertex,
0041                MuonMETInfo& muonMETInfo) {
0042 
0043 
0044   
0045   double mu_p     = muonP4.P();
0046   double mu_pt    = muonP4.Pt();
0047   double mu_phi   = muonP4.Phi();
0048   double mu_eta   = muonP4.Eta();
0049   double mu_vz    = muonVertex.z()/100.;
0050   double mu_pz    = muonP4.Pz();
0051     
0052   
0053   //add in the muon's pt
0054   metx -= mu_pt*cos(mu_phi);
0055   mety -= mu_pt*sin(mu_phi);
0056 
0057   double ecalPhi, ecalTheta;
0058   double hcalPhi, hcalTheta;
0059   double hoPhi, hoTheta;
0060   
0061 
0062   //should always be false for FWLite
0063   //unless you want to supply co-ordinates at 
0064   //the calorimeter sub-detectors yourself
0065   if(muonMETInfo.useTkAssociatorPositions) {
0066     ecalPhi   = muonMETInfo.ecalPos.Phi();
0067     ecalTheta = muonMETInfo.ecalPos.Theta();
0068     hcalPhi   = muonMETInfo.hcalPos.Phi();
0069     hcalTheta = muonMETInfo.hcalPos.Theta();
0070     hoPhi     = muonMETInfo.hoPos.Phi();
0071     hoTheta   = muonMETInfo.hoPos.Theta();
0072   } else {
0073     
0074     /*
0075       use the analytical solution for the
0076       intersection of a helix with a cylinder
0077       to find the positions of the muon
0078       at the various calo surfaces
0079     */
0080     double rEcal = 1.290;
0081     double rHcal = 1.9;
0082     double rHo   = 3.82;
0083     if(abs(mu_eta) > 0.3) rHo = 4.07;
0084     //distance from the center of detector to face of Ecal
0085     double zFaceEcal = 3.209;
0086     if(mu_eta < 0 ) zFaceEcal = -1*zFaceEcal;
0087     //distance from the center of detector to face of Hcal
0088     double zFaceHcal = 3.88;
0089     if(mu_eta < 0 ) zFaceHcal = -1*zFaceHcal;    
0090     
0091     //now we have to get Phi
0092     //bending radius of the muon (units are meters)
0093     double bendr = mu_pt*1000/(300*bfield);
0094 
0095     double tb_ecal = TMath::ACos(1-rEcal*rEcal/(2*bendr*bendr)); //helix time interval parameter
0096     double tb_hcal = TMath::ACos(1-rHcal*rHcal/(2*bendr*bendr)); //helix time interval parameter
0097     double tb_ho   = TMath::ACos(1-rHo*rHo/(2*bendr*bendr));     //helix time interval parameter
0098     double xEcal,yEcal,zEcal;
0099     double xHcal,yHcal,zHcal; 
0100     double xHo, yHo,zHo;
0101     //Ecal
0102     //in the barrel and if not a looper
0103     if(fabs(mu_pz*bendr*tb_ecal/mu_pt+mu_vz) < fabs(zFaceEcal) && rEcal < 2*bendr) { 
0104       xEcal = bendr*(TMath::Sin(tb_ecal+mu_phi)-TMath::Sin(mu_phi));
0105       yEcal = bendr*(-TMath::Cos(tb_ecal+mu_phi)+TMath::Cos(mu_phi));
0106       zEcal = bendr*tb_ecal*mu_pz/mu_pt + mu_vz;
0107     } else { //endcap 
0108       if(mu_pz > 0) {
0109         double te_ecal = (fabs(zFaceEcal) - mu_vz)*mu_pt/(bendr*mu_pz);
0110         xEcal = bendr*(TMath::Sin(te_ecal+mu_phi) - TMath::Sin(mu_phi));
0111         yEcal = bendr*(-TMath::Cos(te_ecal+mu_phi) + TMath::Cos(mu_phi));
0112         zEcal = fabs(zFaceEcal);
0113       } else {
0114         double te_ecal = -(fabs(zFaceEcal) + mu_vz)*mu_pt/(bendr*mu_pz);
0115         xEcal = bendr*(TMath::Sin(te_ecal+mu_phi) - TMath::Sin(mu_phi));
0116         yEcal = bendr*(-TMath::Cos(te_ecal+mu_phi) + TMath::Cos(mu_phi));
0117         zEcal = -fabs(zFaceEcal);
0118       }
0119     }
0120 
0121     //Hcal
0122     if(fabs(mu_pz*bendr*tb_hcal/mu_pt+mu_vz) < fabs(zFaceHcal) && rEcal < 2*bendr) { //in the barrel
0123       xHcal = bendr*(TMath::Sin(tb_hcal+mu_phi)-TMath::Sin(mu_phi));
0124       yHcal = bendr*(-TMath::Cos(tb_hcal+mu_phi)+TMath::Cos(mu_phi));
0125       zHcal = bendr*tb_hcal*mu_pz/mu_pt + mu_vz;
0126     } else { //endcap 
0127       if(mu_pz > 0) {
0128         double te_hcal = (fabs(zFaceHcal) - mu_vz)*mu_pt/(bendr*mu_pz);
0129         xHcal = bendr*(TMath::Sin(te_hcal+mu_phi) - TMath::Sin(mu_phi));
0130         yHcal = bendr*(-TMath::Cos(te_hcal+mu_phi) + TMath::Cos(mu_phi));
0131         zHcal = fabs(zFaceHcal);
0132       } else {
0133         double te_hcal = -(fabs(zFaceHcal) + mu_vz)*mu_pt/(bendr*mu_pz);
0134         xHcal = bendr*(TMath::Sin(te_hcal+mu_phi) - TMath::Sin(mu_phi));
0135         yHcal = bendr*(-TMath::Cos(te_hcal+mu_phi) + TMath::Cos(mu_phi));
0136         zHcal = -fabs(zFaceHcal);
0137       }
0138     }
0139     
0140     //Ho - just the barrel
0141     xHo = bendr*(TMath::Sin(tb_ho+mu_phi)-TMath::Sin(mu_phi));
0142     yHo = bendr*(-TMath::Cos(tb_ho+mu_phi)+TMath::Cos(mu_phi));
0143     zHo = bendr*tb_ho*mu_pz/mu_pt + mu_vz;  
0144     
0145     ecalTheta = TMath::ACos(zEcal/sqrt(pow(xEcal,2) + pow(yEcal,2)+pow(zEcal,2)));
0146     ecalPhi   = atan2(yEcal,xEcal);
0147     hcalTheta = TMath::ACos(zHcal/sqrt(pow(xHcal,2) + pow(yHcal,2)+pow(zHcal,2)));
0148     hcalPhi   = atan2(yHcal,xHcal);
0149     hoTheta   = TMath::ACos(zHo/sqrt(pow(xHo,2) + pow(yHo,2)+pow(zHo,2)));
0150     hoPhi     = atan2(yHo,xHo);
0151 
0152     //2d radius in x-y plane
0153     double r2dEcal = sqrt(pow(xEcal,2)+pow(yEcal,2));
0154     double r2dHcal = sqrt(pow(xHcal,2)+pow(yHcal,2));
0155     double r2dHo   = sqrt(pow(xHo,2)  +pow(yHo,2));
0156     
0157     /*
0158       the above prescription is for right handed helicies only
0159       Positively charged muons trace a left handed helix
0160       so we correct for that 
0161     */
0162     if(muonCharge > 0) {
0163          
0164       //Ecal
0165       double dphi = mu_phi - ecalPhi;
0166       if(fabs(dphi) > TMath::Pi()) 
0167         dphi = 2*TMath::Pi() - fabs(dphi);
0168       ecalPhi = mu_phi - fabs(dphi);
0169       if(fabs(ecalPhi) > TMath::Pi()) {
0170         double temp = 2*TMath::Pi() - fabs(ecalPhi);
0171         ecalPhi = -1*temp*ecalPhi/fabs(ecalPhi);
0172       }
0173       xEcal = r2dEcal*TMath::Cos(ecalPhi);
0174       yEcal = r2dEcal*TMath::Sin(ecalPhi);
0175       
0176       //Hcal
0177       dphi = mu_phi - hcalPhi;
0178       if(fabs(dphi) > TMath::Pi()) 
0179         dphi = 2*TMath::Pi() - fabs(dphi);
0180       hcalPhi = mu_phi - fabs(dphi);
0181       if(fabs(hcalPhi) > TMath::Pi()) {
0182         double temp = 2*TMath::Pi() - fabs(hcalPhi);
0183         hcalPhi = -1*temp*hcalPhi/fabs(hcalPhi);
0184       }
0185       xHcal = r2dHcal*TMath::Cos(hcalPhi);
0186       yHcal = r2dHcal*TMath::Sin(hcalPhi);
0187          
0188       //Ho
0189       dphi = mu_phi - hoPhi;
0190       if(fabs(dphi) > TMath::Pi())
0191         dphi = 2*TMath::Pi() - fabs(dphi);
0192       hoPhi = mu_phi - fabs(dphi);
0193       if(fabs(hoPhi) > TMath::Pi()) {
0194         double temp = 2*TMath::Pi() - fabs(hoPhi);
0195         hoPhi = -1*temp*hoPhi/fabs(hoPhi);
0196       }
0197       xHo = r2dHo*TMath::Cos(hoPhi);
0198       yHo = r2dHo*TMath::Sin(hoPhi);
0199       
0200     }
0201   }
0202   
0203   if(!muonMETInfo.useHO) muonMETInfo.hoE = 0;
0204   
0205   //for isolated muons
0206   if(!muonMETInfo.useAverage) {
0207     
0208     double mu_Ex =  muonMETInfo.ecalE*sin(ecalTheta)*cos(ecalPhi)
0209       + muonMETInfo.hcalE*sin(hcalTheta)*cos(hcalPhi)
0210       + muonMETInfo.hoE*sin(hoTheta)*cos(hoPhi);
0211     double mu_Ey =  muonMETInfo.ecalE*sin(ecalTheta)*sin(ecalPhi)
0212       + muonMETInfo.hcalE*sin(hcalTheta)*sin(hcalPhi)
0213       + muonMETInfo.hoE*sin(hoTheta)*sin(hoPhi);
0214 
0215     metx += mu_Ex;
0216     mety += mu_Ey;
0217     
0218   } else { //non-isolated muons - derive the correction
0219     
0220     //dE/dx in matter for iron:
0221     //-(11.4 + 0.96*fabs(log(p0*2.8)) + 0.033*p0*(1.0 - pow(p0, -0.33)) )*1e-3
0222     //from http://cmslxr.fnal.gov/lxr/source/TrackPropagation/SteppingHelixPropagator/src/SteppingHelixPropagator.ccyes,
0223     //line ~1100
0224     //normalisation is at 50 GeV
0225     double dEdx_normalization = -(11.4 + 0.96*fabs(log(50*2.8)) + 0.033*50*(1.0 - pow(50, -0.33)) )*1e-3;
0226     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;
0227     
0228     double temp = 0.0;
0229     
0230     if(muonMETInfo.useHO) {
0231       //for the Towers, with HO
0232       if(fabs(mu_eta) < 0.2)
0233         temp = 2.75*(1-0.00003*mu_p);
0234       if(fabs(mu_eta) > 0.2 && fabs(mu_eta) < 1.0)
0235         temp = (2.38+0.0144*fabs(mu_eta))*(1-0.0003*mu_p);
0236       if(fabs(mu_eta) > 1.0 && fabs(mu_eta) < 1.3)
0237         temp = 7.413-5.12*fabs(mu_eta);
0238       if(fabs(mu_eta) > 1.3)
0239         temp = 2.084-0.743*fabs(mu_eta);
0240     } else {
0241       if(fabs(mu_eta) < 1.0)
0242         temp = 2.33*(1-0.0004*mu_p);
0243       if(fabs(mu_eta) > 1.0 && fabs(mu_eta) < 1.3)
0244         temp = (7.413-5.12*fabs(mu_eta))*(1-0.0003*mu_p);
0245       if(fabs(mu_eta) > 1.3)
0246         temp = 2.084-0.743*fabs(mu_eta);
0247     }
0248       
0249     double dep = temp*dEdx_normalization/dEdx_numerator;
0250     if(dep < 0.5 ) dep = 0.0;
0251     
0252     //use the average phi of the 3 subdetectors
0253     if(fabs(mu_eta) < 1.3) {
0254       metx += dep*cos((ecalPhi+hcalPhi+hoPhi)/3);
0255       mety += dep*sin((ecalPhi+hcalPhi+hoPhi)/3);
0256     } else {
0257       metx += dep*cos( (ecalPhi+hcalPhi)/2);
0258       mety += dep*cos( (ecalPhi+hcalPhi)/2);
0259     }
0260   }
0261 
0262   
0263 }
0264 //----------------------------------------------------------------------------
0265
[ 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!