Merge pull request cms-sw#1 from erc-asymow/WmassNanoProd_10_6_26

Wmass nano prod 10 6 26

Configuration/Applications/python/ConfigBuilder.py

@@ -1715,6 +1715,9 @@ def prepare_USERNANO(self, sequence = "nanoAOD"):
         if sequence in ["nanotpSequence", "nanotpSequenceMC"]:
             self._options.customisation_file.append("PhysicsTools/NanoAOD/nanoTP_cff.customizeNANOTP")
 
+        if sequence in ["nanotpSequenceMC"]:
+            self._options.customisation_file.append("PhysicsTools/NanoAOD/nano_cff.customizeGenLeptonPrecision")
+
     def prepare_EI(self, sequence = None):
         ''' Enrich the schedule with event interpretation '''
         from Configuration.StandardSequences.EventInterpretation import EventInterpretation

MagneticField/ParametrizedEngine/src/ParametrizedMagneticFieldFactory.cc

@@ -10,6 +10,7 @@
 #include "OAEParametrizedMagneticField.h"
 #include "ParabolicParametrizedMagneticField.h"
 #include "PolyFit2DParametrizedMagneticField.h"
+#include "PolyFit3DParametrizedMagneticField.h"
 
 #include "FWCore/Utilities/interface/Exception.h"
 
@@ -32,7 +33,8 @@ ParametrizedMagneticFieldFactory::get(string version, const ParameterSet& parame
     return result;
   } else if (version=="PolyFit3D") {
     // V. Maroussov polynomial fit to mapping data
-    throw cms::Exception("InvalidParameter")<<"PolyFit3D is not supported anymore";
+    std::auto_ptr<MagneticField> result( new PolyFit3DParametrizedMagneticField(parameters));
+    return result;
   } else if (version=="Parabolic"){
     // FIXME implement configurable parameters to be passed to ctor
 //   vector<double> params =  parameters.getParameter<vdouble>("parameters");

MagneticField/ParametrizedEngine/src/PolyFit3DParametrizedMagneticField.cc

@@ -0,0 +1,92 @@
+/** \file
+ *
+ *  \author N. Amapane
+ */
+
+#include "PolyFit3DParametrizedMagneticField.h"
+#include <FWCore/ParameterSet/interface/ParameterSet.h>
+#include <FWCore/MessageLogger/interface/MessageLogger.h>
+
+#include "BFit3D.h"
+
+
+using namespace std;
+using namespace magfieldparam;
+
+PolyFit3DParametrizedMagneticField::PolyFit3DParametrizedMagneticField(double bVal) : 
+  theParam(new BFit3D())
+{
+  theParam->SetField(bVal);
+}
+
+
+PolyFit3DParametrizedMagneticField::PolyFit3DParametrizedMagneticField(const edm::ParameterSet& parameters) : theParam(new BFit3D()) {
+  theParam->SetField(parameters.getParameter<double>("BValue"));
+
+  // Additional options (documentation by Vassili):
+
+  // By default, the package accepts signed value of "r". That means,
+  // one can cross r=0 and orientation of the coordinate "orts"
+  // e_r and e_phi will not be flipped over.
+  // In other words for an r<0 the e_r points inward, in the direction r=0.
+  // This is a "natural" mode. However, the default behavior may be
+  // changed by the call
+
+  //  theParam->UseSignedRad(false);
+
+  // In that case with crossing of r=0 e_r and e_phi will be flipped in
+  // such a way that e_r always points outward. In other words instead of 
+  // (r<0, phi) the (abs(r), phi+PI) will be used in this mode.
+
+  // The expansion coefficients for a nominal field in between measurement
+  // field values (2.0T, 3.5T, 3.8T and 4.0T) by default are calculated by
+  // means of a linear piecewise interpolation. Another provided
+  // interpolation mode is cubic spline. This mode can be switched
+  // on by the call:
+
+  //  theParam->UseSpline(true);
+
+  // From practical point of view the use of spline interpolation doesn't
+  // change much, but it makes the coefficients' behavior a bit more
+  // physical at very low/high field values.
+}
+
+
+PolyFit3DParametrizedMagneticField::~PolyFit3DParametrizedMagneticField() {
+  delete theParam;
+}
+
+
+GlobalVector
+PolyFit3DParametrizedMagneticField::inTesla(const GlobalPoint& gp) const {
+
+  if (isDefined(gp)) {
+    return inTeslaUnchecked(gp);
+  } else {
+    edm::LogWarning("MagneticField|FieldOutsideValidity") << " Point " << gp << " is outside the validity region of PolyFit3DParametrizedMagneticField";
+    return GlobalVector();
+  }
+}
+
+GlobalVector
+PolyFit3DParametrizedMagneticField::inTeslaUnchecked(const GlobalPoint& gp) const {
+  double Br, Bz, Bphi;
+  theParam->GetField(gp.perp()/100., gp.z()/100., gp.phi(),
+		     Br, Bz, Bphi);
+
+  double cosphi = cos(gp.phi());
+  double sinphi = sin(gp.phi());
+
+  return GlobalVector(Br*cosphi - Bphi*sinphi,
+		      Br*sinphi + Bphi*cosphi, 
+		      Bz);  
+}
+
+bool
+PolyFit3DParametrizedMagneticField::isDefined(const GlobalPoint& gp) const {
+  double z = fabs(gp.z());
+  double r = gp.perp();
+  //"rectangle" |z|<3.5, r<1.9 _except_ the "corners" |z|+2.5*r>6.7, everything in meters
+  if (z>350. || r>190 || z+2.5*r>670.) return false;
+  return true;
+}

MagneticField/ParametrizedEngine/src/PolyFit3DParametrizedMagneticField.h

@@ -0,0 +1,40 @@
+#ifndef PolyFit3DParametrizedMagneticField_h
+#define PolyFit3DParametrizedMagneticField_h
+
+/** \class PolyFit3DParametrizedMagneticField
+ *
+ *  Magnetic Field engine wrapper for V. Maroussov's 3D parametrization
+ *  of the MT data.
+ *
+ *  \author N. Amapane
+ */
+
+#include "MagneticField/Engine/interface/MagneticField.h"
+
+namespace edm { class ParameterSet; }
+namespace magfieldparam { class BFit3D; }
+
+
+class PolyFit3DParametrizedMagneticField : public MagneticField {
+ public:
+  /// Constructor. Fitted bVal for the nominal currents are:
+  /// 2.0216; 3.5162;  3.8114; 4.01242188708911
+  PolyFit3DParametrizedMagneticField(double bVal = 3.8114);
+
+  /// Constructor. Parameters taken from a PSet
+  PolyFit3DParametrizedMagneticField(const edm::ParameterSet& parameters);
+
+  /// Destructor
+  virtual ~PolyFit3DParametrizedMagneticField();
+
+  GlobalVector inTesla (const GlobalPoint& gp) const;
+
+  GlobalVector inTeslaUnchecked (const GlobalPoint& gp) const;
+
+  bool isDefined(const GlobalPoint& gp) const;
+
+ private:
+  magfieldparam::BFit3D* theParam;
+};
+#endif
+

PhysicsTools/NanoAOD/interface/FlattenedValueMapVectorTableProducer.h

@@ -84,10 +84,10 @@ class FlattenedValueMapVectorTableProducer : public edm::stream::EDProducer<> {
                 iEvent.getByToken(intVecMaps_[i], vmap);
                 const auto& results = readVals(*vmap, objs, sizes);
 
-                auto intsizetab = std::make_unique<nanoaod::FlatTable>(objs.size(), this->name_+intNames_[i]+"_Counts", false, false);
+                auto intsizetab = std::make_unique<nanoaod::FlatTable>(objs.size(), this->name_ + "_" + intNames_[i]+"_Counts", false, false);
                 intsizetab->template addColumn<int>("", sizes, "Number of entries per object", nanoaod::FlatTable::IntColumn, -1);
                 intsizetab->setDoc(doc_);
-                auto intvectab = std::make_unique<nanoaod::FlatTable>(results.size(), this->name_+intNames_[i], false, false);
+                auto intvectab = std::make_unique<nanoaod::FlatTable>(results.size(), this->name_+ "_" + intNames_[i], false, false);
                 intvectab->template addColumn<int>("Vals", results, intDocs_[i], nanoaod::FlatTable::IntColumn, intPrecisions_[i]);
                 intvectab->setDoc(doc_);
 
@@ -100,10 +100,10 @@ class FlattenedValueMapVectorTableProducer : public edm::stream::EDProducer<> {
                 iEvent.getByToken(floatVecMaps_[i], vmap);
                 const auto& results = readVals(*vmap, objs, sizes);
 
-                auto floatsizetab = std::make_unique<nanoaod::FlatTable>(objs.size(), this->name_+floatNames_[i]+"_Counts", false, false);
+                auto floatsizetab = std::make_unique<nanoaod::FlatTable>(objs.size(), this->name_ + "_" + floatNames_[i]+"_Counts", false, false);
                 floatsizetab->template addColumn<int>("", sizes, "Number of entries per object", nanoaod::FlatTable::IntColumn, -1);
                 floatsizetab->setDoc(doc_);
-                auto floatvectab = std::make_unique<nanoaod::FlatTable>(results.size(), this->name_+floatNames_[i], false, false);
+                auto floatvectab = std::make_unique<nanoaod::FlatTable>(results.size(), this->name_ + "_" + floatNames_[i], false, false);
                 floatvectab->template addColumn<float>("Vals", results, floatDocs_[i], nanoaod::FlatTable::FloatColumn, floatPrecisions_[i]);
                 floatvectab->setDoc(doc_);
 

PhysicsTools/NanoAOD/interface/SimpleFlatTableProducer.h

@@ -10,6 +10,7 @@
 
 #include <vector>
 #include <boost/ptr_container/ptr_vector.hpp>
+#include <iostream>
 
 template<typename T, typename TProd>
 class SimpleFlatTableProducerBase : public edm::stream::EDProducer<> {
@@ -89,10 +90,13 @@ class SimpleFlatTableProducerBase : public edm::stream::EDProducer<> {
                     void fill(std::vector<const T *> selobjs, nanoaod::FlatTable & out) const override {
                         std::vector<ValType> vals(selobjs.size());
                         for (unsigned int i = 0, n = vals.size(); i < n; ++i) {
-                            if(this->precision_ == -2){
-                            vals[i] = MiniFloatConverter::reduceMantissaToNbitsRounding(func_(*selobjs[i]),precisionFunc_(*selobjs[i]));
+                            ValType val =  func_(*selobjs[i]);
+                            if(this->precision_ == -2) {
+                                auto prec =  precisionFunc_(*selobjs[i]);
+                                vals[i] = prec > 0 ? MiniFloatConverter::reduceMantissaToNbitsRounding(val, prec) : val;
                             }
-                            else vals[i] = func_(*selobjs[i]);
+                            else 
+                                vals[i] = val;
                         }
                         out.template addColumn<ValType>(this->name_, vals, this->doc_, this->type_,this->precision_);
                     }

PhysicsTools/NanoAOD/python/generalTracks_cff.py

@@ -20,6 +20,9 @@
         trackOriginalAlgo = Var('originalAlgo()', 'int', precision=-1, doc='Track original algo enum'),
         qualityMask = Var('qualityMask()', 'int', precision=-1, doc='Quality mask of the track.'),
         extraIdx = Var('extra().key()', 'int', precision=-1, doc='Index of the TrackExtra in the original collection'),
+        vx = Var('vx', 'float', precision=-1, doc='Track X position'),
+        vy = Var('vy', 'float', precision=-1, doc='Track Y position'),
+        vz = Var('vz', 'float', precision=-1, doc='Track Z position'),
     ),
 )