Go to the documentation of this file.00001
00002
00003
00005
00006
00008
00009 #include "FWCore/PluginManager/interface/ModuleDef.h"
00010 #include "FWCore/Framework/interface/MakerMacros.h"
00011
00012 #include "FWCore/Framework/interface/EDProducer.h"
00013 #include "FWCore/Framework/interface/Event.h"
00014 #include "FWCore/Framework/interface/ESHandle.h"
00015 #include "FWCore/Framework/interface/EventSetup.h"
00016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00017
00018 #include "FWCore/Utilities/interface/InputTag.h"
00019 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00020 #include "FWCore/ServiceRegistry/interface/Service.h"
00021
00023
00024 #include "DataFormats/Common/interface/Handle.h"
00025 #include "DataFormats/Common/interface/EDProduct.h"
00026 #include "DataFormats/Common/interface/Ref.h"
00027 #include "DataFormats/Provenance/interface/ProcessHistory.h"
00028 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
00029
00030 #include "DataFormats/Math/interface/LorentzVector.h"
00031 #include "DataFormats/Math/interface/Vector3D.h"
00032
00033 #include <string>
00034
00036
00037 using namespace edm;
00038 using namespace std;
00039
00041
00042
00043
00045
00046 class BeamSpotFromSimProducer : public edm::EDProducer
00047 {
00048 public:
00049
00050
00051 typedef math::XYZPoint Point;
00052 enum { dimension = 7 };
00053 typedef math::Error<dimension>::type CovarianceMatrix;
00054
00055
00057 explicit BeamSpotFromSimProducer(const edm::ParameterSet& iConfig);
00058 virtual ~BeamSpotFromSimProducer();
00059
00060 protected:
00061
00062 private:
00063
00066 virtual void beginRun( edm::Run& run, const edm::EventSetup& iSetup );
00067 virtual void endRun( edm::Run& run, const edm::EventSetup& iSetup );
00068 virtual void produce( edm::Event& iEvent, const edm::EventSetup& iSetup );
00069
00070 double meanX_;
00071 double meanY_;
00072 double meanZ_;
00073
00074 double sigmaX_;
00075 double sigmaY_;
00076 double sigmaZ_;
00077
00078 double dxdz_;
00079 double dydz_;
00080
00081 Point point_;
00082 CovarianceMatrix error_;
00083
00084 };
00085
00086
00088
00089 BeamSpotFromSimProducer::BeamSpotFromSimProducer(edm::ParameterSet const& iConfig)
00090 {
00091 produces< reco::BeamSpot >( "BeamSpot" ).setBranchAlias("BeamSpot");
00092 }
00093
00095
00096 BeamSpotFromSimProducer::~BeamSpotFromSimProducer()
00097 {
00098 }
00099
00101
00102 void BeamSpotFromSimProducer::endRun(edm::Run& run, const edm::EventSetup& iSetup)
00103 {
00105 }
00106
00108
00109 void BeamSpotFromSimProducer::beginRun(edm::Run& run, const edm::EventSetup& iSetup )
00110 {
00111 }
00112
00114
00115 void BeamSpotFromSimProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
00116 {
00117
00118
00119 static bool gotProcessParameterSet=false;
00120
00121 if(!gotProcessParameterSet){
00122
00123 gotProcessParameterSet=true;
00124
00125 edm::ParameterSet ps;
00126
00127
00128
00129
00130 unsigned int iProcess=iEvent.processHistory().size();
00131
00132
00133
00134
00135
00136
00137 if (iProcess>0) iProcess=0;
00138
00139 std::string nameProcess = iEvent.processHistory()[iProcess].processName();
00140
00141
00142
00143
00144 iEvent.getProcessParameterSet(nameProcess,ps);
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155 string type = ps.getParameterSet("VtxSmeared").getParameter<string>("@module_type");
00156
00157 if (type=="HLLHCEvtVtxGenerator") {
00158
00159 double SigmaX = ps.getParameterSet("VtxSmeared").getParameter<double>("SigmaX");
00160 double SigmaY = ps.getParameterSet("VtxSmeared").getParameter<double>("SigmaY");
00161 double SigmaZ = ps.getParameterSet("VtxSmeared").getParameter<double>("SigmaZ");
00162
00163 meanX_ = ps.getParameterSet("VtxSmeared").getParameter<double>("MeanX");
00164 meanY_ = ps.getParameterSet("VtxSmeared").getParameter<double>("MeanY");
00165 meanZ_ = ps.getParameterSet("VtxSmeared").getParameter<double>("MeanZ");
00166
00167 double HalfCrossingAngle = ps.getParameterSet("VtxSmeared").getParameter<double>("HalfCrossingAngle");
00168 double CrabAngle = ps.getParameterSet("VtxSmeared").getParameter<double>("CrabAngle");
00169
00170 static const double sqrtOneHalf=sqrt(0.5);
00171
00172 double sinbeta=sin(CrabAngle);
00173 double cosbeta=cos(CrabAngle);
00174 double cosalpha=cos(HalfCrossingAngle);
00175 double sinamb=sin(HalfCrossingAngle-CrabAngle);
00176 double cosamb=cos(HalfCrossingAngle-CrabAngle);
00177
00178 sigmaX_=sqrtOneHalf*hypot(SigmaZ*sinbeta,SigmaX*cosbeta)/cosalpha;
00179
00180 sigmaY_=sqrtOneHalf*SigmaY;
00181
00182 sigmaZ_=sqrtOneHalf/hypot(sinamb/SigmaX,cosamb/SigmaZ);
00183
00184 dxdz_=0.0;
00185 dydz_=0.0;
00186
00187 point_=Point(meanX_,meanY_,meanZ_);
00188
00189 for (unsigned int j=0; j<7; j++) {
00190 for (unsigned int k=j; k<7; k++) {
00191 error_(j,k) = 0.0;
00192 }
00193 }
00194
00195
00196
00197 error_(0,0)=0.1*sigmaX_;
00198 error_(1,1)=0.1*sigmaY_;
00199 error_(2,2)=0.1*sigmaZ_;
00200
00201
00202 error_(3,3)=0.1*sigmaZ_;
00203 error_(6,6)=0.1*sigmaX_;
00204
00205
00206
00207 error_(4,4)=0.01*sigmaX_/sigmaZ_;
00208 error_(5,5)=error_(4,4);
00209 }
00210 else if (type=="GaussEvtVtxGenerator") {
00211
00212 sigmaX_ = ps.getParameterSet("VtxSmeared").getParameter<double>("SigmaX");
00213 sigmaY_ = ps.getParameterSet("VtxSmeared").getParameter<double>("SigmaY");
00214 sigmaZ_ = ps.getParameterSet("VtxSmeared").getParameter<double>("SigmaZ");
00215
00216 meanX_ = ps.getParameterSet("VtxSmeared").getParameter<double>("MeanX");
00217 meanY_ = ps.getParameterSet("VtxSmeared").getParameter<double>("MeanY");
00218 meanZ_ = ps.getParameterSet("VtxSmeared").getParameter<double>("MeanZ");
00219
00220 dxdz_=0.0;
00221 dydz_=0.0;
00222
00223 point_=Point(meanX_,meanY_,meanZ_);
00224
00225 for (unsigned int j=0; j<7; j++) {
00226 for (unsigned int k=j; k<7; k++) {
00227 error_(j,k) = 0.0;
00228 }
00229 }
00230
00231
00232
00233 error_(0,0)=0.1*sigmaX_;
00234 error_(1,1)=0.1*sigmaY_;
00235 error_(2,2)=0.1*sigmaZ_;
00236
00237
00238 error_(3,3)=0.1*sigmaZ_;
00239 error_(6,6)=0.1*sigmaX_;
00240
00241
00242
00243 error_(4,4)=0.01*sigmaX_/sigmaZ_;
00244 error_(5,5)=error_(4,4);
00245
00246 }
00247 else {
00248 LogError("BeamSpotFromSimProducer") <<"In BeamSpotFromSimProducer type="<<type
00249 <<" don't know what to do!"<<std::endl;
00250
00251 }
00252
00253 }
00254
00255
00257 std::auto_ptr< reco::BeamSpot > BeamSpotForOutput( new reco::BeamSpot(point_,
00258 sigmaZ_,
00259 dxdz_,
00260 dydz_,
00261 sigmaX_,
00262 error_,
00263 reco::BeamSpot::Fake));
00264
00265
00266 iEvent.put( BeamSpotForOutput, "BeamSpot");
00267
00268 }
00269
00270
00271
00272
00273 DEFINE_FWK_MODULE(BeamSpotFromSimProducer);
00274