CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_6_1_2_SLHC4_patch1/src/RecoVertex/BeamSpotProducer/plugins/BeamSpotFromSimProducer.cc

Go to the documentation of this file.
00001 
00002 //  Producer by Anders  //
00003 //    Jan 2013 @ CU    //
00005 
00006 
00008 // FRAMEWORK HEADERS
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 // DATA FORMATS HEADERS
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 // NAMESPACES
00037 using namespace edm;
00038 using namespace std;
00039 
00041 //                          //
00042 //     CLASS DEFINITION     //
00043 //                          //
00045 
00046 class BeamSpotFromSimProducer : public edm::EDProducer
00047 {
00048 public:
00049 
00050   // point in the space
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 // CONSTRUCTOR
00089 BeamSpotFromSimProducer::BeamSpotFromSimProducer(edm::ParameterSet const& iConfig) // :   config(iConfig)
00090 {
00091   produces< reco::BeamSpot >( "BeamSpot" ).setBranchAlias("BeamSpot");
00092 }
00093 
00095 // DESTRUCTOR
00096 BeamSpotFromSimProducer::~BeamSpotFromSimProducer()
00097 {
00098 }  
00099 
00101 // END JOB
00102 void BeamSpotFromSimProducer::endRun(edm::Run& run, const edm::EventSetup& iSetup)
00103 {
00105 }
00106 
00108 // BEGIN JOB
00109 void BeamSpotFromSimProducer::beginRun(edm::Run& run, const edm::EventSetup& iSetup )
00110 {
00111 }
00112 
00114 // PRODUCE
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     //fetch the real process name using the processHistory accessor of 
00127     //the edm::Event, then find the     
00128     // last-1 entry and ask for it's processName
00129 
00130     unsigned int iProcess=iEvent.processHistory().size();
00131     //std::cout << "Process History size = "<<iProcess<<std::endl;
00132     //for (unsigned int i=0;i<iProcess;i++){
00133     //  std::cout << "Process "<<i<<" name: "
00134     //          <<iEvent.processHistory()[i].processName()<<std::endl;
00135     //}
00136 
00137     if (iProcess>0) iProcess=0;
00138 
00139     std::string nameProcess = iEvent.processHistory()[iProcess].processName();
00140     
00141     // Now ask the edm::Event for the top level parameter set of this process, 
00142     //return it in ps
00143 
00144     iEvent.getProcessParameterSet(nameProcess,ps);
00145 
00146     // The VtxSmeared has a parameter name "SigmaZ" this will 
00147     // retrieve it.
00148 
00149     //cout << "Will print out ps:" << endl;
00150     //cout << ps << endl;
00151     //cout << "Done printing ps" << endl;
00152 
00153     //    double SigmaZ = ps.getParameterSet("HLLHCVtxSmearingParameters").getParameter<double>("SigmaZ");
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       //arbitrarily set position errors to 1/10 of width.
00197       error_(0,0)=0.1*sigmaX_;
00198       error_(1,1)=0.1*sigmaY_;
00199       error_(2,2)=0.1*sigmaZ_;
00200     
00201       //arbitrarily set width errors to 1/10 of width.
00202       error_(3,3)=0.1*sigmaZ_;
00203       error_(6,6)=0.1*sigmaX_;
00204       
00205       //arbitrarily set error on beam axis direction to 1/100 of
00206       //beam aspect ratio    
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       //arbitrarily set position errors to 1/10 of width.
00233       error_(0,0)=0.1*sigmaX_;
00234       error_(1,1)=0.1*sigmaY_;
00235       error_(2,2)=0.1*sigmaZ_;
00236     
00237       //arbitrarily set width errors to 1/10 of width.
00238       error_(3,3)=0.1*sigmaZ_;
00239       error_(6,6)=0.1*sigmaX_;
00240       
00241       //arbitrarily set error on beam axis direction to 1/100 of
00242       //beam aspect ratio    
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 // // DEFINE THIS AS A PLUG-IN
00273 DEFINE_FWK_MODULE(BeamSpotFromSimProducer);
00274