Merge pull request #48594 from tvami/UpCosmics_master

Add new functionality for upward going case in the Cosmics generator

Author: cmsbuild

Configuration/Generator/python/UndergroundCosmicMu_cfi.py

@@ -38,6 +38,7 @@
     MultiMuonFileName = cms.string("CORSIKAmultiMuon.root"),
     MultiMuonFileFirstEvent = cms.int32(1),
     MultiMuonNmin = cms.int32(2),              
+    SurfaceDepth = cms.double(0.0)
 )
 
 

Configuration/Generator/python/UndergroundCosmicSPLooseMu_cfi.py

@@ -18,6 +18,7 @@
     MaxEnu = cms.double(10000.),                
     NuProdAlt = cms.double(7.5e6),                       
     AcptAllMu = cms.bool(False), 
+    SurfaceDepth = cms.double(0.0),
     ElossScaleFactor = cms.double(1.0),
     RadiusOfTarget = cms.double(8000.0),
     ZDistOfTarget = cms.double(15000.0),

GeneratorInterface/CosmicMuonGenerator/interface/CosMuoGenProducer.h

@@ -77,6 +77,9 @@ namespace edm {
 
     bool AllMu;  //Accepting All Muons regardeless of direction
 
+    /// depth of the surface where upgoing muon is generated relative to bottom of CMS
+    double SurfaceDepthToCMS;
+
     // external cross section and filter efficiency
     double extCrossSect;
     double extFilterEff;

GeneratorInterface/CosmicMuonGenerator/interface/CosmicMuonGenerator.h

@@ -220,6 +220,9 @@ class CosmicMuonGenerator {
 
   bool AcptAllMu;  //Accepting All Muons regardeless of direction
 
+  /// depth of the surface where upgoing muon is generated relative to bottom of CMS
+  double SurfaceDepth;
+
   // random number generator
   CLHEP::HepRandomEngine* RanGen;
   bool delRanGen;
@@ -270,6 +273,7 @@ class CosmicMuonGenerator {
   void setMaxEnu(double MaxEn);
   void setNuProdAlt(double NuPrdAlt);
   void setAcptAllMu(bool AllMu);
+  void setSurfaceDepth(double SurfaceDepthToCMS);
 
   // initialize the generator
   void setRandomEngine(CLHEP::HepRandomEngine* v);

GeneratorInterface/CosmicMuonGenerator/python/CMSCGENproducer_cfi.py

@@ -41,5 +41,6 @@
     MaxEnu = cms.double(10000.),                
     NuProdAlt = cms.double(7.5e6),                       
     AcptAllMu = cms.bool(False),                
+    SurfaceDepth = cms.double(0.0),
     Verbosity = cms.bool(False)
 )

GeneratorInterface/CosmicMuonGenerator/src/CosMuoGenProducer.cc

@@ -43,6 +43,7 @@ edm::CosMuoGenProducer::CosMuoGenProducer(const ParameterSet& pset)
       MaxEn(pset.getParameter<double>("MaxEnu")),
       NuPrdAlt(pset.getParameter<double>("NuProdAlt")),
       AllMu(pset.getParameter<bool>("AcptAllMu")),
+      SurfaceDepthToCMS(pset.getParameter<double>("SurfaceDepth")),
       extCrossSect(pset.getUntrackedParameter<double>("crossSection", -1.)),
       extFilterEff(pset.getUntrackedParameter<double>("filterEfficiency", -1.)),
       cmVerbosity_(pset.getParameter<bool>("Verbosity")),
@@ -91,6 +92,7 @@ edm::CosMuoGenProducer::CosMuoGenProducer(const ParameterSet& pset)
   CosMuoGen->setMaxEnu(MaxEn);
   CosMuoGen->setNuProdAlt(NuPrdAlt);
   CosMuoGen->setAcptAllMu(AllMu);
+  CosMuoGen->setSurfaceDepth(SurfaceDepthToCMS);
   produces<HepMCProduct>("unsmeared");
   produces<GenEventInfoProduct>();
   produces<GenRunInfoProduct, edm::Transition::EndRun>();
@@ -142,17 +144,21 @@ void edm::CosMuoGenProducer::produce(Event& e, const edm::EventSetup& es) {
     std::cout << "CosMuoGen->Id_at=" << CosMuoGen->Id_at << "  CosMuoGen->Vx_at=" << CosMuoGen->Vx_at
               << "  CosMuoGen->Vy_at=" << CosMuoGen->Vy_at << "  CosMuoGen->Vz_at=" << CosMuoGen->Vz_at
               << "  CosMuoGen->T0_at=" << CosMuoGen->T0_at << std::endl;
-    std::cout << "  Px=" << CosMuoGen->Px_at << "  Py=" << CosMuoGen->Py_at << "  Pz=" << CosMuoGen->Pz_at << std::endl;
+    std::cout << "  Px=" << CosMuoGen->Px_at << "  Py=" << CosMuoGen->Py_at << "  Pz=" << CosMuoGen->Pz_at
+              << "  E_at=" << CosMuoGen->E_at << std::endl;
+
     for (unsigned int i = 0; i < CosMuoGen->Id_sf.size(); ++i) {
       std::cout << "Id_sf[" << i << "]=" << CosMuoGen->Id_sf[i] << "  Vx_sf[" << i << "]=" << CosMuoGen->Vx_sf[i]
                 << "  Vy_sf=" << CosMuoGen->Vy_sf[i] << "  Vz_sf=" << CosMuoGen->Vz_sf[i]
                 << "  T0_sf=" << CosMuoGen->T0_sf[i] << "  Px_sf=" << CosMuoGen->Px_sf[i]
-                << "  Py_sf=" << CosMuoGen->Py_sf[i] << "  Pz_sf=" << CosMuoGen->Pz_sf[i] << std::endl;
+                << "  Py_sf=" << CosMuoGen->Py_sf[i] << "  Pz_sf=" << CosMuoGen->Pz_sf[i]
+                << "  E_sf=" << CosMuoGen->E_sf[i] << std::endl;
       std::cout << "phi_sf=" << atan2(CosMuoGen->Px_sf[i], CosMuoGen->Pz_sf[i]) << std::endl;
       std::cout << "Id_ug[" << i << "]=" << CosMuoGen->Id_ug[i] << "  Vx_ug[" << i << "]=" << CosMuoGen->Vx_ug[i]
                 << "  Vy_ug=" << CosMuoGen->Vy_ug[i] << "  Vz_ug=" << CosMuoGen->Vz_ug[i]
                 << "  T0_ug=" << CosMuoGen->T0_ug[i] << "  Px_ug=" << CosMuoGen->Px_ug[i]
-                << "  Py_ug=" << CosMuoGen->Py_ug[i] << "  Pz_ug=" << CosMuoGen->Pz_ug[i] << std::endl;
+                << "  Py_ug=" << CosMuoGen->Py_ug[i] << "  Pz_ug=" << CosMuoGen->Pz_ug[i]
+                << "  E_ug=" << CosMuoGen->E_ug[i] << std::endl;
       std::cout << "phi_ug=" << atan2(CosMuoGen->Px_ug[i], CosMuoGen->Pz_ug[i]) << std::endl;
       ;
     }

GeneratorInterface/CosmicMuonGenerator/src/CosmicMuonGenerator.cc

@@ -49,7 +49,7 @@ void CosmicMuonGenerator::initialize(CLHEP::HepRandomEngine* rng) {
       std::cout << "  radius of sphere  around  target = " << Target3dRadius << " mm" << std::endl;
 
     if (MinTheta > 90. * Deg2Rad)  //upgoing muons from neutrinos
-      SurfaceRadius = (RadiusCMS) * (-tan(MinTheta)) + MinStepSize;
+      SurfaceRadius = (RadiusCMS + SurfaceDepth) * (-tan(MinTheta)) + MinStepSize;
     else
       SurfaceRadius = (SurfaceOfEarth + PlugWidth + RadiusTargetEff) * tan(MaxTheta) + Target3dRadius;
     if (Debug)
@@ -174,7 +174,7 @@ void CosmicMuonGenerator::nextEvent() {
 
     double Vy;
     if (MinTheta > 90. * Deg2Rad)  //upgoing muons from neutrinos
-      Vy = -RadiusCMS;
+      Vy = -(RadiusCMS + SurfaceDepth);
     else
       Vy = SurfaceOfEarth + PlugWidth;  // [mm]
 
@@ -195,7 +195,7 @@ void CosmicMuonGenerator::nextEvent() {
     // if angles are ok, propagate to target
     if (goodOrientation()) {
       if (MinTheta > 90. * Deg2Rad)  //upgoing muons from neutrinos
-        OneMuoEvt.propagate(0., RadiusOfTarget, ZDistOfTarget, ZCentrOfTarget, TrackerOnly, MTCCHalf);
+        OneMuoEvt.propagate(ElossScaleFactor, RadiusOfTarget, ZDistOfTarget, ZCentrOfTarget, TrackerOnly, MTCCHalf);
       else
         OneMuoEvt.propagate(ElossScaleFactor, RadiusOfTarget, ZDistOfTarget, ZCentrOfTarget, TrackerOnly, MTCCHalf);
     }
@@ -990,7 +990,7 @@ bool CosmicMuonGenerator::goodOrientation() {
 
   double rVY;
   if (MinTheta > 90. * Deg2Rad)  //upgoing muons from neutrinos
-    rVY = -sqrt(RxzV * RxzV + RadiusCMS * RadiusCMS);
+    rVY = -sqrt(RxzV * RxzV + (RadiusCMS + SurfaceDepth) * (RadiusCMS + SurfaceDepth));
   else
     rVY = sqrt(RxzV * RxzV + (SurfaceOfEarth + PlugWidth) * (SurfaceOfEarth + PlugWidth));
 
@@ -1011,6 +1011,13 @@ bool CosmicMuonGenerator::goodOrientation() {
   double dTheta = Pi;
   if (std::fabs(rVY) > Target3dRadius)
     dTheta = asin(Target3dRadius / std::fabs(rVY));
+
+  if (MinTheta > 90. * Deg2Rad) {
+    // upgoing muon's vertex could be outside the target sphere
+    ThetaV = -ThetaV;
+    Theta = Pi - Theta;
+  }
+
   //std::cout << "    dPhi = " <<   dPhi << "  (" <<   Phi << " <p|V> " <<   PhiV << ")" << std::endl;
   //std::cout << "  dTheta = " << dTheta << "  (" << Theta << " <p|V> " << ThetaV << ")" << std::endl;
 
@@ -1278,3 +1285,8 @@ void CosmicMuonGenerator::setAcptAllMu(bool AllMu) {
   if (NotInitialized)
     AcptAllMu = AllMu;
 }
+
+void CosmicMuonGenerator::setSurfaceDepth(double SurfaceDepthToCMS) {
+  if (NotInitialized)
+    SurfaceDepth = SurfaceDepthToCMS;
+}